lognormal3.html

lognormal3.mws

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

> restart;

> with(LinearAlgebra,Eigenvalues);

> with(VectorCalculus,Jacobian);

The differential equations for the first three cumulants are denoted eq1, eq2, eq3prel:

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

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

> eq3prel:=(a1*N-kappa1)*kappa1-6*f2*kappa1*kappa2-6*kappa1*kappa3+(3*f1*a1*N-1)*kappa2-6*kappa2^2+3*(a1*N-f2)*kappa3-3*kappa4;

The log-normality assumption means that kappa4 is expressed in terms of kappa1 and kappa2 as follows (see kappa34.mws):

> eq3:=subs(kappa4=16*kappa2^3/kappa1^2+15*kappa2^4/kappa1^4+6*kappa2^5/kappa1^6+kappa2^6/kappa1^8,eq3prel);

Solve the first equation for kappa2. This gives kappa2 as a function of kappa1.

> kap2:=solve(eq1,kappa2);

Substitute this value of kappa2 into the second equation.

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

Solve this equation for kappa3. This gives kappa3 as a function of kappa1.

> kap3:=solve(eq2a,kappa3);

Insert the values of kappa2 and kappa3 into eq3; then multiply by -1/3. This gives an expression in kappa1.

> eq3a:=-map(simplify,collect(subs(kappa2=kap2,kappa3=kap3,eq3),kappa1))/3;

The problem now is to solve this forth-degree equation in kappa1. (Some minor (cosmetic) simplification can be achieved by introducing K and sigma^2.) The equation to be solved can be written

> eq3b:=collect(kappa1^2*eq3a,kappa1);

Look for asymptotic approximation! Put kappa1 = AN + B + C/N.

> eq3c:=subs(kappa1=A*N+B+C/N,eq3b);

Form 3 equations (for A, B, C) by putting the coefficients of N^6, N^5, N^4 equal to zero.

> eqA:=coeff(eq3c,N,6);

> eqB:=coeff(eq3c,N,5);

> eqC:=coeff(eq3c,N,4);

To study eqA: Scaling, put A/a1 = xx.

> eqAa:=expand(subs(A=a1*xx,eqA/a1^6));

> plot(eqAa,xx=0..1.2);

The equation eqAa=0 has two real roots. It is easy to see that xx=1 is a root.

> subs(xx=1,eqAa);

> XX:=fsolve(eqAa,xx);

The corresponding A-values are:

> AA:=a1*XX[1],a1;

Now solve equations eqB and eqC for B and C, corresponding to these two A-values.

> eqB1:=subs(A=AA[1],eqB);

> BB[1]:=solve(eqB1,B);

> eqC1:=simplify(subs(A=AA[1],B=BB[1],eqC));

> Csol:=solve(eqC1,C);

> CC[1]:=expand(op(1,Csol)*op(2,Csol))*op(3,Csol);

This can be simplified a bit by noting that f1-f2 equals 2*sigma^2/K. Note:

> 0.6371823775+0.2415082880;

So the numerator of CC[1] can be written (f1-f2)*(0.40*f1+0.24*(f1-f2)).

An asymptotic approximation for the first solution for kappa1 has thus been found. It has the form

AA[1]*N + BB[1] + CC[1]/N. Use it to derive asymptotic approximations of kappa2 and kappa3, using kap2 and kap3.

Write kappa2 sim DD[1]*N^2+EE[1]*N and kappa3 sim FF[1]*N^3.

Recalling kap2:

> kap2;

Substitute the first kappa1-value into kap2:

> kap21:=expand(convert(asympt(subs(kappa1=AA[1]*N+BB[1]+CC[1]/N,kap2),N,1),polynom));

> DD[1]:=coeff(kap21,N,2);

> EE[1]:=simplify(coeff(kap21,N,1));

> kap31:=map(expand,convert(asympt(subs(kappa1=AA[1]*N+BB[1]+CC[1]/N,kap3),N,1),polynom));

> FF[1]:=coeff(kap31,N,3);

Now proceed to find approximations of the three components of the remaining critical point, with subscript 2. This is the critical point that I expect to be stable.

> eqB2:=expand(subs(A=AA[2],eqB));

> BB[2]:=solve(eqB2,B);

> eqC2:=simplify(subs(A=AA[2],B=BB[2],eqC));

> CC[2]:=solve(eqC2,C);

> kap22:=expand(convert(asympt(subs(kappa1=AA[2]*N+BB[2]+CC[2]/N,kap2),N,2),polynom));

> DD[2]:=factor(coeff(kap22,N,1));

> EE[2]:=factor(simplify(coeff(kap22,N,0)));

So kappa2 is sim DD[2]*N + EE[2].

> kap32:=map(simplify,convert(asympt(subs(kappa1=AA[2]*N+BB[2]+CC[2]/N,K=a1*N,kap3),N,2),polynom));

> FF[2]:=coeff(kap32,N,1);

Here, kappa3 is sim FF[2]*N.

So I have derived asymptotic approximatons of the first three cumulants for the two critical points...

I would also like to study the local stability of the 2 critical points. I do this by evaluating the eigenvalues of the jacobian matrices evaluated at each of the four critical points...

> F:=[eq1,eq2,eq3];

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

Determine the Jacobian in each of the 2 critical points:

> J1:=subs(kappa1=AA[1]*N+BB[1]+CC[1]/N,kappa2=DD[1]*N^2+EE[1]*N,kappa3=FF[1]*N^3,J):

Modify the matrix J1, replacing each term by the first 2 terms in its asymptotic approximation. The modified matrix is denoted J1a.

> J1a:=J1:

> J1a[1,1]:=convert(asympt(J1[1,1],N,1),polynom);

> J1a[1,2];

> J1a[1,3];

> J1a[2,1]:=convert(asympt(J1[2,1],N,0),polynom);

> J1a[2,2]:=convert(asympt(J1[2,2],N,1),polynom);

> J1a[2,3];

> J1a[3,1]:=simplify(convert(asympt(J1[3,1],N,11),polynom));

> J1a[3,2]:=simplify(convert(asympt(J1[3,2],N,10),polynom));

> J1a[3,3]:=convert(asympt(J1[3,3],N,1),polynom);

Now determine the eigenvalues of the modified matrix J1a.

> eg1:=Eigenvalues(J1a):

> a1:='a1':
eg11a:=convert(asympt(eg1[1],N,3),polynom):

> assume(a1>0);
eg11b:=simplify(evalc(eg11a));

> a1:='a1':
eg12a:=convert(asympt(eg1[2],N,3),polynom):

> assume(a1>0);
eg12b:=simplify(evalc(eg12a));

> a1:='a1':
eg13a:=convert(asympt(eg1[3],N,3),polynom):

> assume(a1>0);
eg13b:=simplify(evalc(eg13a));

> a1:='a1':

The eigenvalues of J1 are asymptotically 1.37K, -4.80K, -K. Hence the corresponding critical point is unstable.

It now remains to determine J2 and its eigenvalues.

> J2:=subs(kappa1=AA[2]*N+BB[2]+CC[2]/N,kappa2=DD[2]*N+EE[2],kappa3=FF[2]*N,J):

Maple can determine the eigenvalues of this matrix, but the expressions for the eigenvalues are very long.

Proceed to modify the elements of J2, replacing each of them by the first two terms in its asymptotic approximation. The modified matrix is called J2a.

> J2a:=J2:

> J2a[1,1]:=convert(asympt(J2[1,1],N,1),polynom);