C---------------------------------------------------------------------- C---------------------------------------------------------------------- C SPECTRUM.f : Stability spectrum calculation for the C predator prey reaction-diffusion model C C C Part of the Electronic Supplementary Material for C Sherratt, J.A. & Smith M.J. (2008) C Periodic Travelling Waves in Cyclic Populations: C Field studies and reaction-diffusion models C J. R. Soc. Interface 5, 483-505. C C C This is the AUTO code used to determine wave stability at given C positions along the wave family in Figure 2(a). C Instructions for running this code are given in the online appendix. C C SUBROUTINE FUNC again defines the four ODEs that have travelling wave C solutions; these are equations A4 in the Electronic Supplementary Material C Note that in C this case we have used z=x-ct as the travelling wave coordinate. C These equations are followed C by the eigenfunction equations (A6). C We first of all converted these equations C into coupled first order ODEs (to give four equations, and then C we split these C into their real and imaginary parts (since AUTO cannot directly handle C complex numbers). This means that the two eigenfunction equations (A6) are C actually represented by eight first order ODEs. C C Note that the wave speed and wavelength need to be accurately defined for C the given solution whose stability you want to calculate. C C SUBROUTINE STPNT defines the parameter values involved in continuation, C in this case, gamma and the real and imaginary parts of the eigenvalue. C Note that in this example we will be performing the continuation from the C zero eigenvalue, which is always part of the stability spectrum C of periodic travelling waves (corresponding to translation of the wave). C C SUBROUTINE BCND defines the boundary conditions for the problem. See the C main text for details. C C SUBROUTINE ICND defines the integral constraints for the problem. See C Rademacher et al (2007) for details. C C---------------------------------------------------------------------- C---------------------------------------------------------------------- C SUBROUTINE FUNC(NDIM,U,ICP,PAR,IJAC,F,DFDU,DFDP) C ---------- ---- C C Evaluates the algebraic equations or ODE right hand side C C Input arguments : C NDIM : Dimension of the ODE system C U : State variables C ICP : Array indicating the free parameter(s) C PAR : Equation parameters C C Values to be returned : C F : ODE right hand side values C C Normally unused Jacobian arguments : IJAC, DFDU, DFDP (see manual) C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION U(NDIM), PAR(*), F(NDIM), ICP(*) C C Define the parameters psigma=0.15d0 !prey to predator conversion rate pmu=0.05d0 !predator death rate pkappa=0.2d0 !half-saturation constant in prey consumption by predators c=par(1) !wave speed alpha=par(2) !ratio of predator to prey dispersal rate period=par(3) !wavelength gamma=par(4) !normally varied between 0 and 2pi, see the online appendix for why we vary this parameter evalr=par(5) !the real part of the eigenvalue evali=par(6) !the imaginary part of the eigenvalue C u1ptw=u(1) u2ptw=u(2) v1ptw=u(3) v2ptw=u(4) u1r=u(5) u1i=u(6) u2r=u(7) u2i=u(8) v1r=u(9) v1i=u(10) v2r=u(11) v2i=u(12) fkinet=(1-u1ptw)*u1ptw-(u1ptw*v1ptw)/(u1ptw+pkappa) !the prey kinetic equation gkinet=(psigma*u1ptw*v1ptw)/(u1ptw+pkappa)-pmu*v1ptw !the predator kinetic equation Du=alpha**0.5 !prey dispersal rate Dv=(1/alpha)**0.5 !predator dispersal rate C These are the first derivatives of the kinetic equations with respect to predator and prey density a11=1-2*u1ptw-(pkappa*v1ptw)/((u1ptw+pkappa)**2) a12=(-u1ptw)/(u1ptw+pkappa) a21=(psigma*pkappa*v1ptw)/((u1ptw+pkappa)**2) a22=((psigma*u1ptw)/(u1ptw+pkappa))-pmu C C These first four equations accurately calculate the periodic travelling wave solutions to () f(1)=period*u2ptw f(2)=period*((1/Du)*(-c*u2ptw-fkinet)) f(3)=period*v2ptw f(4)=period*((1/Dv)*(-c*v2ptw-gkinet)) C These eight equations calculate the eigenfunction equations () f(5)=period*u2r f(6)=period*u2i f(7)=period*((1/Du)*(evalr*u1r-evali*u1i & -(a11*u1r+a12*v1r)-c*u2r)) f(8)=period*((1/Du)*(evali*u1r+evalr*u1i & -(a11*u1i+a12*v1i)-c*u2i)) f(9)=period*v2r f(10)=period*v2i f(11)=period*((1/Dv)*(evalr*v1r-evali*v1i & -(a21*u1r+a22*v1r)-c*v2r)) f(12)=period*((1/Dv)*(evali*v1r+evalr*v1i & -(a21*u1i+a22*v1i)-c*v2i)) C RETURN END C---------------------------------------------------------------------- C---------------------------------------------------------------------- SUBROUTINE STPNT(NDIM,U,PAR) C ---------- ----- IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION U(NDIM), PAR(*) C C Initialize the equation parameters C par(2)=1.0d0 !ratio of predator to prey dispersal rate par(4)=0.0d0 !gamma par(5)=0.0d0 !real part of eigenvalue par(6)=0.0d0 !imaginary part of the eigenvalue C C We don't need to define initial solutions since they are read in C as datafiles. C C Here are three different wave speed and wavelength combinations for use C as initial data C COMBINATION 1 - COPY SPECTRUM.DAT.1 TO SPECTRUM.DAT FOR INITIAL DATA par(1)=0.5d0 !wave speed par(3)=2.1192137072E+01 !wavelength C COMBINATION 2 - COPY SPECTRUM.DAT.2 TO SPECTRUM.DAT FOR INITIAL DATA C par(1)=1.0d0 !wave speed C par(3)= 4.7320624885E+01!wavelength C COMBINATION 3 - COPY SPECTRUM.DAT.3 TO SPECTRUM.DAT FOR INITIAL DATA C par(1)=5.0d0 !wave speed C par(3)= 4.1933418094E+02 !wavelength RETURN END C---------------------------------------------------------------------- C---------------------------------------------------------------------- SUBROUTINE BCND(NDIM,PAR,ICP,NBC,U0,U1,FB,IJAC,DBC) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION PAR(*),ICP(*),U0(NDIM),U1(NDIM),FB(NBC) C u1ptw0=u0(1) u2ptw0=u0(2) v1ptw0=u0(3) v2ptw0=u0(4) u1r0=u0(5) u1i0=u0(6) u2r0=u0(7) u2i0=u0(8) v1r0=u0(9) v1i0=u0(10) v2r0=u0(11) v2i0=u0(12) C u1ptw1=u1(1) u2ptw1=u1(2) v1ptw1=u1(3) v2ptw1=u1(4) u1r1=u1(5) u1i1=u1(6) u2r1=u1(7) u2i1=u1(8) v1r1=u1(9) v1i1=u1(10) v2r1=u1(11) v2i1=u1(12) C gamma=par(4) cg=cos(gamma) sg=sin(gamma) C C C THESE EQUATIONS ENSURE THAT THE LEFT AND RIGHT HAND BOUNDARIES C OF THE TRAVELLING WAVE SOLUTION ARE THE SAME C fb(1)=u1ptw0-u1ptw1 fb(2)=u2ptw0-u2ptw1 fb(3)=v1ptw0-v1ptw1 fb(4)=v2ptw0-v2ptw1 C C THESE EQUATIONS ENSURE THAT THE LEFT HAND BOUNDARY C OF THE EIGENFUNCTION CALCULATION MATCHES THE RIGHT HAND BOUNDARY C WHEN MODIFIED BY SOME PHASE CHANGE C fb(5)=-u1r1+u1r0*cg-u1i0*sg fb(6)=-u1i1+u1i0*cg+u1r0*sg fb(7)=-u2r1+u2r0*cg-u2i0*sg fb(8)=-u2i1+u2i0*cg+u2r0*sg fb(9)=-v1r1+v1r0*cg-v1i0*sg fb(10)=-v1i1+v1i0*cg+v1r0*sg fb(11)=-v2r1+v2r0*cg-v2i0*sg fb(12)=-v2i1+v2i0*cg+v2r0*sg C RETURN END C---------------------------------------------------------------------- C---------------------------------------------------------------------- SUBROUTINE ICND(NDIM,PAR,ICP,NINT,U,UOLD,UDOT,UPOLD,FI,IJAC,DINT) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION U(NDIM),UOLD(NDIM),UDOT(NDIM),UPOLD(NDIM) DIMENSION FI(NINT),DINT(NINT,*),ICP(*),PAR(*) C fi(1)=0.0d0 do 100 j=1,4 fi(1)=fi(1)+upold(j)*(uold(j)-u(j)) 100 continue C u1r=u(5) u1i=u(6) u2r=u(7) u2i=u(8) v1r=u(9) v1i=u(10) v2r=u(11) v2i=u(12) u1rold=uold(5) u1iold=uold(6) u2rold=uold(7) u2iold=uold(8) v1rold=uold(9) v1iold=uold(10) v2rold=uold(11) v2iold=uold(12) C fi(2)=-1.0d0 & +u1r*u1r+u1i*u1i & +u2r*u2r+u2i*u2i & +v1r*v1r+v1i*v1i & +v2r*v2r+v2i*v2i fi(3)=u1rold*u1i-u1iold*u1r & +u2rold*u2i-u2iold*u2r & +v1rold*v1i-v1iold*v1r & +v2rold*v2i-v2iold*v2r C RETURN END C---------------------------------------------------------------------- C---------------------------------------------------------------------- C The following subroutines are not used here, C but they must be supplied as dummy routines C C SUBROUTINE FOPT RETURN END C SUBROUTINE PVLS RETURN END C---------------------------------------------------------------------- C----------------------------------------------------------------------