binomial2.mws  030128

Determine the critical points of the first two moment closure equations, and their stabilities, under assumption of binomial distribution!

>

restart;

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

The rhss of the differential equations for the first two cumulants are eq1 and eq2prel. I reserve the name eq2 for the rhs of 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 binomial distribution has kappa3=2*kappa2^2/kappa1-kappa2.  Inserting this expression for kappa3 into eq2prel gives:

>	eq2:=subs(kappa3=2*kappa2^2/kappa1-kappa2,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 point by (x[2],y[2]).

>	x:=solve(eq2a,kappa1);

The first of these x-values is not a critical point. The reason is that the second rhs (eq2) is not defined for kappa1=0.

The kappa2-value at the one critical point is:

>	y[2]:=normal(subs(kappa1=x[2],kap2));

Asymptotic approximations of the coordinates of the critical point:

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

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

Study the stability of the critical point:

>	F:=[eq1,eq2];

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

Determine the Jacobian at the critical point:

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

The eigenvalues are:

>	eg2:=Eigenvalues(J2);

Asymptotic approximations of the eigenvalues:

>	assume(a1>0);

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

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

>

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