phaseportpoisson2.html

Plot phase portrait for moment closure equations for Verhulst logistic model, assuming Poisson distribution.

This Maple work sheet uses Maple8, and the new package LinearAlgebra. Date 030402.

> restart;
with(DEtools,DEplot);
with(plots,display);
with(LinearAlgebra,Eigenvectors);
with(VectorCalculus,Jacobian);
XP:=Vector(2);

In the original Verhulst model, take alpha1=alpha2=1, R0=4, N=100. Then a1=3/5, f1=5/3, f2=3/5. Take r/K=1 (by rescaling time). Put kappa1=X[1], kappa2=X[2].

The equations are given by the procedure G:

>    G:=proc(N,t,X,XP)
  XP[1]:=(60-X[1])*X[1]-X[2];
  XP[2]:=(100-0.6*X[1])*X[1]+119.4*X[2]-4*X[1]*X[2]-2*X[1];
end;

Run G to determine XP:

> G(1,t,X,XP):

We determine the isoclines where the derivatives are equal to zero:

> isocline1:=solve(XP[1],X[2]);
isocline2:=solve(XP[2],X[2]);

Plot the isoclines:

> plot({isocline1,isocline2},X[1]=-10..100,y=-100..1000,discont=true);

> expand(subs(X[2]=isocline1,XP[2]));

> plot(%,X[1]=-10..80);

There are three critical points! Solve for them:

> critp:=solve({XP[1],XP[2]},{X[1],X[2]});

I use P1, P2, P3 to denote the three critical points. They are numbered from the left.
CHECK the numbering here!

> P1prel:=critp[1];
P2prel:=critp[2];
P3prel:=critp[3];

CHECK the ordering here!

> P1:=[rhs(P1prel[1]),rhs(P1prel[2])];
P2:=[rhs(P2prel[2]),rhs(P2prel[1])];
P3:=[rhs(P3prel[2]),rhs(P3prel[1])];

Now determine the Jacobian of the right-hand sides of the given system of equations:

> J:=Jacobian(XP,[X[1],X[2]]);

Determine the jacobian in each of the critical points:

> J1:=subs(P1prel,J);

> J2:=subs(P2prel,J);

> J3:=subs(P3prel,J);

The eigenvalues and eigenvectors of these three matrices:

> eg1:=Eigenvectors(J1,output=list);

> eg2:=Eigenvectors(J2,output=list);

> eg3:=Eigenvectors(J3,output=list);

Notice that P1 is an unstable node, P2 is a saddle point, P3 is a stable node. Furthermore,

CHECK the numbering here!

eg1[1] : slow manifold, eg1[2]: fast manifold of unstable node P0. Both manifolds have negative slopes.

eg2[1]: unstable manifold, eg2[2]: stable manifold of saddle point P1.

eg3[1]: slow manifold, eg3[2]: fast manifold of stable node P2.

A Maple procedure for determining initial values on each manifold near each critical point...

Think of it as transforming eg-to-init (eg2init):

> eg2init:=proc(P,eg,d)
local u,lu,v,k;
u:=eg[3]; # This is an eigenvector.
lu:=sqrt(u[1]^2+u[2]^2); # This is the length of the eigenvector u.
v:=d*u/lu; # This is a vector of length d in the direction of u.
[seq([x(0)=P[1]+(-1)^k*v[1],y(0)=P[2]+(-1)^k*v[2]],k=1..2)];
#These are two initial values near P in the direction of the eigenvector.
end proc:

CHECK the numbering here!

> init2unstable:=eg2init(P2,eg2[1],10);
init2stable:=eg2init(P2,eg2[2],10);
init3fast:=eg2init(P3,eg3[2],20);

Plot two orbits on the unstable manifold of P2:

> plot2unstable:=DEplot(G,[x(t),y(t)],t=0..0.2,number=2,init2unstable,x=-10..100,y=-100..1000,arrows=none,stepsize=0.01,linecolor=red):
plot2unstable;

Plot two orbits on the stable manifold of the saddle P2:

> plot2stable:=DEplot(G,[x(t),y(t)],t=-0.2..0,number=2,init2stable,x=-10..100,y=-100..1000,arrows=none,stepsize=0.01,linecolor=red):
plot2stable;

Plot two orbits on the fast manifold of P3:

> plot3fast:=DEplot(G,[x(t),y(t)],t=-0.05..0,number=2,init3fast,x=-10..100,y=-100..1000,arrows=none,linecolor=red,stepsize=0.01):
plot3fast;

Plot one orbit on the slow manifold of P3:

> plot3slow:=DEplot(G,[x(t),y(t)],t=0..0.05,number=2,{[x(0)=500,y(0)=20]},x=-10..100,y=-100..1000,stepsize=0.001,arrows=none,linecolor=red):
plot3slow;

NOTE: The arrows on the final phase portrait become less pronounced if the commend "color=red" below is replaced by "color=yellow". Plot the direction field:

> plotarrows:=DEplot(G,[x(t),y(t)],t=0..0.1,number=2,{[x(0)=60,y(0)=31]},x=-10..100,y=-100..1000,stepsize=0.01,color=red,linecolor=red):
plotarrows;

Plot three additional orbits:

> init1:=[x(0)=10,y(0)=0]:
plot1:=DEplot(G,[x(t),y(t)],t=0..0.1,number=2,[init1],x=-10..100,y=-100..1000,stepsize=0.01,arrows=none,linecolor=red):
plot1;

> init2:=[x(0)=50,y(0)=1000]:
plot2:=DEplot(G,[x(t),y(t)],t=0..0.1,number=2,[init2],x=-10..100,y=-100..1000,stepsize=0.01,arrows=none,linecolor=red):
plot2;

> init3:=[x(0)=100,y(0)=1000]:
plot3:=DEplot(G,[x(t),y(t)],t=0..0.05,number=2,[init3],x=-10..100,y=-100..1000,stepsize=0.002,arrows=none,linecolor=red):
plot3;

Final plot of all orbits:

> display({plot2stable,plot2unstable,plot3fast,plot3slow,plotarrows,plot1,plot2,plot3},labels=[Mean,Variance],title=Poisson_Assumption);

>