An Autonomous Nonlinear System

Here is an example of a nonlinear, autonomous 2X2 systems of ordinary differential equations.

This notebook was created by T. Shelton, Southwestern University to supplement a handout for  Differential Equations II, updated for Mathematica 11.1.

reference  Differential Equations: Theory, Technique, and Practice, Simmons & Krantz, McGraw-Hill Co, Inc. 2007

dx/dt=f[x,y]
dy/dt=g[x,y]

Nullclines, Direction Field, Equilibria

Clear[f,g]
f[x_,y_]:= x - x^2 - x y
g[x_,y_]:= - 1/4  y - 1/4  y^2

x-nullclines

Solve[f[x,y]==0]

Solve[f[x,y]==0,y]

Above, we see that one of the x-nullclines cannot be expressed for y as a function of x. 
Below, either method of solving yields the same y-nullclines.

y-nullclines

Solve[g[x,y]==0]

Solve[g[x,y]==0,y]

equilibria

Clear[CPRules,numCPs, CPs]
CPRules = Solve[{f[x,y]==0,g[x,y]==0},{x,y}]
numCPs = Length[CPRules]
CPs = Table[{x,y}/.CPRules[[i]],{i,1,numCPs}]

Thus, we have four isolated critical points or equilibria.

We choose a graphing window containing interesting places, and we create a plot of nullclines and equilibria.

Clear[xmin,xmax,ymin,ymax]
{xmin,xmax,ymin,ymax} = {-2.0,4.0,-4.0,4.0};
Clear[npx1,npx2,npy1, npy2,np,pts]
npx1 =Plot[1-x,{x, xmin, xmax},PlotStyle->{Red,Dashed,Thick}];
npx2 =ParametricPlot[{0,t},{t, ymin, ymax},PlotStyle->{Red,Dashed,Thick}];
npy1 = Plot[-1,{x, xmin, xmax},PlotStyle->{Purple,DotDashed,Thick}];
npy2 = Plot[0,{x, xmin, xmax},PlotStyle->{Purple,DotDashed,Thick}];
np = Show[{npx1,npx2,npy1,npy2},PlotRange->{{xmin,xmax},{ymin,ymax}}];
pts = Graphics[{PointSize[Large],Point[CPs ]}];

When a nullcline is along an axis, the graph is clearer with a frame and no axes.

{Show[np,pts], Show[{np,pts}, Axes->False, Frame->True]}

the direction field

Clear[vPlot]
vPlot = VectorPlot[{f[x,y],g[x,y]},{x,xmin,xmax},{y,ymin,ymax}];

all together

Show[{np,pts,vPlot}]

The scaling of the vectors above is more accurate, but below we can actually see vectors even when the magnitude is close to 0.  Each can be beneficial.

Clear[vPlot2]
vPlot2 = VectorPlot[{f[x,y],g[x,y]},{x,xmin,xmax},{y,ymin,ymax},VectorScale->{Small,Small,None}];

Show[{np,pts,vPlot2    }, Axes->False, Frame->True]

We show a Mathematica StreamPlot instead of a vector field.

Show[{np,pts,StreamPlot[{f[x,y],g[x,y]},{x,xmin,xmax},{y,ymin,ymax}]   }, Axes->False, Frame->True]

We show two unsatisfactory graphs.  We really need to redo the vector plots if we change the graphing window.

Clear[xmin,xmax,ymin,ymax]
{xmin,xmax,ymin,ymax} = {-.5, 1.5, -.5, .8};{Show[{np,pts,vPlot    },PlotRange->{{xmin,xmax},{ymin,ymax}}, Axes->False, Frame->True],Show[{np,pts,vPlot2    },PlotRange->{{xmin,xmax},{ymin,ymax}}, Axes->False, Frame->True]}

Clear[xmin,xmax,ymin,ymax]
{xmin,xmax,ymin,ymax} = {-.5, 1.5, -.5, .8};
Show[{np,pts,
VectorPlot[{f[x,y],g[x,y]},{x,xmin,xmax},{y,ymin,ymax},VectorScale->{Small,Small,None}]
},PlotRange->{{xmin,xmax},{ymin,ymax}}, Axes->False, Frame->True]

