%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import math

def euler(func, initial, deltas):
    ys = [initial]
    for i in range(1, len(xs)):
        cur = np.array(ys[-1][:])
        dt = xs[i] - xs[i - 1]
        cur = cur + np.array(func(cur, xs[i])) * dt
        ys.append(cur)
    return np.array(ys)



from scipy.integrate import odeint

def growth(u, t):
    return u[0] * 0.1

def f(k, r, t):
    return k * math.e ** (r * t)

init = [4]
xs = np.linspace(0, 20, 20)

ys = euler(growth, init, xs)
ybs = odeint(growth, init, xs)

plt.plot(xs, ys[:, 0])
plt.plot(xs, ybs[:, 0])
plt.plot(xs, f(4, 0.1, xs))

old = 1
for i in range(20):
    dt = 10.0 / 2 ** i
    end = 20
    xs = np.linspace(0, end, int(end / dt))
    ys = euler(growth, init, xs)
    newerr = abs(ys[-1][0] - f(4, 0.1, end))
    print("%f \t%f\t%f\t%f\t%f\t" % \
          (dt, ys[-1][0], f(4, 0.1, end), newerr, newerr / old))
    old = newerr

def rk2(func, initial, deltas):
    ys = [initial]
    for i in range(1, len(xs)):
        old = np.array(ys[-1][:])
        dt = xs[i] - xs[i - 1]
        d1 = np.array(func(old, xs[i]))
        cur = old + np.array(func(old, xs[i])) * dt
        d2 = np.array(func(cur, xs[i] + dt))
        better = old + ((d1 + d2) / 2.0) * dt
        ys.append(better)
    return np.array(ys)



init = [4]
xs = np.linspace(0, 20, 5)

ys = euler(growth, init, xs)
ybs = rk2(growth, init, xs)
ycs = odeint(growth, init, xs)

plt.plot(xs, ys[:, 0])
plt.plot(xs, ybs[:, 0])
plt.plot(xs, ycs[:, 0])

old = 1
for i in range(20):
    dt = 10.0 / 2 ** i
    end = 20
    xs = np.linspace(0, end, int(end / dt))
    ys = rk2(growth, init, xs)
    newerr = abs(ys[-1][0] - f(4, 0.1, end))
    print("%f \t%f\t%f\t%f\t%f\t" % \
          (dt, ys[-1][0], f(4, 0.1, end), newerr, newerr / old))
    old = newerr

def rk4(func, initial, deltas):
    ys = [initial]
    for i in range(1, len(xs)):
        old = np.array(ys[-1][:])
        dt = xs[i] - xs[i - 1]
        d1 = np.array(func(old, xs[i])) * dt
        d2 = np.array(func(old + 0.5 * d1, xs[i] + 0.5 * dt)) * dt
        d3 = np.array(func(old + 0.5 * d2, xs[i] + 0.5 * dt)) * dt
        d4 = np.array(func(old + d3, xs[i] + dt)) * dt
        better = old + ((d1 + 2 * d2 + 2 * d3 + d4) / 6.0)
        ys.append(better)
    return np.array(ys)



old = 1
for i in range(20):
    dt = 10.0 / 2 ** i
    end = 20
    xs = np.linspace(0, end, int(end / dt))
    ys = rk4(growth, init, xs)
    newerr = abs(ys[-1][0] - f(4, 0.1, end))
    print("%f \t%f\t%f\t%f\t%f\t" % \
          (dt, ys[-1][0], f(4, 0.1, end), newerr, newerr / old))
    old = newerr