Math 2401 Differential Equations

Notebook Author/Course Instructor: Barry McQuarrie
Date: Nov 2015
Lecture topic: Locally Linear Systems

Example 1

x'(t)=(Subscript[f, 1])(x,y)=4-2y
y'(t)=(Subscript[f, 2])(x,y)=12-3(x^2)

Autonomous Solution (loses t dependence)

The system is autonomous, and when written as a ratio (dy/dt)/(dx/dt) it is straightforward to separate and integrate to get:

family=12 x-(x^3)==4y-(y^2)+c;

familylist=family/.c->Range[-20,20,4];

lines=ContourPlot[Evaluate[familylist],{x,-5,5},{y,-5,10},PlotRange->{{-5,5},{-5,10}},FrameLabel->{"x","y"},ContourStyle->Purple,BaseStyle->20]

Critical Points

We can solve for critical points, which is where the derivatives are zero:

f1[x_,y_]=4-2y
f2[x_,y_]=12-3(x^2)
Solve[{f1[x,y]==0,f2[x,y]==0},{x,y}]

Note that these are the points we saw interesting things happening before. We can expand everything in Taylor series about each of these points in turn, and effectively linearize the nonlinear system when we are close to a given stationary point. We could then look at eigenvalues of the matrix A for the linearized problem to classify the points of interest.

Linearizing about Equilibrium Point (2,2)

jacobian[x_,y_]={{D[f1[x,y],x],D[f1[x,y],y]},{D[f2[x,y],x],D[f2[x,y],y]}};
MatrixForm[jacobian[x,y]]
MatrixForm[jacobian[2,2]]

The Jacobian matrix has eigenvalues that are real and differ in sign.

Eigenvalues[jacobian[2,2]]

We expect the linear system to have a saddle point at (2,2), and the nonlinear system to behave similarly near (2,2).

CP1=Graphics[{Red,Disk[{2,2},0.1]}];

stream=StreamPlot[{f1[x,y],f2[x,y]},{x,0,4},{y,0,4},VectorScale->{Tiny,Small,None},FrameLabel->{"x","y"},BaseStyle->{26},StreamPoints->40,StreamStyle->{{Orange,Thick}},Axes->True,AxesOrigin->{0,0}];

streamlinear=StreamPlot[jacobian[2,2].{x-2,y-2},{x,0,4},{y,0,4},VectorScale->{Tiny,Small,None},FrameLabel->{"x","y"},BaseStyle->{26},StreamPoints->40,StreamStyle->{{Green,Thick}},Axes->True,AxesOrigin->{0,0}];

Show[stream,streamlinear,CP1]

Linearizing about Equilibrium Point (-2,2)

jacobian[x_,y_]={{D[f1[x,y],x],D[f1[x,y],y]},{D[f2[x,y],x],D[f2[x,y],y]}};
MatrixForm[jacobian[x,y]]
MatrixForm[jacobian[-2,2]]

The Jacobian matrix has eigenvalues that are complex conjugates.

Eigenvalues[jacobian[-2,2]]

We expect the solution to have a center or spiral point at (-2,2). We can't determine if we have a center or a spiral point since the locally linear system is very sensitive to small changes about critical points where the eigenvalues are purely complex valued.

CP1=Graphics[{Red,Disk[{-2,2},0.1]}];

stream=StreamPlot[{f1[x,y],f2[x,y]},{x,-4,0},{y,0,4},VectorScale->{Tiny,Small,None},FrameLabel->{"x","y"},BaseStyle->{26},StreamPoints->40,StreamStyle->{{Orange,Thick}},Axes->True,AxesOrigin->{0,0}];

streamlinear=StreamPlot[jacobian[-2,2].{x+2,y-2},{x,-4,0},{y,0,4},VectorScale->{Tiny,Small,None},FrameLabel->{"x","y"},BaseStyle->{26},StreamPoints->40,StreamStyle->{{Green,Thick}},Axes->True,AxesOrigin->{0,0}];

