general3.mws

general3.mws  030128

Determine the critical points of the first three moment closure equations for the Verhulst model, and their stabilities, under assumption that kappa4 = a4*N!

This Maple worksheet uses Maple 8, and the "new" package LinearAlgebra. 030226.

>    restart;

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

[Jacobian]

[Eigenvalues]

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

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

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

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;

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 assumption is that kappa4=a4*N:

>    eq3:=subs(kappa4=a4*N,eq3prel);

eq3 := (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*a4*N

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

>    kap2:=solve(eq1,kappa2);

kap2 := -(-a1*N+kappa1)*kappa1

Substitute this value of kappa2 into the second equation.

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

eq2a := (f1*a1*N-f2*kappa1)*kappa1-(2*a1*N-f2)*(-a1*N+kappa1)*kappa1+4*kappa1^2*(-a1*N+kappa1)-2*kappa3

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

>    kap3:=solve(eq2a,kappa3);

kap3 := 1/2*kappa1*f1*a1*N+kappa1*a1^2*N^2-3*a1*N*kappa1^2-1/2*kappa1*f2*a1*N+2*kappa1^3

Insert the values of kappa2 and kappa3 into eq3. Divide by -3. This gives an expression in kappa1.

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

eq3a := 6*kappa1^4-12*a1*N*kappa1^3+a1*N*(7*a1*N-2*f2+2*f1)*kappa1^2-1/2*a1*N*(3*f1*a1*N+2*a1^2*N^2-3*f2*a1*N-f2*f1+f2^2)*kappa1+a4*N
eq3a := 6*kappa1^4-12*a1*N*kappa1^3+a1*N*(7*a1*N-2*f2+2*f1)*kappa1^2-1/2*a1*N*(3*f1*a1*N+2*a1^2*N^2-3*f2*a1*N-f2*f1+f2^2)*kappa1+a4*N

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.)  Assume there is an asymptotic approximation of the form

kappa1 = A*N + B + C/N:

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

eq3b := 6*(A*N+B+C/N)^4-12*a1*N*(A*N+B+C/N)^3+a1*N*(7*a1*N-2*f2+2*f1)*(A*N+B+C/N)^2-1/2*a1*N*(3*f1*a1*N+2*a1^2*N^2-3*f2*a1*N-f2*f1+f2^2)*(A*N+B+C/N)+a4*N
eq3b := 6*(A*N+B+C/N)^4-12*a1*N*(A*N+B+C/N)^3+a1*N*(7*a1*N-2*f2+2*f1)*(A*N+B+C/N)^2-1/2*a1*N*(3*f1*a1*N+2*a1^2*N^2-3*f2*a1*N-f2*f1+f2^2)*(A*N+B+C/N)+a4*N

Determine the coefficients of N^4, N^3, N^2!  Putting them equal to zero gives equations for A, B, C.

>    eqA:=coeff(eq3b,N,4);

eqA := 6*A^4-12*a1*A^3+7*a1^2*A^2-a1^3*A

>    eqB:=coeff(eq3b,N,3);

eqB := 24*B*A^3-36*a1*B*A^2+a1*(-2*f2+2*f1)*A^2+14*a1^2*B*A-1/2*a1*(-3*f2*a1+3*f1*a1)*A-a1^3*B

>    eqC:=coeff(eq3b,N,2);

eqC := 12*(2*C*A+B^2)*A^2+24*B^2*A^2-12*a1*(C*A^2+2*B^2*A+A*(2*C*A+B^2))+2*a1*(-2*f2+2*f1)*B*A+7*a1^2*(2*C*A+B^2)-1/2*a1*(f2^2-f2*f1)*A-1/2*a1*(-3*f2*a1+3*f1*a1)*B-a1^3*C
eqC := 12*(2*C*A+B^2)*A^2+24*B^2*A^2-12*a1*(C*A^2+2*B^2*A+A*(2*C*A+B^2))+2*a1*(-2*f2+2*f1)*B*A+7*a1^2*(2*C*A+B^2)-1/2*a1*(f2^2-f2*f1)*A-1/2*a1*(-3*f2*a1+3*f1*a1)*B-a1^3*C

