>
CALCULATION OF NORMAL FORM FOR PREDATOR PREY SYSTEM NEAR HO PF BIFURCATION
THIS MAPLE WORKSHEET COMPLEMENTS THE PAPER "PERIODIC TRAVELLING WAVES IN
CYCLIC PREDATOR-PREY SYSTEMS" BY J.A. SHERRATT
INITIAL SETUPS FOR MAPLE
> restart;
> assume(A>1);
> assume(B>0);
> interface(showassumed= 0);
DEFINE THE EQUATIONS
> f:=q->q/(1+q);
q
f := q -> -----
1 + q
> fp :=P*(A*f(H*C)-1)/(B*A):
> fh:=H*(1-H)-P*f(H*C):
(these are the nondimensionalised population kinetics terms: as in equation (A.3) in the paper)
STEP 2 IN APPENDIX OF PAPER: STEADY STATES AND LINEAR STABILITY ANALYSIS
> Hs:=solve(A*f(H*C)=1,H);
1
Hs := ---------
C (A - 1)
> Ps:=A*Hs*(1-Hs);
/ 1 \
A |1 - ---------|
\ C (A - 1)/
Ps := -----------------
C (A - 1)
> tr:=1-2*Hs-Ps*C*subs(q=Hs*C,diff(f(q),q)):
> Ccrit:=solve(tr=0,C);
A + 1
Ccrit := -----
A - 1
> det:=C*Ps*subs(q=Hs*C,diff(f(q),q))/(B*A):
> real_evalue:=tr/2:
> imag_evalue:=sqrt(det-tr*tr/4):
> lambda0:=mu*subs(C=Ccrit,diff(real_evalue,C)):
(here mu=C-Ccrit)
> omega0:=subs(C=Ccrit,imag_evalue):
NOW SET C=Ccrit
> C:=Ccrit:
> simplify(real_evalue); #CHECK ON VALUE -- SHOULD BE ZERO
0
> simplify(imag_evalue); #CHECK ON VALUE
sqrt(A (A - 1) (A + 1) B)
-------------------------
A (A + 1) B
STEP 3 IN APPENDIX OF PAPER: CONVERT LINEAR PART TO NORMAL FORM
> P:=Ps-phat*sqrt(A*(1-2*Hs)/B):
> H:=hhat+Hs:
> curlyF:=-fp/sqrt(A*(1-2*Hs)/B)+omega0*hhat:
> curlyG:=fh-omega0*phat:
NOW RUN SOME CHECKS ON curlyF AND curlyG -- there should be no linear terms
> simplify(subs(phat=0,hhat=0,curlyF));
0
> simplify(subs(phat=0,hhat=0,curlyG));
0
> simplify(subs(phat=0,hhat=0,diff(curlyF,phat)));
0
> simplify(subs(phat=0,hhat=0,diff(curlyF,hhat)));
0
> simplify(subs(phat=0,hhat=0,diff(curlyG,phat)));
0
> simplify(subs(phat=0,hhat=0,diff(curlyG,hhat)));
0
IF ALL SIX OF THE ABOVE GIVE ZERO, WE CAN CONTINUE.
STEP 4: CALCULATE LAMBDA1 AND OMEGA1
(using the formula given in the paper)
> Fpp:= subs(phat=0,hhat=0,diff(curlyF,phat$2)):
> Fph:= subs(phat=0,hhat=0,diff(curlyF,phat,hhat)):
> Fhh:= subs(phat=0,hhat=0,diff(curlyF,hhat$2)):
> Gpp:= subs(phat=0,hhat=0,diff(curlyG,phat$2)):
> Gph:= subs(phat=0,hhat=0,diff(curlyG,phat,hhat)):
> Ghh:= subs(phat=0,hhat=0,diff(curlyG,hhat$2)):
> Fppp:=subs(phat=0,hhat=0,diff(curlyF,phat$3)):
> Fpph:=subs(phat=0,hhat=0,diff(curlyF,phat$2,hhat)):
> Fphh:=subs(phat=0,hhat=0,diff(curlyF,phat,hhat$2)):
> Fhhh:=subs(phat=0,hhat=0,diff(curlyF,hhat$3)):
> Gppp:=subs(phat=0,hhat=0,diff(curlyG,phat$3)):
> Gpph:=subs(phat=0,hhat=0,diff(curlyG,phat$2,hhat)):
> Gphh:=subs(phat=0,hhat=0,diff(curlyG,phat,hhat$2)):
> Ghhh:=subs(phat=0,hhat=0,diff(curlyG,hhat$3)):
> lambda1:=-(Fppp+Fphh+Gpph+Ghhh)/16+
> (Fpp*Gpp-Fhh*Ghh
> -Fph*(Fpp+Fhh)+Gph*(Gpp+Ghh))
> /(16*omega0):
> omega1:=(Gppp+Gphh-Fpph-Fhhh)/16+
> (Fpp*(Gph-Fhh)+Ghh*(Fph-Gpp)
> -3*Fhh*(Gph+Fhh)-3*Gpp*(Fph+Gpp)
> -2*(Fpp+Fhh-Gph)**2-2*(Gpp+Ghh-Fph)**2)
> /(48*omega0):
END OF CALCULATION: NOW SIMPLIFY THE RESULTS
> simplify(lambda0);
mu (A - 1)
1/2 ----------
A (A + 1)
> simplify(lambda1);
A + 1
1/4 -----
A
> simplify(omega0);
sqrt(A (A - 1) (A + 1) B)
-------------------------
A (A + 1) B
> simplify(omega1);
5 4 3 2 2 2
sqrt(A + 1) (A B + A + 4 A B - 2 A + 4 B A - 5 B A + 1)
- 1/24 -------------------------------------------------------------
(7/2) (3/2)
A sqrt(A - 1) B
NOW USE THE RESULTS TO TEST THE STABILITY OF PERIODIC TRAVELLING WAVES
GENERATED BY INVASION
> stability_test:=simplify((omega1/lambda1)**2);
> # Wavetrain is stable iff stability_test<1.1479
5 4 3 2 2 2 2
(A B + A + 4 A B - 2 A + 4 B A - 5 B A + 1)
stability_test := 1/36 --------------------------------------------------
2 5 3
(A - 1) A B
NOW VISUALISE THE STABLE/UNSTABLE PARAMETER REGIONS
(stable means stable periodic waves generated behind invasion;
unstable means spatiotemporal chaos generated behind invasion).
> with(plots):
> istable:=0: iunstable:=0:
> Amin:=1.01: Amax:=6.9: Bmin:=1.01: Bmax:=10:
> npoints:=30: Astep:=(Amax-Amin)/(npoints-1): Bstep:=(Bmax-Bmin)/(npoints-1):
> for aa from Amin to Amax by Astep do
> for bb from Bmin to Bmax by Bstep do
> if eval(subs(A=aa,B=bb,stability_test))<1.1479 then
> istable:=istable+1:
> Astable[istable]:=aa:
> Bstable[istable]:=bb
> else iunstable:=iunstable+1:
> Aunstable[iunstable]:=aa:
> Bunstable[iunstable]:=bb
> fi
> od
> od:
> stablepoints:={seq([Astable[i],Bstable[i]],i=1..istable)}:
> unstablepoints:={seq([Aunstable[i],Bunstable[i]],i=1..iunstable)}:
> stableplot:=pointplot(stablepoints,color=red,symbol=circle,thickness=2):
> unstableplot:=pointplot(unstablepoints,color=blue,symbol=box,thickness=2):
>
neutrallystableplot:=implicitplot(stability_test=1.1479,A=1.01..6.9,B=1.01..10,
color=green,thickness=3,numpoints=1000):
> unassign('ps');
> finalplot1:=display({stableplot,unstableplot,neutrallystableplot},
> labels=["A","B"],title="CIRCLE=stable, BOX=unstable"):
> plotsetup(default):
> finalplot1;
NOW INVESTIGATE AMPLITUDE AND SPEED OF PERIODIC TRAVELLING WAVES
STARTING WITH THE RESULTS FOR THE DIMENSIONLESS EQUATIONS
> rsq:=2*lambda0*(sqrt((lambda1**2+omega1**2))-lambda1)/(omega1**2):
(here rsq is the square of the amplitude of the periodic waves behind invasion)
> speed:=(omega0+omega1*rsq)/sqrt(lambda0-lambda1*rsq):
(here speed is the speed of the periodic waves behind invasion)
> dimlessplot:=array(1..2):
> densityplot((subs(mu=1.0, rsq)), A=1..7, B=1..10, grid=[30,30], style=PATCHNOGRID, colorstyle=HUE, axes=BOXED, labels=["A","B"], scaling=UNCONSTRAINED, title="AMPLITUDE OF PERIODIC WAVES");
> densityplot((subs(mu=1.0, speed)), A=1..7, B=1..10, grid=[30,30], style=PATCHNOGRID, colorstyle=HUE, axes=BOXED, labels=["A","B"], scaling=UNCONSTRAINED, title="SPEED OF PERIODIC WAVES");
NOW CONSIDER DIMENSIONAL VALUES, BASED ON THE PARAMETERS IN NISBET ET AL
(theor pop biol 40, 125-147, 1991) FOR A PHYTOPLANKTON--ZOOPLANKTON INTERACTION
> preymin:=h0*(-sqrt(rsq)+Hs):
> preymax:=h0*(sqrt(rsq)+Hs):
> preymin_plankton:=subs(mu=0.1*C, h0=C/6.1, A=0.5/b, B=0.5, preymin):
> preymax_plankton:=subs(mu=0.1*C, h0=C/6.1, A=0.5/b, B=0.5, preymax):
> minplot:=plot(preymin_plankton, b=0.1..0.45):
> maxplot:=plot(preymax_plankton, b=0.1..0.45):
> finalplot2:=display({minplot,maxplot},labels=["b ( /day)"," "], title="Prey min and max (mg/l) vs b"):
> plotsetup(default):
> finalplot2;
> dimensionalspeed:=speed*sqrt(delta*rparam):
> speed_plankton:=subs(mu=0.1*C, A=0.5/0.14, B=rparam/0.5, delta=8.6, dimensionalspeed):
> finalplot3:=plot(speed_plankton, rparam=0.3..2.0, labels=["r ( /day)"," "], title="Wave speed (cm/day) vs r"):
> plotsetup(default):
> finalplot3;
>
THE END