Show[stream,streamlinear,CP1]

Example 2 Van der Pol Oscillator

x'=y
 y'=-x+(1-(x^2))y

Critical Points

We can solve for critical points, which is where the derivatives are zero:

f1[x_,y_]=y
f2[x_,y_]=-x+(1-(x^2))y
Solve[{f1[x,y]==0,f2[x,y]==0},{x,y}]

Linearizing about Equilibrium Point (0,0)

jacobian[x_,y_]={{D[f1[x,y],x],D[f1[x,y],y]},{D[f2[x,y],x],D[f2[x,y],y]}};
MatrixForm[jacobian[x,y]]
MatrixForm[jacobian[0,0]]

The Jacobian matrix has eigenvalues that are real and differ in sign.

Eigenvalues[jacobian[0,0]]//N

We expect the linear system to have a spiral point at (0,0).

CP1=Graphics[{Red,Disk[{0,0},0.05]}];

stream=StreamPlot[{f1[x,y],f2[x,y]},{x,-1,1},{y,-1,1},VectorScale->{Tiny,Small,None},FrameLabel->{"x","y"},BaseStyle->{26},StreamPoints->40,StreamStyle->{{Orange,Thick}},Axes->True,AxesOrigin->{0,0}];

streamlinear=StreamPlot[jacobian[0,0].{x-0,y-0},{x,-1,1},{y,-1,1},VectorScale->{Tiny,Small,None},FrameLabel->{"x","y"},BaseStyle->{26},StreamPoints->40,StreamStyle->{{Green,Thick}},Axes->True,AxesOrigin->{0,0}];

Show[stream,streamlinear,CP1]

Example 3

x'(t)=(Subscript[f, 1])(x,y)=4-2 (y^2)-((y x)/4)
y'(t)=(Subscript[f, 2])(x,y)=12-3 (x^2)-y

Critical Points

We can solve for critical points, which is where the derivatives are zero:

f1[x_,y_]=4-2 (y^2)-y*x/4
f2[x_,y_]=12-3 (x^2)-y
NSolve[{f1[x,y]==0,f2[x,y]==0},{x,y}]

There are four critical points for the system.

Linearizing about Equilibrium Point ( 2.1255, -1.55328)

jacobian[x_,y_]={{D[f1[x,y],x],D[f1[x,y],y]},{D[f2[x,y],x],D[f2[x,y],y]}};
MatrixForm[jacobian[x,y]]
MatrixForm[jacobian[2.1255,-1.55328]]

The Jacobian matrix has eigenvalues that are real and differ in sign.

Eigenvalues[jacobian[2.1255,-1.55328]]

We expect the linear system to have a stable spiral point at ( 2.1255, -1.55328).

CP1=Graphics[{Red,Disk[{2.1255,-1.55328},0.02]}];

stream=StreamPlot[{f1[x,y],f2[x,y]},{x,1.5,2.5},{y,-2,-1},VectorScale->{Tiny,Small,None},FrameLabel->{"x","y"},BaseStyle->{26},StreamPoints->15,StreamStyle->{{Orange,Thick}},Axes->True,AxesOrigin->{0,0}];

streamlinear=StreamPlot[jacobian[2.1255,-1.55328].{x-2.1255,y+1.55328},{x,1.5,2.5},{y,-2,-1},VectorScale->{Tiny,Small,None},FrameLabel->{"x","y"},BaseStyle->{26},StreamPoints->15,StreamStyle->{{Green,Thick}},Axes->True,AxesOrigin->{0,0}];

Show[stream,streamlinear,CP1]

Linearizing about Equilibrium Point ( -1.8676410938338404`,1.5357502338694027`)

jacobian[x_,y_]={{D[f1[x,y],x],D[f1[x,y],y]},{D[f2[x,y],x],D[f2[x,y],y]}};
MatrixForm[jacobian[x,y]]
MatrixForm[jacobian[-1.86764,1.53575]]

