Lotka Volterra - Predator Prey

The system we are interested in is the following

(dx/dt)=ax-bxy
(dy/dt)=-cy+dxy

Null-Cline Analysis

As an example lets set a=1, b=2, c=3, d=4:

a=1;
b=2;
c=3;
d=4;

The null-clines can be found analytically:

Solve[a*x-b*x*y==0,{x,y}]
Solve[-c*y+d*x*y==0,{x,y}]

Lets plot the vector field and the null-clines:

plot1=VectorPlot[{a*x-b*x*y,-c*y+d*x*y},{x,0,2*c/d},{y,0,2*a/b}];
plot2=ParametricPlot[{{0,t},{t,a/b},{c/d,t},{t,0}},{t,0,2},PlotStyle->Thickness[.01]];
Show[plot1,plot2]

From the null-cline analysis we can see that the solution curves will rotate around the fixed point x=c/d, y=a/b.

Fixed Point Analysis

Solve[{a*x-b*x*y==0,-c*y+d*x*y==0},{x,y}]

Lets analyze the stability of these fixed points by determing the Jacobian matrix and its eigenvalues.

f[x_,y_]:=a*x-b*x*y;
g[x_,y_]:=-c*y+d*x*y;

MatrixForm[({{D[f[x,y],x], D[f[x,y],y]}, {D[g[x,y],x], D[g[x,y],y]}})]

J[x_,y_]:=(({{1-2 y, -2 x}, {4 y, -3+4 x}}))

Eigenvalues[J[0,0]]

Therefore the fixed point (0,0) is a saddle point.

Eigenvalues[J[3/4,1/2]]

According to linear stability analysis the fixed point is a center. This needs to be analyzed further.

Analysis of Center

Noticing that (dy/dx)=(-c*y+d*x*y)/(a*x-b*x*y) we can intergrate:

DSolve[y'[x]==(-c*y[x]+d*x*y[x])/(a*x-b*x*y[x]),y[x],x]

Solve[y==-(1/2) ProductLog[-((2 (E^(4 x-con)))/(x^3))],con]

So the solution curves are given by the contours of the function h(x,y)=(x^3)y*(e^(-2y))*(e^(-4x)).

ContourPlot[x^3*y*Exp[-2*y]*Exp[-4*x],{x,0,2*c/d},{y,0,2*a/b}]

The contours are actually closed curves so this proves that the orbits are periodic.

Numerical Analysis

s=NDSolve[{x'[t]==a*x[t]-b*x[t]*y[t],y'[t]==-c*y[t]+d*x[t]*y[t],x[0]==y[0]==0},{x,y},{t,20}];
plot3=ParametricPlot[Evaluate[{x[t],y[t]}/.s],{t,0,20}];
For[i=1,i<5,i++,
x0=1/i;
s=NDSolve[{x'[t]==a*x[t]-b*x[t]*y[t],y'[t]==-c*y[t]+d*x[t]*y[t],x[0]==y[0]==x0},{x,y},{t,20}];
plottemp=ParametricPlot[Evaluate[{x[t],y[t]}/.s],{t,0,20}];
plot3=Show[plottemp,plot3]];
Show[plot1,plot2,plot3]

Lets plot how x and y evolve as a function of time.

s=NDSolve[{x'[t]==a*x[t]-b*x[t]*y[t],y'[t]==-c*y[t]+d*x[t]*y[t],x[0]==y[0]==.1},{x,y},{t,20}];
Plot[Evaluate[{x[t],y[t]}/.s],{t,0,20},PlotStyle->Automatic]