var('x,y,t')
xp = x*(1-x-y)
yp = y*(3/4-y-x/2)
a = -1.5
b = 2
c = -1.5
d = 1.5
G = plot_vector_field([xp,yp],(x,a,b),(y,c,d))
# plot x-nullclines
G += parametric_plot((0,t),(t,c,d),color='red',thickness=5)
G += plot(1-x,(x,a,b),color='red',thickness=5,ymin=c,ymax=d)
# plot y=nullclines
G += parametric_plot((t,0),(t,a,b),color='purple',thickness=5)
G += plot(3/4-x/2,(x,a,b),color='purple',thickness=5,ymin=c,ymax=d)
# find equlibria (where red and purple cross)
equs = solve([xp,yp],(x,y),soln_dict=true)
show("Equilibria locate at the following:")
for equ in equs:
show(equ[0], ", ", equ[1])
X = function('X')(t)
Y = function('Y')(t)
de1 = diff(X,t) - xp(x=X,y=Y) == 0
de2 = diff(Y,t) - yp(x=X,y=Y) == 0
npts = 50
@interact
def _(x0 = slider(a,b,0.1,1),y0 = slider(c,d,0.1,1)):
soln = desolve_system_rk4([xp,yp], [x,y], ics=[0,x0,y0],
ivar=t, end_points=npts)
H = Graphics()
for k in range(npts):
H += points((soln[k][1],soln[k][2]))
(G+H).show(figsize=(10,10))
Click to the left again to hide and once more to show the dynamic interactive window |
|||||||||||||||||||||
plot_vector_field([xp,yp],(x,0.4,0.6),(y,0.4,0.6)).show()
|
var('x,y,t')
xp = x*(1-x-y)
yp = y*(3/4-y-x/2)
A1 = jacobian(xp,(x,y))
A2 = jacobian(yp,(x,y))
show(A1)
show(A2)
Jacob = matrix([[A1[0][0],A1[0][1]],[A2[0][0],A2[0][1]]])
show(Jacob)
for equ in equs:
A = Jacob(x=equ[0].rhs(),y=equ[1].rhs())
show(equ," has matrix ",A," with eigenvalues ",A.eigenvalues())
|
|
var('x,y,t')
xp = x*(1-x)-x*y
yp = x*y-y/2
a = -0.5
b = 1.3
c = -0.1
d = 1.1
G = plot_vector_field([xp,yp],(x,a,b),(y,c,d))
# plot x-nullclines
# G += parametric_plot((0,t),(t,c,d),color='red',thickness=5)
# G += plot(1-x,(x,a,b),color='red',thickness=5,ymin=c,ymax=d)
# plot y=nullclines
# G += parametric_plot((t,0),(t,a,b),color='purple',thickness=5)
# G += plot(3/4-x/2,(x,a,b),color='purple',thickness=5,ymin=c,ymax=d)
# find equlibria (where red and purple cross)
equs = solve([xp,yp],(x,y),soln_dict=true)
show("Equilibria locate at the following:")
for equ in equs:
show(equ[0], ", ", equ[1])
X = function('X')(t)
Y = function('Y')(t)
de1 = diff(X,t) - xp(x=X,y=Y) == 0
de2 = diff(Y,t) - yp(x=X,y=Y) == 0
npts = 200
@interact
def _(x0 = slider(a,b,0.1,1),y0 = slider(c,d,0.1,1)):
soln = desolve_system_rk4([xp,yp], [x,y], ics=[0,x0,y0],
ivar=t, end_points=npts)
H = Graphics()
for k in range(npts):
H += points((soln[k][1],soln[k][2]))
(G+H).show(figsize=(10,10))
Click to the left again to hide and once more to show the dynamic interactive window |
|||||||||||||||||||||
|
|