Clear[xmin,xmax,ymin,ymax]
{xmin,xmax,ymin,ymax} = {1.5, 2.5, -1.5, 0.};
Show[{np,pts,
VectorPlot[{f[x,y],g[x,y]},{x,xmin,xmax},{y,ymin,ymax},VectorScale->{Small,Small,None}]
},PlotRange->{{xmin,xmax},{ymin,ymax}}, Axes->False, Frame->True]

Clear[xmin,xmax,ymin,ymax]
{xmin,xmax,ymin,ymax} = {-1.5, 2.5, -1.5,-.5};
Show[{np,pts,
VectorPlot[{f[x,y],g[x,y]},{x,xmin,xmax},{y,ymin,ymax},VectorScale->{Small,Small,None}]
},PlotRange->{{xmin,xmax},{ymin,ymax}}, Axes->False, Frame->True]

Clear[xmin,xmax,ymin,ymax]
{xmin,xmax,ymin,ymax} = {-0.5, 0.5, -1.2,-.8};
Show[{np,pts,
VectorPlot[{f[x,y],g[x,y]},{x,xmin,xmax},{y,ymin,ymax},VectorScale->{Small,Small,None}]
},PlotRange->{{xmin,xmax},{ymin,ymax}}, Axes->False, Frame->True]

Clear[xmin,xmax,ymin,ymax]
{xmin,xmax,ymin,ymax} = {1.5, 2.5, -1.2,-.8};
Show[{np,pts,
VectorPlot[{f[x,y],g[x,y]},{x,xmin,xmax},{y,ymin,ymax},VectorScale->{Small,Small,None}]
},PlotRange->{{xmin,xmax},{ymin,ymax}}, Axes->False, Frame->True]

Solution Technique:  Linearization and the Jacobian Eigensystem Method

dx/dt=f[x,y], dy/dt=g[x,y]  where f[x_, y_] := x-(x^2)-x y, g[x_, y_] := -(y/4)-((y^2)/4).
We compute the Jacobian.

Clear[J]
J[x_,y_] = {   {  D[f[x,y],x],  D[f[x,y],y] }, {  D[g[x,y],x],  D[g[x,y],y] }  }

CP 1

Clear[i,A, x0, y0]
i=1;
{x0,y0} = {x,y}/.CPRules[[i]]
A = J[x0,y0]

Clear[eval1,eval2,evec1,evec2]
{{eval1,eval2},{evec1,evec2}} = Eigensystem[A]

evals are -1 and -1/4.
Our eigenvalues are real, not equal, and have the same algebraic sign.  Therefore, this critical point is an improper node.

Clear[xmin,xmax,ymin,ymax]
{xmin,xmax,ymin,ymax} = {x0-.25,x0+.25, y0-.25, y0+.25}

Show[{VectorPlot[{f[x,y],g[x,y]},{x,xmin,xmax},{y,ymin,ymax}],
ListPlot[{{x0,y0}},PlotStyle->{PointSize[.02],Red}] }]

Show[{StreamPlot[{f[x,y],g[x,y]},{x,xmin,xmax},{y,ymin,ymax}],
VectorPlot[{f[x,y],g[x,y]},{x,xmin,xmax},{y,ymin,ymax},VectorStyle->{Black}],
ListPlot[{{x0,y0}},PlotStyle->{PointSize[.02],Red}] }]

Show[{Plot[1-x,{x, xmin, xmax},PlotStyle->{Red,Dashed,Thick}],ParametricPlot[{0,t},{t, ymin, ymax},PlotStyle->{Red,Dashed,Thick}], Plot[-1,{x, xmin, xmax},PlotStyle->{Purple,DotDashed,Thick}],Plot[0,{x, xmin, xmax},PlotStyle->{Purple,DotDashed,Thick}],

VectorPlot[{f[x,y],g[x,y]},{x,xmin,xmax},{y,ymin,ymax},VectorStyle->{Black}],
ListPlot[{{x0,y0}},PlotStyle->{PointSize[.02],Red}] }]

solution to the associated linearized system near the critical point

check for linear independence of the first solution, written as a row vector,...

Clear[x1]
x1=Expand[ (E^(eval1 * t))evec1]

... and the second solution, written as a row vector,...

Clear[x2]
x2=Expand[(E^(eval2*t))evec2]

here is the matrix whose determinant is the Wronskian of the two solutions

MatrixForm[Transpose[ {x1,x2}]]

Det[ Transpose[ {x1,x2}]]

Since this non-zero constant multiple of an exponential function is never zero, our two solutions are linearly independent.  

