CoCalc -- Lab09_slides.ipynb

Hosted by CoCalc

Download

Kernel: SageMath (stable)

Finding and classifying equilibrium points

Steps:

Find equilibrium point
Get the Jacobian
Test the eigenvalues

Solve Equilibrium Points for Continuous Time System

Example: Shark Tuna model:

$S' = 0.01ST-0.2S$

$T' = 0.05T-0.01ST$

In [3]:

_=var('S,T')
S_prime = 0.01*S*T-0.2*S
T_prime = 0.05*T-0.01*S*T

sol = solve([S_prime==0,T_prime==0],(S,T))
show(sol)

Out[3]:

How to retrieve the result?

In [4]:

sol

Out[4]:

[[S == 0, T == 0], [S == 5, T == 20]]

In [5]:

sol[0]

Out[5]:

[S == 0, T == 0]

In [7]:

sol[0][0]
sol[0][1]

Out[7]:

T == 0

In [9]:

sol[0][1].rhs()
sol[0][1].rhs()

Out[9]:

0

In [16]:

t = srange(0,100,0.1)
sol = desolve_odeint([S_prime,T_prime],ics =[4,18],times=t,dvars=[S,T])

fig1 = plot_vector_field((S_prime,T_prime),(S,3,7),(T,12,25),title='S=5,T=20')
fig2= point((5,20),color='red')
fig3 = list_plot(zip(sol[:,0],sol[:,1]),axes_labels=['S','R'])
show(fig1+fig2+fig3)

j = jacobian([S_prime,T_prime],[S,T])
j2 = j.subs({S:5,T:20})
show(j2.eigenvalues())

Out[16]:

In [17]:

t = srange(0,100,0.1)
sol = desolve_odeint([S_prime,T_prime],ics =[1,1],times=t,dvars=[S,T])

fig1 = plot_vector_field((S_prime,T_prime),(S,0,70),(T,0,100),title='S=0,T=0')
fig2= point((0,0),color='red')
fig3 = list_plot(zip(sol[:,0],sol[:,1]),axes_labels=['S','R'])
show(fig1+fig2+fig3)

j = jacobian([S_prime,T_prime],[S,T])
j1 = j.subs({S:0,T:0})
show(j1.eigenvalues())

Out[17]:

In [98]:

Out[98]:

Solve Does not Always Work

Model we are going to use today:

$H' = \frac{1}{1+G^n}-0.2H$

$P'=H-0.2P$

$G'=P-0.2G$

In [11]:

_ = var('H,P,G')
n=12
H_prime=1/(1+G^n)-0.2*H
P_prime=H-0.2*P
G_prime=P-0.2*G

sol = solve([H_prime,P_prime,G_prime],(H,P,G))
show(sol)

Out[11]:

None of them are real number!!! But there is at least one root!!!

What should we do when solve does not work?

Use: find_root!!!

find_root Can find "One" equilibrium point

For a "one" variable equation

Example: $x^7+cos(x)+2x-300$

In [90]:

_=var('x')

In [91]:

plot(x^7+cos(x)+2*x-300,(x,0,3.5))

Out[91]:

There is definitely a root between 2 and 3

find_root(fun,x_min,x_max,tolerlance)

In [100]:

find_root(x^7+cos(x)+2*x-300,2,3,rtol=0.001)

Out[100]:

2.2545915040920135

But there are 3 variables, not one!

$H' = \frac{1}{1+G^n}-0.2H$

$P'=H-0.2P$

$G'=P-0.2G$

However, if we only need to find equilibrium point, we only need to solve one variable!!!

$0 = \frac{1}{1+G^n}-0.2H$

$0=H-0.2P\implies 5H = P$

$0=P-0.2G\implies 5P=G$

$\implies 25H = G$

Finding equilibrium point is just solving a function of H!!!!

In [0]: