import numpy as np
from scipy.optimize import least_squares

def bisectie(f, a, b, tol=1e-8):
    fa, fb = f(a), f(b)
    if fa == 0: return a
    if fb == 0: return b
    if fa * fb > 0: raise ValueError("geen tekenwisseling")
    while (b-a)/2 > tol:
        m = (a+b)/2
        fm = f(m)
        if fm == 0: return m
        if fa*fm < 0: b = m
        else: a, fa = m, fm
    return (a+b)/2

def rk4(f, y, h):
    k1 = f(y)
    k2 = f(y+h*k1/2)
    k3 = f(y+h*k2/2)
    k4 = f(y+h*k3)
    return y+h*(k1+2*k2+2*k3+k4)/6

# Twee knopen, dezelfde rekenwaarden als i31/e17. Tijd in seconde.
C = np.array([80000.0, 1100000.0])
t = np.arange(0, 8*3600+1, 600)
def model(logp):
    ua, koppeling = np.exp(logp)  # positieve W/K, schaalbaar fitten
    y = np.array([18.0,18.0])
    waarden = [y.copy()]
    for links,rechts in zip(t,t[1:]):
        vermogen = 800.0 if links < 4*3600 else 0.0
        def f(z):
            lucht, buffer = z
            q = koppeling*(lucht-buffer)
            return np.array([vermogen-ua*(lucht-5)-q,q])/C
        for _ in range(10): y = rk4(f,y,60.0)
        waarden.append(y.copy())
    return np.array(waarden)

rng = np.random.default_rng(36)
waar = model(np.log([20.2,30.0]))
metingen = waar + rng.normal(0,0.05,waar.shape)
fit = least_squares(lambda p: (model(p)-metingen).ravel()/0.05,
                    np.log([18,26]), bounds=(np.log([5,5]),np.log([80,100])))
ua, koppeling = np.exp(fit.x)
print("UA, H", *(f"{v:.3f}" for v in (ua,koppeling)), "W/K")
print("residu rms", f"{np.sqrt(np.mean((model(fit.x)-metingen)**2)):.4f}", "K")
print("nulpunt x^2-2", f"{bisectie(lambda x:x*x-2,1,2):.8f}")
for h in (1800,900,450):
    euler = 20.0
    r = 20.0
    for _ in range(int(10800/h)):
        euler += h*(-euler/10800)
        r = rk4(lambda z:-z/10800,r,h)
    exact = 20*np.exp(-1)
    print("stap", h, "Eulerfout", f"{abs(euler-exact):.6f}",
          "RK4-fout", f"{abs(r-exact):.6f}")
assert abs(ua-20.2)<0.2 and abs(koppeling-30)<0.5