Here is our linearized solution near this critical point:

Clear[LocalLinSoln]
LocalLinSoln = Expand[c1 x1 + c2 x2+{x0,y0}]

Clear[TrajectoryTable]
TrajectoryTable = Table[LocalLinSoln/.{c1->a,c2->b},{a,-.2,.2,.2},{b,-.2,.2,.2}]

Show[{ParametricPlot[TrajectoryTable,{t,0,6}],
ListPlot[TrajectoryTable/.{t->0},PlotStyle->{PointSize[.02],Black}] ,
ListPlot[{{x0,y0}},PlotStyle->{PointSize[.02],Red}]  }]

Show[{ParametricPlot[TrajectoryTable,{t,0,6}],
ListPlot[TrajectoryTable/.{t->0},PlotStyle->{PointSize[.02],Black}] ,
VectorPlot[{f[x,y],g[x,y]},{x,xmin,xmax},{y,ymin,ymax}],
ListPlot[{{x0,y0}},PlotStyle->{PointSize[.02],Red}] }]

Change of variables

Clear[f1,g1,u,v]
f1[u_,v_] =Expand[ f[u+x0, v+y0]]
g1[u_,v_] = Expand[ g[u+x0, v+y0] ]

Solve[{f1[u,v]==0, g1[u,v]==0}]

Note that {u,v} = {0,0} is an isolated critical point of the transformed system.

Clear[J1]
J1[u_,v_] = {   {  D[f1[u,v],u],  D[f1[u,v],v] }, {  D[g1[u,v],u],  D[g1[u,v],v] }  }

Clear[A1]
A1 = J1[0,0]

{f1[u,v], g1[u,v]} - A1.{u,v}

Clear[r,s]
{r[u_,v_] , s[u_,v_] } = {f1[u,v], g1[u,v]} - A1.{u,v}

Eigensystem[A1]

Clear[umin,umax,vmin,vmax]
{umin,umax,vmin,vmax} = {xmin,xmax,ymin,ymax} + {x0,x0,y0,y0}

{  Show[{StreamPlot[{f[x,y],g[x,y]},{x,xmin,xmax},{y,ymin,ymax}],
VectorPlot[{f[x,y],g[x,y]},{x,xmin,xmax},{y,ymin,ymax},VectorStyle->{Black}],
ListPlot[{{x0,y0}},PlotStyle->{PointSize[.02],Red}] }],
Show[{StreamPlot[{f1[u,v],g1[u,v]},{u,umin,umax},{v,vmin,vmax}],
VectorPlot[{f1[u,v],g1[u,v]},{u,umin,umax},{v,vmin,vmax},VectorStyle->{Black}],
ListPlot[{{0,0}},PlotStyle->{PointSize[.02],Red}] }]    }

CP 2

Clear[i,A, x0, y0]
i=2;
{x0,y0} = {x,y}/.CPRules[[i]]
A = J[x0,y0]

Clear[eval1,eval2,evec1,evec2]
{{eval1,eval2},{evec1,evec2}} = Eigensystem[A]

evals are real, not equal, and have opposite algebraic signs.  Therefore, this critical point is a saddle point.

Clear[xmin,xmax,ymin,ymax]
{xmin,xmax,ymin,ymax} = {x0-.25,x0+.25, y0-.25, y0+.25}

Show[{VectorPlot[{f[x,y],g[x,y]},{x,xmin,xmax},{y,ymin,ymax}],
ListPlot[{{x0,y0}},PlotStyle->{PointSize[.02],Red}] }]

Show[{StreamPlot[{f[x,y],g[x,y]},{x,xmin,xmax},{y,ymin,ymax}],
VectorPlot[{f[x,y],g[x,y]},{x,xmin,xmax},{y,ymin,ymax},VectorStyle->{Black}],
ListPlot[{{x0,y0}},PlotStyle->{PointSize[.02],Red}] }]

Show[{Plot[1-x,{x, xmin, xmax},PlotStyle->{Red,Dashed,Thick}],ParametricPlot[{0,t},{t, ymin, ymax},PlotStyle->{Red,Dashed,Thick}], Plot[-1,{x, xmin, xmax},PlotStyle->{Purple,DotDashed,Thick}],Plot[0,{x, xmin, xmax},PlotStyle->{Purple,DotDashed,Thick}],

VectorPlot[{f[x,y],g[x,y]},{x,xmin,xmax},{y,ymin,ymax},VectorStyle->{Black}],
ListPlot[{{x0,y0}},PlotStyle->{PointSize[.02],Red}] }]

