import math
import numpy as np

# Populaties x,y zijn aantallen in een denkbeeldige proef; niet Habitat-data.
r=.60; a=.030; eta=.40; d=.30; cap=100.
xs=d/(eta*a); ys=r/a; ys_cap=r/a*(1-xs/cap)
def f(z,logistisch=False):
    x,y=z
    groei=r*x*(1-x/cap) if logistisch else r*x
    return np.array([groei-a*x*y,eta*a*x*y-d*y])
def integreer(h,eind=40.,methode="Euler",logistisch=False,z0=(40.,10.)):
    n=round(eind/h)
    if not math.isclose(n*h,eind): raise ValueError("niet gehele stapreeks")
    z=np.array(z0,float); zs=[z.copy()]
    for _ in range(n):
        if methode=="Euler": z=z+h*f(z,logistisch)
        elif methode=="RK4":
            k1=f(z,logistisch); k2=f(z+h*k1/2,logistisch)
            k3=f(z+h*k2/2,logistisch); k4=f(z+h*k3,logistisch)
            z=z+h*(k1+2*k2+2*k3+k4)/6
        else: raise ValueError("methode")
        if not np.all(np.isfinite(z)) or np.min(z)<=0:
            raise ValueError("negatieve/niet-eindige populatie; stap verkleinen")
        zs.append(z.copy())
    return np.arange(n+1)*h,np.array(zs)
def invariant(z):
    x,y=np.asarray(z).T
    return eta*a*(x-xs-xs*np.log(x/xs))+a*(y-ys-ys*np.log(y/ys))
J=np.array([[0.,-a*xs],[eta*a*ys,0.]])
Jcap=np.array([[-r*xs/cap,-a*xs],[eta*a*ys_cap,0.]])
print("LV evenwicht",xs,ys,"eigenwaarden",np.linalg.eigvals(J))
print("kleine periode",2*math.pi/math.sqrt(r*d),"dag")
print("logistische prooi",xs,ys_cap,"eigenwaarden",np.linalg.eigvals(Jcap))
ref=integreer(.0025,methode="RK4")
iv0=float(invariant(ref[1][0])); print("referentie invariantdrift",float(np.max(np.abs(invariant(ref[1])-iv0))))
rows=[]
for h in [.2,.1,.025]:
    t,z=integreer(h); drift=float(invariant(z[-1])-iv0)
    err=float(np.max(np.abs(z[-1]-ref[1][-1])))
    rows.append((h,z[-1,0],z[-1,1],drift,err))
    print("Euler h,x40,y40,Ddrift,max eindfout",*[round(v,6) for v in rows[-1]])
print("RK4 x40,y40",ref[1][-1])
try: integreer(.5)
except ValueError as fout: print("Euler h=0.5 geweigerd:",fout)
else: raise AssertionError("ongeldige grove proef niet geweigerd")
capref=integreer(.01,eind=80,methode="RK4",logistisch=True)
print("logistisch dag80",capref[1][-1])
assert abs(xs-25)<1e-12 and abs(ys-20)<1e-12
assert np.max(np.abs(invariant(ref[1])-iv0))<1e-8