Contact
CoCalc Logo Icon
StoreFeaturesDocsShareSupport News AboutSign UpSign In
| Download
Views: 111
#1 var("M","H") Mprime=0.36*M-0.02*M^2-0.03*M*H Hprime=0.01*M*H-0.1*H-0.005*H^2 solve([Mprime,Hprime],[M,H])
(M, H) [[M == 0, H == 0], [M == 18, H == 0], [M == 0, H == -20], [M == 12, H == 4]]
#2 j=jacobian([Mprime,Hprime],[M,H]) show(j)
(0.0300000000000000H0.0400000000000000M+0.3600000000000000.0300000000000000M0.0100000000000000H0.0100000000000000H+0.0100000000000000M0.100000000000000)\displaystyle \left(\begin{array}{rr} -0.0300000000000000 \, H - 0.0400000000000000 \, M + 0.360000000000000 & -0.0300000000000000 \, M \\ 0.0100000000000000 \, H & -0.0100000000000000 \, H + 0.0100000000000000 \, M - 0.100000000000000 \end{array}\right)
#3 j.subs({M:0,H:0})
[ 0.360000000000000 0] [ 0 -0.100000000000000]
j.subs({M:18,H:0})
[-0.360000000000000 -0.540000000000000] [ 0 0.0800000000000000]
j.subs({M:0,H:-20})
[ 0.960000000000000 0] [-0.200000000000000 0.100000000000000]
j.subs({M:12,H:4})
[ -0.240000000000000 -0.360000000000000] [ 0.0400000000000000 -0.0200000000000000]
#4 var("M","H") t=srange(0,100,0.1) sol=desolve_odeint([0.36*M-0.02*M^2-0.03*M*H,0.01*M*H-0.1*H-0.005*H^2], ics=[1,1], dvars=[M,H],times=t) plot_vector_field([0.36*M-0.02*M^2-0.03*M*H,0.01*M*H-0.1*H-0.005*H^2], (M,0,18), (H,0,5))+point([1,1], size=30)+point([0,0], size=30, color="red")+list_plot(sol, color="orange",plotjoined=True, thickness=2)+point([18,0], size=30, color="purple")+point([12,4],size=30,color="green")
(M, H)
#5 var("N","P") r1=0.4 r2=0.03 j=70 w=300 d=1000 k=3000 Nprime=r1*N*(1-(N/k))-w*((N)/(d+N))*P Pprime=r2*P*(1-((j*P)/(N))) solve([Nprime,Pprime],[N,P])
(N, P) [[N == 0, P == 0], [N == -1000, P == 0], [N == 3000, P == 0], [N == -500/7*sqrt(45109) - 105500/7, P == -50/49*sqrt(45109) - 10550/49], [N == 500/7*sqrt(45109) - 105500/7, P == 50/49*sqrt(45109) - 10550/49]]
#6 p=jacobian([Nprime,Pprime],[N,P]) show(p)
(0.000266666666666667N300PN+1000+300NP(N+1000)2+0.400000000000000300NN+10002.10000000000000P2N24.20000000000000PN+0.0300000000000000)\displaystyle \left(\begin{array}{rr} -0.000266666666666667 \, N - \frac{300 \, P}{N + 1000} + \frac{300 \, N P}{{\left(N + 1000\right)}^{2}} + 0.400000000000000 & -\frac{300 \, N}{N + 1000} \\ \frac{2.10000000000000 \, P^{2}}{N^{2}} & -\frac{4.20000000000000 \, P}{N} + 0.0300000000000000 \end{array}\right)
r=p.subs({N:500/7*sqrt(45109) - 105500/7,P:50/49*sqrt(45109) - 10550/49}) show(r)
(0.01904761904761904510930(45109211)7(45109197)+30(45109211)27(45109197)2+4.41904761904762300(45109211)451091970.0004285714285714290.0300000000000000)\displaystyle \left(\begin{array}{rr} -0.0190476190476190 \, \sqrt{45109} - \frac{30 \, {\left(\sqrt{45109} - 211\right)}}{7 \, {\left(\sqrt{45109} - 197\right)}} + \frac{30 \, {\left(\sqrt{45109} - 211\right)}^{2}}{7 \, {\left(\sqrt{45109} - 197\right)}^{2}} + 4.41904761904762 & -\frac{300 \, {\left(\sqrt{45109} - 211\right)}}{\sqrt{45109} - 197} \\ 0.000428571428571429 & -0.0300000000000000 \end{array}\right)
r.eigenvalues()
[1/1400*(1185703*sqrt(45109) - sqrt(-601256621013526*sqrt(45109) + 127700168580419822) - 251829341)/(197*sqrt(45109) - 41959), 1/1400*(1185703*sqrt(45109) + sqrt(-601256621013526*sqrt(45109) + 127700168580419822) - 251829341)/(197*sqrt(45109) - 41959)]
#It is a saddle point
#7 t=srange(0,100,0.1) sol=desolve_odeint([r1*N*(1-(N/k))-w*((N)/(d+N))*P,r2*P*(1-((j*P)/(N)))], ics=[1,1], dvars=[N,P],times=t) plot_vector_field([r1*N*(1-(N/k))-w*((N)/(d+N))*P,r2*P*(1-((j*P)/(N)))], (N,0,3000), (P,0,5))+ list_plot(sol,plotjoined=True)
#8 var("N,P") r1=0.4 r2=0.03 j=70 w=12 d=1000 k=3000 Nprime=r1*N*(1-(N/k))-w*((N)/(d+N))*P Pprime=r2*P*(1-((j*P)/(N))) solve([Nprime,Pprime],[N,P])
(N, P) [[N == 0, P == 0], [N == -1000, P == 0], [N == 3000, P == 0], [N == -500/7*sqrt(613) + 2500/7, P == -50/49*sqrt(613) + 250/49], [N == 500/7*sqrt(613) + 2500/7, P == 50/49*sqrt(613) + 250/49]]
"A change in w resulted in a change in stability"
'A change in w resulted in a change in stability'
boo=jacobian([Nprime,Pprime],[N,P]) show(boo) boo1=boo.subs(N=500/7*sqrt(613) + 2500/7,P=50/49*sqrt(613) + 250/49) boo1.eigenvalues()
(0.000266666666666667N12PN+1000+12NP(N+1000)2+0.40000000000000012NN+10002.10000000000000P2N24.20000000000000PN+0.0300000000000000)\displaystyle \left(\begin{array}{rr} -0.000266666666666667 \, N - \frac{12 \, P}{N + 1000} + \frac{12 \, N P}{{\left(N + 1000\right)}^{2}} + 0.400000000000000 & -\frac{12 \, N}{N + 1000} \\ \frac{2.10000000000000 \, P^{2}}{N^{2}} & -\frac{4.20000000000000 \, P}{N} + 0.0300000000000000 \end{array}\right)
[-1/1400*(3679*sqrt(613) + sqrt(134733626*sqrt(613) + 4000114862) + 65827)/(19*sqrt(613) + 487), -1/1400*(3679*sqrt(613) - sqrt(134733626*sqrt(613) + 4000114862) + 65827)/(19*sqrt(613) + 487)]
#9 sol=desolve_odeint([r1*N*(1-(N/k))-12*((N)/(d+N))*P,r2*P*(1-((j*P)/(N)))], ics=[1,1], dvars=[N,P],times=t) plot_vector_field([r1*N*(1-(N/k))-3000*((N)/(d+N))*P,r2*P*(1-((j*P)/(N)))], (N,0,3000), (P,0,5))+ list_plot(sol,plotjoined=True)
"The arrows are pointing the opposite way. That change in stability is the result of a Hopf Birurcation"
'The arrows are pointing the opposite way. That change in stability is the result of a Hopf Birurcation'
#10 "The simplified equation is 0.0008*G*(1+0.008*G^12)-1=0"
'The simplified equation is 0.0008*G*(1+0.008*G^12)-1=0'
#11 var("G") find_root(0.0008*G*(1+0.008*G^12)-1,0,10)
G 2.508750670737234
#P=0.2G 0.2*2.508750670737234
0.501750134147447
#H=0.2*P 0.2*0.501750134147447
0.100350026829489
#12 "It is an unstabe eq point"
'It is an unstabe eq point'
#13 H_Prime=(1/(1+G^n))-0.2*H P_prime=H-0.2*P G_prime=P-0.2*G t=srange(0,400,0.1) sol=desolve_odeint([H_prime,P_prime,G_prime], ics=[1,2,3], dvars=[H,P,G], times=t) list_plot(sol, plotjoined=True) #I do not know why it came out as an error
Error in lines 1-1 Traceback (most recent call last): File "/cocalc/lib/python2.7/site-packages/smc_sagews/sage_server.py", line 1188, in execute flags=compile_flags) in namespace, locals File "", line 1, in <module> TypeError: unsupported operand type(s) for ** or pow(): 'sage.symbolic.expression.Expression' and 'function'
#14 var("G") find_root(0.0008*G*(1+0.008*G^6)-1,0,10)
G 5.516960270495793
#P=0.2G 0.2*5.516960270495793
1.10339205409916
#H=0.2*P 0.2*1.10339205409916
0.220678410819832
"It is a stable spiral eq point"
'It is a stable spiral eq point'
#15 @interact def model(n=(2.0,12.0)): var("G,H,P") Hprime=(1/(1+G^n))-0.2*H Pprime=H-0.2*P Gprime=P-0.2*G t=srange(0,1000) sol=desolve_odeint([Hprime,Pprime,Gprime], ics=[1,1,1], dvars=[H,P,G], times=t) list_plot(sol,plotjoined=True)
Interact: please open in CoCalc
"The birfurcation ocuurs when n equals about 7.5"
'The birfurcation ocuurs when n equals about 7.5'
#16 "It is much easier to find when a Hopf birurcation occurs on Sage because it is less time consuming than doing it by hand and it does not leave much room for mistakes or inconsistencies"
'It is much easier to find when a Hopf birurcation occurs on Sage because it is less time consuming than doing it by hand and it does not leave much room for mistakes or inconsistencies'
#17 "The eq point should undergo a change in stability in order for a Hopf Bifurcation to occcur. Detecting it numerically, although very time consuming, can be accomplished by seeing how a change in parameter values can affect the (negative or positive) signs of the eigenvalues. If an eigenvalue changes from being negative to positive or vice versa, a Hopf bifucation has occured"
'The eq point should undergo a change in stability in order for a Hopf Bifurcation to occcur. Detecting it numerically, although very time consuming, can be accomplished by seeing how a change in parameter values can affect the (negative or positive) signs of the eigenvalues. If an eigenvalue changes from being negative to positive or vice versa, a Hopf bifucation has occured'
#18 @interact def function (n=(2,12,0.01)): var("H,P,G") Hprime= Pprime= Gprime= h1=find_root((1/(1+(25*H)^n))-0.2*H,0,10) p1=5*0.0795469 g1=5*p1 j=jacobian([Hprime,Pprime,Gprime],[H,P,G]) M=matrix(RDF,M) M.eigenvalues()
Error in lines 1-12 Traceback (most recent call last): File "/cocalc/lib/python2.7/site-packages/smc_sagews/sage_server.py", line 1188, in execute flags=compile_flags) in namespace, locals File "<string>", line 4 Hprime= ^ SyntaxError: invalid syntax