normal3.mws

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

Written with Maple8, on 030224.  I use the new linear algebra package LinearAlgebra.   

>    restart;

>    with(LinearAlgebra,Eigenvalues);

>    with(VectorCalculus,Jacobian);

[Eigenvalues]

[Jacobian]

>    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

In the normal distribution,  kappa4=0.  Inserting this value for kappa4 into eq3prel leads to moment closure and the rhs eq3:

>    eq3:=subs(kappa4=0,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

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 funciton 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*(-2*f2+7*a1*N+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
eq3a := 6*kappa1^4-12*a1*N*kappa1^3+a1*N*(-2*f2+7*a1*N+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

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.)  One root is equal to 0. The corresponding values for kappa2 and kappa3 are also 0:

>    kap2a:=subs(kappa1=0,kap2);

kap2a := 0

>    kap3a:=subs(kappa1=0,kap3);

kap3a := 0

The remaining 3 roots of eq3b satisfy the following third-degree equation:

>    eq3b:=collect(eq3a/kappa1,kappa1);

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

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

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

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

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

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

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

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

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

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

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

>    AA:=solve(eqA,A);

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

>    evalf(%);

a1, .7886751347*a1, .2113248653*a1

The A values of the remaining 3 critical points are given by AA[1]-AA[3]. I expect AA[2], AA[3] to correspond to spurious solutions.  

Determine asymptotic approx of the 3 coordinates of each of the remaining 3 critical points:

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

eqB1 := a1^2*B+a1^2*(-2*f2+2*f1)-1/2*a1*(-3*f2*a1+3*f1*a1)

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

BB[1] := 1/2*f2-1/2*f1

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

eqC1 := a1^2*C-1/2*a1*f2*f1+1/2*a1*f1^2

>    CC[1]:=solve(eqC1,C);

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

An asymptotic approximation for the first non-zero solution for kappa1 has thus been found. It is

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+EE[1]  and  kappa3 sim FF[1]*N.

>    kap2;

(a1*N-kappa1)*kappa1

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

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

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

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

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

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

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

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

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

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

Now proceed to find approximations of the remaining 2 critical points, with subscripts 2 and 3.

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

eqB2 := a1^2*B-a1^2*B*3^(1/2)+1/2*a1^2*f2-1/3*a1^2*f2*3^(1/2)-1/2*a1^2*f1+1/3*a1^2*f1*3^(1/2)

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

BB[2] := 1/6*(-3*f1+2*f1*3^(1/2)+3*f2-2*f2*3^(1/2))/(-1+3^(1/2))

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

eqC2 := 1/12*a1*(24*C*a1*3^(1/2)+9*f1^2-48*C*a1-24*f2*f1+15*f2^2+14*f1*3^(1/2)*f2-10*f2^2*3^(1/2)-4*f1^2*3^(1/2))/(-1+3^(1/2))

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

CC[2] := 1/24*(-9*f1^2+24*f2*f1-15*f2^2-14*f1*3^(1/2)*f2+10*f2^2*3^(1/2)+4*f1^2*3^(1/2))/a1/(3^(1/2)-2)

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

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

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

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

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

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

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

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

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

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

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

And kappa3 = FF[2]*N^3.  Now work on the last of the critical points!

>    eqB3:=simplify(subs(A=AA[3],eqB));

eqB3 := a1^2*B+a1^2*B*3^(1/2)+1/2*a1^2*f2+1/3*a1^2*f2*3^(1/2)-1/2*a1^2*f1-1/3*a1^2*f1*3^(1/2)

>    BB[3]:=solve(eqB3,B);

BB[3] := 1/6*(3*f1+2*f1*3^(1/2)-3*f2-2*f2*3^(1/2))/(1+3^(1/2))

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

eqC3 := 1/12*a1*(24*C*a1*3^(1/2)-9*f1^2+48*C*a1+24*f2*f1-15*f2^2+14*f1*3^(1/2)*f2-10*f2^2*3^(1/2)-4*f1^2*3^(1/2))/(1+3^(1/2))

>    CC[3]:=solve(eqC3,C);

CC[3] := 1/24*(15*f2^2+9*f1^2+10*f2^2*3^(1/2)-24*f2*f1-14*f1*3^(1/2)*f2+4*f1^2*3^(1/2))/a1/(3^(1/2)+2)

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

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

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

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

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

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

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

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

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

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

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

So I have derived asymptotic approximations of the first three cumulants...

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 = 3599480)

Determine the Jacobian in each of the four critical points:

>    J0:=subs(kappa1=0,kappa2=0,kappa3=0,J);

J0 := Matrix(%id = 3603600)

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

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

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

Determine the eigenvalues of the matrix J0. Expect instability...

>    a1:='a1':
eg0:=Eigenvalues(J0):

>    a1:='a1':

>    eg01a:=convert(asympt(eg0[1],N,3),polynom):

>    assume(a1>0);
eg01b:=simplify(eg01a);

eg01b := -6*f1-3*f2+3*a1*N

>    a1:='a1':
eg02a:=convert(asympt(eg0[2],N,2),polynom):

>    assume(a1>0);
eg02b:=simplify(eg02a);

eg02b := a1*N+f1

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

>    assume(a1>0);
eg03b:=simplify(eg03a);

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

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

Next consider the eigenvalues of the matrix J1. The evaluations are computationally demanding. Therefore, I determine a modified matrix J1a, whose entries contain asymptotic approximations of the corresponding entries of J1.

>    a1:='a1';

a1 := 'a1'

>    J1a:=J1:

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

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

>    J1a[1,2];

-1

>    J1a[1,3];

0

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

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

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

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

>    J1a[2,3];

-2

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

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

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

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

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

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

The eigenvalues of J1a:

>    a1:='a1':
 eg1:=Eigenvalues(J1a):

Asymptotic approximations of the eigenvalues:

>    a1:='a1':

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

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

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

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

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

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

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

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

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

>    a1:='a1';

a1 := 'a1'

The eigenvalues are, asymptotically, -K, -2K, -3K. Hence the critical point P1 is as stable.

Now look at the eigenvalues of J2. Determine the modified matrix J2a, as above.

>    J2a:=J2:

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

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

>    J2a[1,2];

-1

>    J2a[1,3];

0

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

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

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

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

>    J2a[2,3];

-2

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

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

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

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

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

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

The eigenvalues of J2a:

>    a1:='a1':
 eg2:=Eigenvalues(J2a):

Asymptotic approximations of the eigenvalues:

>    a1:='a1':

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

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

eg21b := 1.000000000*a1*N-2.142016166*f1+1.306002310*f2-.4e-9*I*f1+.1e-9*I*f2

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

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

eg22b := -3.464101615*a1*N+.7473369920e-9*I*f2-.1494673984e-9*I*f1+1.616782698*f1-4.889509972*f2
eg22b := -3.464101615*a1*N+.7473369920e-9*I*f2-.1494673984e-9*I*f1+1.616782698*f1-4.889509972*f2

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

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

eg23b := .2491123307e-9*I*f1-.3736684960e-9*I*f2-.7427157252*f1+.8514568540*f2-.9999999993*N*a1
eg23b := .2491123307e-9*I*f1-.3736684960e-9*I*f2-.7427157252*f1+.8514568540*f2-.9999999993*N*a1

>    a1:='a1';

a1 := 'a1'

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

Now look at the eigenvalues of J3. Again, determine the modified matrix J3a:

>    J3a:=J3:

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

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

>    J3a[1,2];

-1

>    J3a[1,3];

0

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

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

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

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

>    J3a[2,3];

-2

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

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

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

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

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

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

The eigenvalues of J3a are:

>    a1:='a1':
eg3:=Eigenvalues(J3a):

Asymptotic approximaitons of the eigenvalues:

>    a1:='a1':

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

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

eg31b := 3.464101614*a1*N-.4982246614e-9*I*f2+.1992898645e-9*I*f1-2.162237242*f1-1.110490029*f2
eg31b := 3.464101614*a1*N-.4982246614e-9*I*f2+.1992898645e-9*I*f1-2.162237242*f1-1.110490029*f2

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

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

eg32b := -.9999999988*a1*N-1.530011548*f1+.6939976917*f2-.74e-9*I*f2-.36e-10*I*f1

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

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

eg33b := 1.000000000*a1*N-.2491123307e-9*I*f1+.9964493227e-9*I*f2-1.039802016*f1+1.148543145*f2
eg33b := 1.000000000*a1*N-.2491123307e-9*I*f1+.9964493227e-9*I*f2-1.039802016*f1+1.148543145*f2

>   

The eigenvalues of J3 are asymptotically 3.46K, -K, K.  Thus P3 is unstable.