solution to the associated linearized system near the critical point

check for linear independence of the first solution, written as a row vector,...

Clear[x1]
x1=Expand[ (E^(eval1 * t))evec1]

... and the second solution, written as a row vector,...

Clear[x2]
x2=Expand[(E^(eval2*t))evec2]

here is the matrix whose determinant is the Wronskian of the two solutions

MatrixForm[Transpose[ {x1,x2}]]

Det[ Transpose[ {x1,x2}]]

Since this non-zero constant multiple of an exponential function is never zero, our two solutions are linearly independent.  

Here is our linearized solution near this critical point:

Clear[LocalLinSoln]
LocalLinSoln = Expand[c1 x1 + c2 x2+{x0,y0}]

Clear[TrajectoryTable]
TrajectoryTable = Table[LocalLinSoln/.{c1->a,c2->b},{a,-.2,.2,.2},{b,-.2,.2,.2}]

Show[{ParametricPlot[TrajectoryTable,{t,0,2}],
ListPlot[TrajectoryTable/.{t->0},PlotStyle->{PointSize[.02],Black}] ,
ListPlot[{{x0,y0}},PlotStyle->{PointSize[.02],Red}]  }]

Show[{ParametricPlot[TrajectoryTable,{t,0,2}],
ListPlot[TrajectoryTable/.{t->0},PlotStyle->{PointSize[.02],Black}] ,
VectorPlot[{f[x,y],g[x,y]},{x,xmin,xmax},{y,ymin,ymax}],
ListPlot[{{x0,y0}},PlotStyle->{PointSize[.02],Red}] }]

CP 3

Clear[i,A, x0, y0]
i=3;
{x0,y0} = {x,y}/.CPRules[[i]]
A = J[x0,y0]

Clear[eval1,eval2,evec1,evec2]
{{eval1,eval2},{evec1,evec2}} = Eigensystem[A]

evals are real, not equal, and have the same sign.  Therefore, this critical point is an unstable node.

Clear[xmin,xmax,ymin,ymax]
{xmin,xmax,ymin,ymax} = {x0-.25,x0+.25, y0-.25, y0+.25}

Show[{VectorPlot[{f[x,y],g[x,y]},{x,xmin,xmax},{y,ymin,ymax}],
ListPlot[{{x0,y0}},PlotStyle->{PointSize[.02],Red}] }]

Show[{StreamPlot[{f[x,y],g[x,y]},{x,xmin,xmax},{y,ymin,ymax}],
VectorPlot[{f[x,y],g[x,y]},{x,xmin,xmax},{y,ymin,ymax},VectorStyle->{Black}],
ListPlot[{{x0,y0}},PlotStyle->{PointSize[.02],Red}] }]

Show[{Plot[1-x,{x, xmin, xmax},PlotStyle->{Red,Dashed,Thick}],ParametricPlot[{0,t},{t, ymin, ymax},PlotStyle->{Red,Dashed,Thick}], Plot[-1,{x, xmin, xmax},PlotStyle->{Purple,DotDashed,Thick}],Plot[0,{x, xmin, xmax},PlotStyle->{Purple,DotDashed,Thick}],

VectorPlot[{f[x,y],g[x,y]},{x,xmin,xmax},{y,ymin,ymax},VectorStyle->{Black}],
ListPlot[{{x0,y0}},PlotStyle->{PointSize[.02],Red}] }, PlotRange->{{xmin,xmax},{ymin,ymax}}]

solution to the associated linearized system near the critical point

check for linear independence of the first solution, written as a row vector,...

Clear[x1]
x1=Expand[ (E^(eval1 * t))evec1]

... and the second solution, written as a row vector,...

Clear[x2]
x2=Expand[(E^(eval2*t))evec2]

here is the matrix whose determinant is the Wronskian of the two solutions

MatrixForm[Transpose[ {x1,x2}]]

Det[ Transpose[ {x1,x2}]]

Since this non-zero constant multiple of an exponential function is never zero, our two solutions are linearly independent.  

Here is our linearized solution near this critical point:

Clear[LocalLinSoln]
LocalLinSoln = Expand[c1 x1 + c2 x2+{x0,y0}]