>    AA:=solve(eqA,A);

AA := 0, a1, (1/2+1/6*3^(1/2))*a1, (1/2-1/6*3^(1/2))*a1

>    evalf(%);

0., a1, .7886751347*a1, .2113248653*a1

Rearrange A-values in increasing order!

>    AA:=AA[1],AA[4],AA[3],AA[2];

AA := 0, (1/2-1/6*3^(1/2))*a1, (1/2+1/6*3^(1/2))*a1, a1

There are four critical points. Determine the corresponding values of A and B!

>    BB:=seq(solve(subs(A=AA[k],eqB),B),k=1..4);

BB := 0, 1/12*(-f2-f2*3^(1/2)+f1+f1*3^(1/2))*3^(1/2), 1/12*(f2-f2*3^(1/2)-f1+f1*3^(1/2))*3^(1/2), 1/2*f2-1/2*f1

>    CC:=seq(solve(subs(A=AA[k],B=BB[k],eqC),C),k=1..4);

CC := 0, 1/24*(-2*f2*f1*3^(1/2)-4*f2*f1+2*f1^2*3^(1/2)-f1^2+5*f2^2)*3^(1/2)/a1, 1/24*(-2*f2*f1*3^(1/2)+4*f2*f1+2*f1^2*3^(1/2)+f1^2-5*f2^2)*3^(1/2)/a1, -1/2/a1*f1*(-f2+f1)
CC := 0, 1/24*(-2*f2*f1*3^(1/2)-4*f2*f1+2*f1^2*3^(1/2)-f1^2+5*f2^2)*3^(1/2)/a1, 1/24*(-2*f2*f1*3^(1/2)+4*f2*f1+2*f1^2*3^(1/2)+f1^2-5*f2^2)*3^(1/2)/a1, -1/2/a1*f1*(-f2+f1)

So now I have found approximations of the kappa1-coordinates of all 4 critical points.

Proceed to determine the kappa2- and kappa3-coordinates! Use kap2 and kap3...

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

kap21 := 0

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

kap31 := 0

Put DD[1]=EE[1]=FF[1]=0:

>    DD[1]:=0; EE[1]:=0; FF[1]:=0;

DD[1] := 0

EE[1] := 0

FF[1] := 0

>    kap2;

-(-a1*N+kappa1)*kappa1

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

kap22 := 1/6*a1^2*N^2+1/12*(-f2-f2*3^(1/2)+f1+f1*3^(1/2))*a1*N

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

DD[2] := 1/6*a1^2

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

EE[2] := 1/12*(1+3^(1/2))*a1*(-f2+f1)

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

kap32 := 1/18*a1^3*3^(1/2)*N^3-1/12*a1^2*(-3*f1+f1*3^(1/2)-f2*3^(1/2)+3*f2)*N^2-1/24*a1*3^(1/2)*(-2*f2*f1+f2^2+f1^2)*N

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

FF[2] := 1/18*a1^3*3^(1/2)

Now work on the third critical point:

>    kap23:=map(simplify,convert(asympt(subs(kappa1=AA[3]*N+BB[3]+CC[3]/N,kap2),N,0),polynom));

kap23 := 1/6*a1^2*N^2-1/12*(f2-f2*3^(1/2)-f1+f1*3^(1/2))*a1*N

>    DD[3]:=coeff(kap23,N,2);

DD[3] := 1/6*a1^2

>    EE[3]:=factor(coeff(kap23,N,1));

EE[3] := -1/12*(-1+3^(1/2))*a1*(-f2+f1)

>    kap33:=map(simplify,convert(asympt(subs(kappa1=AA[3]*N+BB[3]+CC[3]/N,kap3),N,0),polynom));

kap33 := -1/18*a1^3*3^(1/2)*N^3

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