The Jacobian matrix has eigenvalues that are real and differ in sign.

Eigenvalues[jacobian[-1.86764,1.53575]]

We expect the linear system to have a stable spiral point at ( -1.8676410938338404`,1.5357502338694027`).

CP1=Graphics[{Red,Disk[{-1.86764,1.53575},0.02]}];

stream=StreamPlot[{f1[x,y],f2[x,y]},{x,-2.5,-1.5},{y,0,2},VectorScale->{Tiny,Small,None},FrameLabel->{"x","y"},BaseStyle->{26},StreamPoints->15,StreamStyle->{{Orange,Thick}},Axes->True,AxesOrigin->{0,0}];

streamlinear=StreamPlot[jacobian[-1.86764,1.53575].{x+1.86764,y-1.53575},{x,-2.5,-1.5},{y,0,2},VectorScale->{Tiny,Small,None},FrameLabel->{"x","y"},BaseStyle->{26},StreamPoints->15,StreamStyle->{{Green,Thick}},Axes->True,AxesOrigin->{0,0}];

Show[stream,streamlinear,CP1]

Linearizing about Equilibrium Point ( -2.1046596500279016`,-1.288776727366705`)

jacobian[x_,y_]={{D[f1[x,y],x],D[f1[x,y],y]},{D[f2[x,y],x],D[f2[x,y],y]}};
MatrixForm[jacobian[x,y]]
MatrixForm[jacobian[ -2.10466,-1.28878]]

The Jacobian matrix has eigenvalues that are real and differ in sign.

Eigenvalues[jacobian[-2.10466,-1.28878]]

We expect the linear system to have a saddle point at ( -2.1046596500279016`,-1.288776727366705`).

CP1=Graphics[{Red,Disk[{-2.10466,-1.28878},0.02]}];

stream=StreamPlot[{f1[x,y],f2[x,y]},{x,-2.5,-1.5},{y,-2,-1},VectorScale->{Tiny,Small,None},FrameLabel->{"x","y"},BaseStyle->{26},StreamPoints->15,StreamStyle->{{Orange,Thick}},Axes->True,AxesOrigin->{0,0}];

streamlinear=StreamPlot[jacobian[-2.10466,-1.28878].{x+2.10466,y+1.28878},{x,-2.5,-1.5},{y,-2,-1},VectorScale->{Tiny,Small,None},FrameLabel->{"x","y"},BaseStyle->{26},StreamPoints->15,StreamStyle->{{Green,Thick}},Axes->True,AxesOrigin->{0,0}];

Show[stream,streamlinear,CP1]

Linearizing about Equilibrium Point (1.88846,1.3011)

jacobian[x_,y_]={{D[f1[x,y],x],D[f1[x,y],y]},{D[f2[x,y],x],D[f2[x,y],y]}};
MatrixForm[jacobian[x,y]]
MatrixForm[jacobian[ 1.88846,1.3011]]

The Jacobian matrix has eigenvalues that are real and differ in sign.

Eigenvalues[jacobian[1.88846,1.3011]]

We expect the linear system to have a saddle point at (1.88846,1.3011).

CP1=Graphics[{Red,Disk[{1.88846,1.3011},0.02]}];

stream=StreamPlot[{f1[x,y],f2[x,y]},{x,1.5,2.5},{y,1,2},VectorScale->{Tiny,Small,None},FrameLabel->{"x","y"},BaseStyle->{26},StreamPoints->15,StreamStyle->{{Orange,Thick}},Axes->True,AxesOrigin->{0,0}];

streamlinear=StreamPlot[jacobian[1.88846,1.3011].{x-1.88846,y-1.3011},{x,1.5,2.5},{y,1,2},VectorScale->{Tiny,Small,None},FrameLabel->{"x","y"},BaseStyle->{26},StreamPoints->15,StreamStyle->{{Green,Thick}},Axes->True,AxesOrigin->{0,0}];

