352 - Topic 20 - Nonlinear systems and the Phase Plane

20 days ago by Professor352

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)) 
       
x0 
y0 

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)) 
       
x0 
y0 

Click to the left again to hide and once more to show the dynamic interactive window