# voorbeeld chemische oscillaties
from sage.calculus.desolvers import desolve_odeint
from sage.calculus.functions import jacobian
x,y=var('x,y')
a=10
#b=3/5*a-25/a+0.4 # evenwichtspunt is stabiel brandpunt
b=3/5*a-25/a # Hopf-bifurcatie: evenwichtspunt (a/5, 1+(a/5)^2) wordt centrum en er onstaat een limietcyclus
#b=3/5*a-25/a-0.4 # evenwichtspunt is instabiel brandpunt en er is een limietcyclus
f=[a-x-4*x*y/(1+x^2),b*x*(1-y/(1+x^2))]
jac=jacobian(f, (x,y))
p1=plot_vector_field(f,(x,0,4),(y,0,10))
sol=desolve_odeint(f,[2,0.5],srange(0,50,0.1),[x,y])
p2=line(zip(sol[:,0],sol[:,1]))
sol=desolve_odeint(f,[2,1],srange(0,50,0.1),[x,y])
p3=line(zip(sol[:,0],sol[:,1]))
sol=desolve_odeint(f,[2,1.5],srange(0,50,0.1),[x,y])
p4=line(zip(sol[:,0],sol[:,1]))
sol=desolve_odeint(f,[2,2],srange(0,50,0.1),[x,y])
p5=line(zip(sol[:,0],sol[:,1]))
sol=desolve_odeint(f,[2,2.5],srange(0,50,0.1),[x,y])
p6=line(zip(sol[:,0],sol[:,1]))
sol=desolve_odeint(f,[2,3],srange(0,50,0.1),[x,y])
p7=line(zip(sol[:,0],sol[:,1]))
sol=desolve_odeint(f,[2,5.3],srange(0,5000,0.1),[x,y])
p9=line(zip(sol[:,0],sol[:,1]),color='green')
sol=desolve_odeint(f,[2,5.1],srange(0,5000,0.1),[x,y])
p8=line(zip(sol[:,0],sol[:,1]),color='red')
show(p1+p2+p3+p4+p5+p6+p7+p9+p8)
A0=jac.subs(x==a/5,y==1+(a/5)^2)
pretty_print(A0)
pretty_print(A0.eigenvalues())
pretty_print(A0.trace())
pretty_print(A0.determinant())