Clear[TrajectoryTable]
TrajectoryTable = Table[LocalLinSoln/.{c1->a,c2->b},{a,-.2,.2,.2},{b,-.2,.2,.2}]

Show[{ParametricPlot[TrajectoryTable,{t,-.5,.5}],
ListPlot[TrajectoryTable/.{t->-.5},PlotStyle->{PointSize[.02],Black}] ,
ListPlot[{{x0,y0}},PlotStyle->{PointSize[.02],Red}]  }]

Show[{ParametricPlot[TrajectoryTable,{t,-.5,.5}],
ListPlot[TrajectoryTable/.{t->-.5},PlotStyle->{PointSize[.02],Black}] ,
VectorPlot[{f[x,y],g[x,y]},{x,xmin,xmax},{y,ymin,ymax}],
ListPlot[{{x0,y0}},PlotStyle->{PointSize[.02],Red}] }]

CP 4

Clear[i,A, x0, y0]
i=4;
{x0,y0} = {x,y}/.CPRules[[i]]
A = J[x0,y0]

Clear[eval1,eval2,evec1,evec2]
{{eval1,eval2},{evec1,evec2}} = Eigensystem[A]

evals are real, not equal, and have opposite signs.  Therefore, this critical point is a saddle point..

Clear[xmin,xmax,ymin,ymax]
{xmin,xmax,ymin,ymax} = {x0-.25,x0+.25, y0-.25, y0+.25}

Show[{VectorPlot[{f[x,y],g[x,y]},{x,xmin,xmax},{y,ymin,ymax}],
ListPlot[{{x0,y0}},PlotStyle->{PointSize[.02],Red}] }]

Show[{StreamPlot[{f[x,y],g[x,y]},{x,xmin,xmax},{y,ymin,ymax}],
VectorPlot[{f[x,y],g[x,y]},{x,xmin,xmax},{y,ymin,ymax},VectorStyle->{Black}],
ListPlot[{{x0,y0}},PlotStyle->{PointSize[.02],Red}] }]

Show[{Plot[1-x,{x, xmin, xmax},PlotStyle->{Red,Dashed,Thick}],ParametricPlot[{0,t},{t, ymin, ymax},PlotStyle->{Red,Dashed,Thick}], Plot[-1,{x, xmin, xmax},PlotStyle->{Purple,DotDashed,Thick}],Plot[0,{x, xmin, xmax},PlotStyle->{Purple,DotDashed,Thick}],

VectorPlot[{f[x,y],g[x,y]},{x,xmin,xmax},{y,ymin,ymax},VectorStyle->{Black}],
ListPlot[{{x0,y0}},PlotStyle->{PointSize[.02],Red}] }, PlotRange->{{xmin,xmax},{ymin,ymax}}]

solution to the associated linearized system near the critical point

check for linear independence of the first solution, written as a row vector,...

Clear[x1]
x1=Expand[ (E^(eval1 * t))evec1]

... and the second solution, written as a row vector,...

Clear[x2]
x2=Expand[(E^(eval2*t))evec2]

here is the matrix whose determinant is the Wronskian of the two solutions

MatrixForm[Transpose[ {x1,x2}]]

Det[ Transpose[ {x1,x2}]]

Since this non-zero constant multiple of an exponential function is never zero, our two solutions are linearly independent.  

Here is our linearized solution near this critical point:

Clear[LocalLinSoln]
LocalLinSoln = Expand[c1 x1 + c2 x2+{x0,y0}]

Clear[TrajectoryTable]
TrajectoryTable = Table[LocalLinSoln/.{c1->a,c2->b},{a,-.2,.2,.2},{b,-.2,.2,.2}]

Show[{ParametricPlot[TrajectoryTable,{t,0,1.5}],
ListPlot[TrajectoryTable/.{t->0},PlotStyle->{PointSize[.02],Black}] ,
ListPlot[{{x0,y0}},PlotStyle->{PointSize[.02],Red}]  }]

Show[{
VectorPlot[{f[x,y],g[x,y]},{x,xmin-.25,xmax+.25},{y,ymin-.25,ymax+.25},PlotRange->{{xmin-.25,xmax+.25},{ymin-.25,ymax+.25}}],
ListPlot[{{x0,y0}},PlotStyle->{PointSize[.02],Red}],
ParametricPlot[TrajectoryTable,{t,0,1.5}],
ListPlot[TrajectoryTable/.{t->0},PlotStyle->{PointSize[.02],Black}]  }]