FF[3] := -1/18*a1^3*3^(1/2)

NOTE here that FF[3] is the coefficient of N^3!  Proceed to the last critical point.  This one has kappa1 sim K.

I expect different asymptotic orders for D,E,F.

>    kap24:=map(simplify,convert(asympt(subs(kappa1=AA[4]*N+BB[4]+CC[4]/N,kap2),N,1),polynom));

kap24 := 1/2*(-f2+f1)*a1*N-1/4*(-f2+f1)^2+1/2*f1*(-f2+f1)

>    DD[4]:=coeff(kap24,N,1);

DD[4] := 1/2*a1*(-f2+f1)

>    EE[4]:=factor(coeff(kap24,N,0));

EE[4] := 1/4*(-f2+f1)*(f1+f2)

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

kap34 := 1/2*N*a1*f2^2-1/2*N*f1*a1*f2+1/2*f1*f2^2-7/4*f2*f1^2+f1^3+1/4*f2^3

>    FF[4]:=factor(coeff(kap34,N,1));

FF[4] := -1/2*a1*f2*(-f2+f1)

So I have derived asymptotic approximations of the first three cumulants at all 4 critical points. Denote the critical points by P1, P2,P3,P4.

>    P1:=[0,0,0];

P1 := [0, 0, 0]

>    P2:=[AA[2]*N+BB[2]+CC[2]/N,DD[2]*N^2+EE[2]*N,FF[2]*N^3];

