Notebook Author/Course Instructor: Barry McQuarrie
Date: Nov 2015
Lecture topic: Competing Species

In General (some help with eigenvalues and linearizing)

Critical Points

f1[x_,y_]=x((Subscript[ϵ, 1])-(Subscript[σ, 1])x-(Subscript[α, 1])y)
f2[x_,y_]=y((Subscript[ϵ, 2])-(Subscript[σ, 2])y-(Subscript[α, 2])x)

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

The interesting point is the one with nonzero x and y. Let's investigate that one.

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 given by:

Eigenvalues[jac]

So the point (0,0) is an unstable node in all cases.

Linearizing about Equilibrium Point (((Subscript[ϵ, 1])/(Subscript[σ, 1])),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->((Subscript[ϵ, 1])/(Subscript[σ, 1])),y->0};
MatrixForm[jac]

The Jacobian matrix has eigenvalues given by:

Eigenvalues[jac]

Linearizing about Equilibrium Point (0,((Subscript[ϵ, 2])/(Subscript[σ, 2])))

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->((Subscript[ϵ, 2])/(Subscript[σ, 2]))};
MatrixForm[jac]

The Jacobian matrix has eigenvalues given by:

Eigenvalues[jac]

Linearizing about Equilibrium Point (-(((Subscript[α, 1]) (Subscript[ϵ, 2])-(Subscript[ϵ, 1]) (Subscript[σ, 2]))/(-(Subscript[α, 1]) (Subscript[α, 2])+(Subscript[σ, 1]) (Subscript[σ, 2]))),-((-(Subscript[α, 2]) (Subscript[ϵ, 1])+(Subscript[ϵ, 2]) (Subscript[σ, 1]))/((Subscript[α, 1]) (Subscript[α, 2])-(Subscript[σ, 1]) (Subscript[σ, 2]))))

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->(Subscript[X, 0]),y->(Subscript[Y, 0])}
MatrixForm[jac]

We can simplify this equation using the fact that this is an equilibrium point.

jac={{-(Subscript[X, 0]) (Subscript[σ, 1]),-(Subscript[X, 0]) (Subscript[α, 1])},{-(Subscript[Y, 0]) (Subscript[α, 2]),-(Subscript[Y, 0]) (Subscript[σ, 2])}}
MatrixForm[jac]

The Jacobian matrix has eigenvalues given by:

Eigenvalues[jac]

To classify the behaviour, we need to know if the quantity under the radical is positive, negative, or zero.

Example

x'(t)=(Subscript[f, 1])(x,y)=x(1.5-x-0.5y)
y'(t)=(Subscript[f, 2])(x,y)=y(2-y-0.75x)

Critical Points

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

f1[x_,y_]=x(1.5-x-0.5y)
f2[x_,y_]=y(2-y-0.75x)
NSolve[{f1[x,y]==0,f2[x,y]==0},{x,y}]

As expected, there are four critical points for the system. Since (Subscript[σ, 1]) (Subscript[σ, 2])-(Subscript[α, 1])(Subscript[α, 2])=(1)(1)-(0.5)(0.75)>0, we expect a stable population at (0.8,1.4).

Direction field to see if we got the right behaviour

field=VectorPlot[{f1[x,y],f2[x,y]},{x,0,2},{y,0,2},VectorScale->{Tiny,Small,None},VectorPoints->20];

scale=0.05;
points=Graphics[{Red,Disk[{0.8,1.4},scale],Disk[{1.5,0},scale],Disk[{0,2},scale],Disk[{0,0},scale]}];

Show[field,points,AspectRatio->1]

initialpoint={};
delta=0.1;
Do[initialpoint=Append[initialpoint,{0.1,0.1+delta*i}],{i,0,20}];
Do[initialpoint=Append[initialpoint,{2.0,0.1+delta*i}],{i,0,20}];
Do[initialpoint=Append[initialpoint,{0.1+delta*i,0.1}],{i,0,20}];
Do[initialpoint=Append[initialpoint,{0.1+delta*i,2.0}],{i,0,20}];

Clear[f1,f2,t,x1,x2,h,taylor2plot]
SeedRandom[1];
listplots={};
f1[t_,x1_,x2_]=x1(1.5-x1-0.5x2);
f2[t_,x1_,x2_]=x2(2-x2-0.75x1);

h=0.01;
numivp=Length[initialpoint];
Do[Clear[t,x1,x2];t[0]=0;x1[0]=initialpoint[[k,1]];x2[0]=initialpoint[[k,2]];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,400}];taylor2plot[k_]=ListPlot[taylor2,PlotStyle->{Thick,Blue},AxesLabel->{"(Subscript[x, 1])","(Subscript[x, 2])"},PlotRange->{{0,2},{0,2}},Joined->True];listplots=Append[listplots,taylor2plot[k]],{k,1,numivp}]

plot1=Show[listplots,points,field,ImageSize->{600,Automatic},AspectRatio->1]

Example

x'(t)=(Subscript[f, 1])(x,y)=x(1.5-0.5x-y)
y'(t)=(Subscript[f, 2])(x,y)=y(2-y-x)

Critical Points

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

f1[x_,y_]=x(1.5-0.5x-y)
f2[x_,y_]=y(2-y-x)
NSolve[{f1[x,y]==0,f2[x,y]==0},{x,y}]

As expected, there are four critical points for the system. Since (Subscript[σ, 1]) (Subscript[σ, 2])-(Subscript[α, 1])(Subscript[α, 2])=(1)(0.5)-(1)(1)<0, we expect an unstable equilibrium population at (1,1).

Direction field to see if we got the right behaviour

field=VectorPlot[{f1[x,y],f2[x,y]},{x,0,3},{y,0,3},VectorScale->{Tiny,Small,None},VectorPoints->20];

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

Show[field,points,AspectRatio->1]

initialpoint={};
delta=0.1;
Do[initialpoint=Append[initialpoint,{0.1,0.1+delta*i}],{i,0,20}];
Do[initialpoint=Append[initialpoint,{2.0,0.1+delta*i}],{i,0,20}];
Do[initialpoint=Append[initialpoint,{0.1+delta*i,0.1}],{i,0,20}];
Do[initialpoint=Append[initialpoint,{0.1+delta*i,2.0}],{i,0,20}];

Clear[f1,f2,t,x1,x2,h,taylor2plot]
SeedRandom[1];
listplots={};
f1[t_,x1_,x2_]=x1(1.5-0.5x1-x2);
f2[t_,x1_,x2_]=x2(2-x2-x1);

h=0.01;
numivp=Length[initialpoint];
Do[Clear[t,x1,x2];t[0]=0;x1[0]=initialpoint[[k,1]];x2[0]=initialpoint[[k,2]];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,400}];taylor2plot[k_]=ListPlot[taylor2,PlotStyle->{Thick,Blue},AxesLabel->{"(Subscript[x, 1])","(Subscript[x, 2])"},PlotRange->{{0,3},{0,3}},Joined->True];listplots=Append[listplots,taylor2plot[k]],{k,1,numivp}]

plot1=Show[listplots,points,field,ImageSize->{600,Automatic},AspectRatio->1]