Determine the critical points of the first two moment closure equations for the Verhulst model, and their stabilities, under assumption of normality!

>

restart;

>	with(VectorCalculus,Jacobian); with(LinearAlgebra,Eigenvalues);

The differential equations for the first two cumulants are eq1 and eq2prel. I reserve the name eq2 for the second moment closure equation.

>	eq1:=(a1*N-kappa1)*kappa1-kappa2;

>	eq2prel:=(f1*a1*N-f2*kappa1)*kappa1+(2*a1*N-f2)*kappa2-4*kappa1*kappa2-2*kappa3;

The normality assumption means that kappa3=0:

>	eq2:=subs(kappa3=0,eq2prel);

Solve the first equation for kappa2:

>	kap2:=solve(eq1,kappa2);

Substitute this value of kappa2 into the second equation.

>	eq2a:=subs(kappa2=kap2,eq2);

Denote the critical points by (x[1],y[1]), (x[2],y[2]), (x[3],y[3]).

>	x:=solve(eq2a,kappa1);

These are the kappa1-coordinates of the three critical points. The corresponding kappa2-coordinates are:

>	for k from 1 to 3 do y[k]:=expand(subs(kappa1=x[k],kap2)) od;

Asymptotic approximations of the critical points:

>	assume(a1>0);

>	x2asymp:=map(simplify,convert(asympt(x[2],N,2),polynom));

>	x3asymp:=map(simplify,convert(asympt(x[3],N,2),polynom));

>	y2asymp:=map(simplify,convert(asympt(y[2],N,2),polynom));

>	y3asymp:=map(simplify,convert(asympt(y[3],N,1),polynom));

>

a1:='a1':

Study the stability of the critical points:

>	F:=[eq1,eq2];

>	J:=Jacobian(F,[kappa1,kappa2]);

Determine the Jacobian in each of the three critical points:

>	J1:=subs(kappa1=x[1],kappa2=y[1],eval(J));

>	J2:=subs(kappa1=x[2],kappa2=y[2],eval(J));

>	J3:=subs(kappa1=x[3],kappa2=y[3],eval(J));

Determine the eigenvalues of each of these 3 matrices:

>	eg1:=Eigenvalues(J1);

>	eg2:=Eigenvalues(J2);

>	eg3:=Eigenvalues(J3);

Asymptotic approximations of the eigenvalues:

>	assume(a1>0);

>	simplify(convert(asympt(eg1[1],N,1),polynom));

>	simplify(convert(asympt(eg1[2],N,1),polynom));

So both of the eigenvalues of the critical point (0,0) are positive. Hence (0,0) is an unstable node.

>	simplify(convert(asympt(eg2[1],N,1),polynom));

>	simplify(convert(asympt(eg2[2],N,1),polynom));

Both eigenvalues are negative. So this critical point is a stable node.

>	simplify(convert(asympt(eg3[1],N,1),polynom));

>	simplify(convert(asympt(eg3[2],N,1),polynom));

>

The third critical point has one positive and one negative eigenvalue. So it is a saddle point.