# euler.py - de methode van Euler, en de verbeterde Euler (Heun) als voorsmaakje van hoofdstuk 36
import numpy as np


def euler(f, y0, t0, t1, h):
    """Lost dy/dt = f(y, t) op van t0 tot t1 met stap h; geeft de tijden en de waarden terug."""
    n = round((t1 - t0) / h)
    t = t0 + h * np.arange(n + 1)
    y = np.empty(n + 1)
    y[0] = y0
    for k in range(n):
        y[k + 1] = y[k] + h * f(y[k], t[k])                 # een stap langs de raaklijn
    return t, y


def heun(f, y0, t0, t1, h):
    """Gemiddelde van de helling aan het begin en (geschat) aan het eind van elke stap."""
    n = round((t1 - t0) / h)
    t = t0 + h * np.arange(n + 1)
    y = np.empty(n + 1)
    y[0] = y0
    for k in range(n):
        k1 = f(y[k], t[k])
        k2 = f(y[k] + h * k1, t[k] + h)
        y[k + 1] = y[k] + h / 2 * (k1 + k2)
    return t, y


def station(T, t):                     # het station na een stroomuitval: dT/dt in graden per uur
    return -(T + 8) / 3


exact = -8 + 28 * np.exp(-1)          # T(3 h) exact
print("    h     Euler T(3)    fout    fout/h   Heun: fout   fout/h²")
for h in (1.0, 0.5, 0.25, 0.1, 0.01):
    fe = euler(station, 20, 0, 3, h)[1][-1] - exact
    fh = heun(station, 20, 0, 3, h)[1][-1] - exact
    print(f"{h:6.2f} {exact + fe:12.4f} {fe:8.4f} {fe / h:8.3f} {fh:12.2e} {fh / h**2:8.3f}")