Show[stream,streamlinear,CP1]

Sketch using NDSolve to see if we got the right behaviour

Clear[f1,f2,t,x1,x2,h,taylor2plot]
SeedRandom[1];
listplots={};
f1[t_,x1_,x2_]=4-2 (x2^2)-x2*x1/4
f2[t_,x1_,x2_]=12-3 (x1^2)-x2

h=0.01;
numivp=100;
Do[
Clear[t,x1,x2];
t[0]=0;
x1[0]=RandomReal[{-4,4}];
x2[0]=RandomReal[{-4,4}];

t[i_]:=t[i]=t[i-1]+h;x1[i_]:=x1[i]=x1[i-1]+h f1[t[i-1],x1[i-1],x2[i-1]]+(1/2) (h^2) (Derivative[1,0,0][f1][t[i-1],x1[i-1],x2[i-1]]+Derivative[0,1,0][f1][t[i-1],x1[i-1],x2[i-1]] f1[t[i-1],x1[i-1],x2[i-1]]+Derivative[0,0,1][f1][t[i-1],x1[i-1],x2[i-1]] f2[t[i-1],x1[i-1],x2[i-1]]);x2[i_]:=x2[i]=x2[i-1]+h f2[t[i-1],x1[i-1],x2[i-1]]+(1/2) (h^2) (Derivative[1,0,0][f2][t[i-1],x1[i-1],x2[i-1]]+Derivative[0,1,0][f2][t[i-1],x1[i-1],x2[i-1]] f1[t[i-1],x1[i-1],x2[i-1]]+Derivative[0,0,1][f2][t[i-1],x1[i-1],x2[i-1]] f2[t[i-1],x1[i-1],x2[i-1]]);taylor2=Table[{x1[k],x2[k]},{k,0,200}];taylor2plot[k_]=ListPlot[taylor2,PlotStyle->{Purple},AxesLabel->{"(Subscript[x, 1])","(Subscript[x, 2])"},PlotRange->{{-4,4},{-4,4}},Joined->True];
listplots=Append[listplots,taylor2plot[k]],
{k,1,numivp}]

nullclines=ContourPlot[{4-2 (y^2)-y*x/4==0,12-3 (x^2)-y==0},{x,-5,5},{y,-5,5},PlotRange->{{-5,5},{-5,5}},FrameLabel->{"x","y"},ContourStyle->Blue,BaseStyle->20];

CP=Graphics[{Red,Disk[{2.1255,-1.55328},0.07],Disk[{-2.10466,-1.28878},0.07],Disk[{-1.86764,1.53575},0.07],Disk[{1.88846,1.3011},0.07]}];


plot1=Show[listplots,CP,nullclines,ImageSize->{600,Automatic},AspectRatio->1,BaseStyle->20]

Example 4

x'(t)=(Subscript[f, 1])(x,y)=(1+x)sin y
y'(t)=(Subscript[f, 2])(x,y)=1-x-cosy

Critical Points

We can solve for stationary points, which is where the derivatives are zero:

Clear[x,y]
f1[x_,y_]=(1+x)Sin[y]
f2[x_,y_]=1-x-Cos[y]
NSolve[{f1[x,y]==0,f2[x,y]==0},{x,y}]

There are four critical points for the system. Following Mathematica's advice, we try reduce:

Reduce[{f1[x,y]==0,f2[x,y]==0},{x,y}]

So there are two sets of solutions: (0,2n π) and (2,π+2n π).

Linearizing about Equilibrium Point (0,2n π)

jacobian[x_,y_]={{D[f1[x,y],x],D[f1[x,y],y]},{D[f2[x,y],x],D[f2[x,y],y]}};
MatrixForm[jacobian[x,y]]
MatrixForm[jacobian[0,2n π]]

It's useful to set global assumptions, which we can do using the following command:

