from __future__ import print_function
from numpy import *
#import matplotlib.pyplot as plt

def euler(f, x0, t0, h, N):
    t = zeros(N+1)
    x = zeros(N+1)
    t[0] = t0
    x[0] = x0
    for k in range(N):
        x[k+1] = x[k] + h*f(t[k], x[k])
        t[k+1] = t[k] + h
    
    return x, t

def euler_mid(f, x0, t0, h, N):
    t = zeros(N+1)
    x = zeros(N+1)
    t[0] = t0
    x[0] = x0
    for k in range(N):
        x_mid = x[k] + (h/2.0)*f(t[k], x[k])
        x[k+1] = x[k] + h*f(t[k]+h/2.0, x_mid)
        t[k+1] = t[k] + h
    
    return x, t

def runge_kutta4(f, x0, t0, h, N):
    t = zeros(N+1)
    x = zeros(N+1)
    t[0] = t0
    x[0] = x0
    for k in range(N):
        k0 = f(t[k], x[k])
        k1 = f(t[k]+0.5*h, x[k] + 0.5*h*k0)
        k2 = f(t[k]+0.5*h, x[k] + 0.5*h*k1)
        k3 = f(t[k]+h, x[k] + h*k2)
        x[k+1] = x[k] + (h/6.0)*(k0 + 2*k1 + 2*k2 + k3)
        t[k+1] = t[k] + h
    
    return x,t

exact = lambda t: exp(t) - t - 1
f = lambda t, x: x + t

a = 0.0;
b = 1.0;
x0 = 0.0;

for i in range(1,6):
    N = 10**(i)
    h = (b-a)/float(N)
    x_euler, t1     = euler(f, x0, a, h, N)
    x_euler_mid, t2 = euler_mid(f, x0, a, h, N)
    x_rk4, t3 = runge_kutta4(f, x0, a, h, N)
    
    abserr_euler     = abs(x_euler[-1] - exact(b))
    abserr_euler_mid = abs(x_euler_mid[-1] - exact(b))
    abserr_rk4 = abs(x_rk4[-1] - exact(b))
    print("h: %6e: err euler: %12e, err euler mid: %12e" % \
          (h, abserr_euler, abserr_euler_mid))