P2 := [(1/2-1/6*3^(1/2))*a1*N+1/12*(-f2-f2*3^(1/2)+f1+f1*3^(1/2))*3^(1/2)+1/24*(-2*f2*f1*3^(1/2)-4*f2*f1+2*f1^2*3^(1/2)-f1^2+5*f2^2)*3^(1/2)/a1/N, 1/6*a1^2*N^2+1/12*(1+3^(1/2))*a1*(-f2+f1)*N, 1/18*a1^3...
P2 := [(1/2-1/6*3^(1/2))*a1*N+1/12*(-f2-f2*3^(1/2)+f1+f1*3^(1/2))*3^(1/2)+1/24*(-2*f2*f1*3^(1/2)-4*f2*f1+2*f1^2*3^(1/2)-f1^2+5*f2^2)*3^(1/2)/a1/N, 1/6*a1^2*N^2+1/12*(1+3^(1/2))*a1*(-f2+f1)*N, 1/18*a1^3...
P2 := [(1/2-1/6*3^(1/2))*a1*N+1/12*(-f2-f2*3^(1/2)+f1+f1*3^(1/2))*3^(1/2)+1/24*(-2*f2*f1*3^(1/2)-4*f2*f1+2*f1^2*3^(1/2)-f1^2+5*f2^2)*3^(1/2)/a1/N, 1/6*a1^2*N^2+1/12*(1+3^(1/2))*a1*(-f2+f1)*N, 1/18*a1^3...

>    P3:=[AA[3]*N+BB[3]+CC[3]/N,DD[3]*N^2+EE[3]*N,FF[3]*N^3];

P3 := [(1/2+1/6*3^(1/2))*a1*N+1/12*(f2-f2*3^(1/2)-f1+f1*3^(1/2))*3^(1/2)+1/24*(-2*f2*f1*3^(1/2)+4*f2*f1+2*f1^2*3^(1/2)+f1^2-5*f2^2)*3^(1/2)/a1/N, 1/6*a1^2*N^2-1/12*(-1+3^(1/2))*a1*(-f2+f1)*N, -1/18*a1^...
P3 := [(1/2+1/6*3^(1/2))*a1*N+1/12*(f2-f2*3^(1/2)-f1+f1*3^(1/2))*3^(1/2)+1/24*(-2*f2*f1*3^(1/2)+4*f2*f1+2*f1^2*3^(1/2)+f1^2-5*f2^2)*3^(1/2)/a1/N, 1/6*a1^2*N^2-1/12*(-1+3^(1/2))*a1*(-f2+f1)*N, -1/18*a1^...
P3 := [(1/2+1/6*3^(1/2))*a1*N+1/12*(f2-f2*3^(1/2)-f1+f1*3^(1/2))*3^(1/2)+1/24*(-2*f2*f1*3^(1/2)+4*f2*f1+2*f1^2*3^(1/2)+f1^2-5*f2^2)*3^(1/2)/a1/N, 1/6*a1^2*N^2-1/12*(-1+3^(1/2))*a1*(-f2+f1)*N, -1/18*a1^...

>    P4:=[AA[4]*N+BB[4]+CC[4]/N,DD[4]*N+EE[4],FF[4]*N];

P4 := [a1*N+1/2*f2-1/2*f1-1/2/a1*f1*(-f2+f1)/N, 1/2*(-f2+f1)*a1*N+1/4*(-f2+f1)*(f1+f2), -1/2*a1*f2*(-f2+f1)*N]

I would also like to study the local stability of the 4 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];

F := [(a1*N-kappa1)*kappa1-kappa2, (f1*a1*N-f2*kappa1)*kappa1+(2*a1*N-f2)*kappa2-4*kappa1*kappa2-2*kappa3, (a1*N-kappa1)*kappa1-6*f2*kappa1*kappa2-6*kappa1*kappa3+(3*f1*a1*N-1)*kappa2-6*kappa2^2+3*(a1*...
F := [(a1*N-kappa1)*kappa1-kappa2, (f1*a1*N-f2*kappa1)*kappa1+(2*a1*N-f2)*kappa2-4*kappa1*kappa2-2*kappa3, (a1*N-kappa1)*kappa1-6*f2*kappa1*kappa2-6*kappa1*kappa3+(3*f1*a1*N-1)*kappa2-6*kappa2^2+3*(a1*...

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

J := Matrix(%id = 2337152)

Determine the Jacobian in each of the four critical points:

>    J1:=subs(kappa1=P1[1],kappa2=P1[2],kappa3=P1[3],J);

J1 := Matrix(%id = 2610860)

Recall AA:

>    AA;

0, (1/2-1/6*3^(1/2))*a1, (1/2+1/6*3^(1/2))*a1, a1

>    J2:=subs(kappa1=P2[1],kappa2=P2[2],kappa3=P2[3],J):

>    J3:=subs(kappa1=P3[1],kappa2=P3[2],kappa3=P3[3],J):

>    J4:=subs(kappa1=P4[1],kappa2=P4[2],kappa3=P4[3],J):

Determine the eigenvalues of the matrices J1-J4.

The computations are demanding. So I modify each of the matrices  J2-J4 by replacing each of its elements by the first (one or) two terms in its asymptotic approximation as N aproaches infinity.  

The modification of J2 is called J2a:

>    J2a:=J2:

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

J2a[1,1] := 1/3*a1*3^(1/2)*N-1/6*(3+3^(1/2))*(-f2+f1)

>    J2a[1,2];

-1

>    J2a[1,3];

0

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

J2a[2,1] := -2/3*a1^2*N^2-2/3*f2*a1*N+2/3*N*a1*f2*3^(1/2)+2/3*f1*a1*N-1/3*N*a1*f1*3^(1/2)

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

J2a[2,2] := 2/3*a1*3^(1/2)*N+1/3*f2*3^(1/2)-1/3*f1*3^(1/2)-f1

>    J2a[2,3];

-2

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

J2a[3,1] := -1/3*a1^3*3^(1/2)*N^3-f2*a1^2*N^2

>    J2a[3,2]:=simplify(convert(asympt(J2[3,2],N,0),polynom));

J2a[3,2] := -2*a1^2*N^2-2*f2*a1*N+2*N*a1*f2*3^(1/2)+2*f1*a1*N-N*a1*f1*3^(1/2)

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

J2a[3,3] := a1*3^(1/2)*N+1/2*f2*3^(1/2)-3/2*f2-1/2*f1*3^(1/2)-3/2*f1

The modification of J3 is called J3a:

>    J3a:=J3:

>    J3a[1,1]:=simplify(convert(asympt(J3[1,1],N,1),polynom));

J3a[1,1] := -1/3*a1*3^(1/2)*N-1/6*f2*3^(1/2)+1/2*f2+1/6*f1*3^(1/2)-1/2*f1

>    J3a[1,2];

-1

>    J3a[1,3];

0

>    J3a[2,1]:=simplify(convert(asympt(J3[2,1],N,0),polynom));

J3a[2,1] := -2/3*a1^2*N^2-2/3*f2*a1*N-2/3*N*a1*f2*3^(1/2)+2/3*f1*a1*N+1/3*N*a1*f1*3^(1/2)

>    J3a[2,2]:=simplify(convert(asympt(J3[2,2],N,1),polynom));

J3a[2,2] := -2/3*a1*3^(1/2)*N-1/3*f2*3^(1/2)+1/3*f1*3^(1/2)-f1

>    J3a[2,3];

-2

>    J3a[3,1]:=convert(asympt(J3[3,1],N,-1),polynom);

J3a[3,1] := 1/3*a1^3*3^(1/2)*N^3-f2*a1^2*N^2

>    J3a[3,2]:=convert(asympt(J3[3,2],N,0),polynom);

J3a[3,2] := -2*a1^2*N^2+((-1+3^(1/2))*a1*(-f2+f1)+3*f1*a1-6*f2*(1/2+1/6*3^(1/2))*a1)*N

>    J3a[3,3]:=map(factor,convert(asympt(J3[3,3],N,1),polynom));

J3a[3,3] := -a1*3^(1/2)*N+1/2*(3^(1/2)-3)*(-f2+f1)-3*f2

Modify J4:

>    J4a:=J4:

>    J4a[1,1]:=simplify(convert(asympt(J4[1,1],N,1),polynom));

J4a[1,1] := -a1*N-f2+f1

>    J4a[1,2];

-1

>    J4a[1,3];

0

>    J4a[2,1]:=simplify(convert(asympt(J4[2,1],N,1),polynom));

J4a[2,1] := -f1*a1*N+f2*f1-f1^2

>    J4a[2,2]:=simplify(convert(asympt(J4[2,2],N,1),polynom));

J4a[2,2] := -2*a1*N-3*f2+2*f1

>    J4a[2,3];

-2

>    J4a[3,1]:=simplify(convert(asympt(J4[3,1],N,1),polynom));

J4a[3,1] := -a1*N-f2+3/2*f2^3-3/2*f2*f1^2+f1

>    J4a[3,2]:=simplify(convert(asympt(J4[3,2],N,1),polynom));

J4a[3,2] := -3*f1*a1*N-3*f1^2-1+3*f2*f1

>    J4a[3,3]:=convert(asympt(J4[3,3],N,1),polynom);

J4a[3,3] := -3*a1*N-6*f2+3*f1

The eigenvalues of the 4 matrices:

>    eg1:=Eigenvalues(J1):

>    eg2:=Eigenvalues(J2a):

>    eg3:=Eigenvalues(J3a):

>    eg4:=Eigenvalues(J4a):

Determine asymptotic approximations of the eigenvalues!

>    a1:='a1':

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

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

eg11b := 3*a1*N-6*f1-3*f2

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

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

eg12b := a1*N+f1

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

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

eg13b := 2*a1*N-f2+5*f1

The three eigenvalues are asymptotically K, 2K, 3K. All of them positive. So the origin (P1) is unstable.

Next consider the eigenvalues of the matrix J2.

>    a1:='a1':

>    eg21a:=evalc(convert(asympt(eg2[1],N,3),polynom)):

>    assume(a1>0);
eg21b:=evalf(simplify(eg21a));

eg21b := -2.162237242*f1-1.110490029*f2-.3985797291e-9*I*f2-.9964493227e-10*I*f1+3.464101614*a1*N
eg21b := -2.162237242*f1-1.110490029*f2-.3985797291e-9*I*f2-.9964493227e-10*I*f1+3.464101614*a1*N

>    a1:='a1':
 eg22a:=evalc(convert(asympt(eg2[2],N,3),polynom)):

>    assume(a1>0);
eg22b:=evalf(simplify(eg22a));

eg22b := -1.530011547*f1+.6939976900*f2-.9999999993*a1*N+.2491123307e-9*I*f1-.1370117819e-8*I*f2
eg22b := -1.530011547*f1+.6939976900*f2-.9999999993*a1*N+.2491123307e-9*I*f1-.1370117819e-8*I*f2

>    a1:='a1':
eg23a:=evalc(convert(asympt(eg2[3],N,3),polynom)):

>    assume(a1>0);
eg23b:=evalf(simplify(eg23a));

eg23b := -1.039802015*f1+1.148543145*f2+1.000000000*a1*N+.9964493227e-9*I*f2-.4982246614e-9*I*f1
eg23b := -1.039802015*f1+1.148543145*f2+1.000000000*a1*N+.9964493227e-9*I*f2-.4982246614e-9*I*f1

The eigenvalues are, asymptotically, 3.46K, -K, K. Hence the critical point P2 is unstable.

Now look at the eigenvalues of J3.

>    a1:='a1':

>    eg31a:=evalc(convert(asympt(eg3[1],N,3),polynom)):

>    assume(a1>0);
eg31b:=evalf(map(simplify,eg31a));

eg31b := .9999999999*a1*N-2.142016166*f1+1.306002309*f2-.1228727676e-2*I*(.1e-6*f1-.3e-6*f2)
eg31b := .9999999999*a1*N-2.142016166*f1+1.306002309*f2-.1228727676e-2*I*(.1e-6*f1-.3e-6*f2)

>    a1:='a1':
eg32a:=evalc(convert(asympt(eg3[2],N,3),polynom)):

>    assume(a1>0);
eg32b:=evalf(map(simplify,eg32a));

eg32b := -3.464101615*a1*N+1.616782699*f1-4.889509970*f2-.6143638378e-3*I*(.1e-5*f1-.11e-5*f2)
eg32b := -3.464101615*a1*N+1.616782699*f1-4.889509970*f2-.6143638378e-3*I*(.1e-5*f1-.11e-5*f2)

>    a1:='a1':
eg33a:=evalc(convert(asympt(eg3[3],N,3),polynom)):

>    assume(a1>0);
eg33b:=evalf(map(simplify,eg33a));

eg33b := -1.000000001*a1*N-.7427157252*f1+.8514568559*f2+.6143638378e-3*I*(.11e-5*f1-.2e-6*f2)
eg33b := -1.000000001*a1*N-.7427157252*f1+.8514568559*f2+.6143638378e-3*I*(.11e-5*f1-.2e-6*f2)

So the eigenvalues of J3 are: K,  -3.46 K,  -K.  Unstable critical point!  

Now look at the eigenvalues of J4.

>    a1:='a1':
eg41a:=evalc(convert(asympt(eg4[1],N,3),polynom)):

>    assume(a1>0);
eg41b:=map(simplify,eg41a);

eg41b := -a1*N-f2+2*f1

>    a1:='a1':
eg42a:=evalc(convert(asympt(eg4[2],N,3),polynom)):

>    assume(a1>0);
eg42b:=map(simplify,eg42a);

eg42b := -3*a1*N-6*f2-3*f1

>    a1:='a1':
eg43a:=evalc(convert(asympt(eg4[3],N,3),polynom)):

>    assume(a1>0):
eg43b:=map(simplify,eg43a);

eg43b := -2*a1*N-3*f2+7*f1

The eigenvalues of J4 are asymptotically -K, -3K, -2K.  Thus P4 is stable.