$Assumptions={Element[n,Integers]}

Now, Mathematica can simplify since it knows that n is an integer:

MatrixForm[jacobian[0,2n π]]//Simplify

The Jacobian matrix has eigenvalues that are complex valued.

Eigenvalues[jacobian[0,2n π]]//Simplify

We would expect to have either a center or spiral point at (0,2n π).

Linearizing about Equilibrium Point (2,π+2n π)

jacobian[x_,y_]={{D[f1[x,y],x],D[f1[x,y],y]},{D[f2[x,y],x],D[f2[x,y],y]}};
MatrixForm[jacobian[x,y]]
MatrixForm[jacobian[2,π+2n π]]//Simplify

The Jacobian matrix has eigenvalues that are real and differ in sign.

Eigenvalues[jacobian[2,π+2n π]]//Simplify

We would expect to have an unstable saddle point at (2,π+2n π).

Stream Plot to see if we got the right behaviour

CP1=Table[Graphics[{Red,Disk[{0,2n π},0.2]}],{n,-5,5}];
CP2=Table[Graphics[{Blue,Disk[{2,π+2n π},0.2]}],{n,-5,5}];

stream=StreamPlot[{f1[x,y],f2[x,y]},{x,-10,10},{y,-10,10},VectorScale->{Tiny,Small,None},FrameLabel->{"x","y"},BaseStyle->{26},StreamPoints->40,StreamStyle->{{Orange,Thick}},Axes->True,AxesOrigin->{0,0}];

Show[stream,CP1,CP2]

Example That Is Not Locally Linear

x'(t)=(Subscript[f, 1])(x,y)=(y^2)
y'(t)=(Subscript[f, 2])(x,y)=(x^2)

Critical Points

We can solve for stationary points, which is where the derivatives are zero:

f1[x_,y_]=y^2
f2[x_,y_]=x^2
Solve[{f1[x,y]==0,f2[x,y]==0},{x,y}]

Linearizing about Equilibrium Point (0,0)

jacobian[x_,y_]={{D[f1[x,y],x],D[f1[x,y],y]},{D[f2[x,y],x],D[f2[x,y],y]}};
jac=jacobian[x,y]/.{x->0,y->0};
MatrixForm[jac]

The Jacobian matrix has eigenvalues that are real and differ in sign.

Eigenvalues[jac]

There is a problem here. If both eigenvalues are zero, then we can't classify the point. This is an example of a system that is NOT locally linear about (0,0)

Direction field

field=VectorPlot[{f1[x,y],f2[x,y]},{x,-1,1},{y,-1,1},VectorScale->{Tiny,Small,None},VectorPoints->20,StreamPoints->10,StreamStyle->{Orange,Thick,Opacity[0.6]},BaseStyle->20];

scale=0.05;
points=Graphics[{Red,Disk[{0,0},scale]}];

Show[field,points,AspectRatio->1]

Let's zoom in to see what exactly is happening near these points:

field=VectorPlot[{f1[x,y],f2[x,y]},{x,-0.0001,0.0001},{y,-0.0001,0.0001},VectorScale->{Tiny,Small,None},VectorPoints->20,StreamPoints->10,StreamStyle->{Orange,Thick,Opacity[0.6]},BaseStyle->20];

scale=0.000005;
points=Graphics[{Red,Disk[{0,0},scale]}];

Show[field,points,AspectRatio->1]

Notice that the behaviour is not turning into one of the linear cases as we zoom in!

field2[k_]:=VectorPlot[{f1[x,y],f2[x,y]},{x,-(10^(-k)),(10^(-k))},{y,-(10^(-k)),(10^(-k))},VectorScale->{Tiny,Small,None},VectorPoints->20,BaseStyle->20];

Manipulate[Show[field2[k],AspectRatio->1],{k,1,7,1}]

This of course isn't suprising, since there is no linear piece of the (Subscript[f, 1]) and (Subscript[f, 2])!