C---------------------------------------------------------------------- C---------------------------------------------------------------------- C ECKHAUS.f : Determines the Eckhaus instability boundary for the C predator prey reaction-diffusion model C C Part of the Electronic Supplementary Material for C Sherratt, J.A. & Smith M.J. (2007) C Periodic Travelling Waves in Cyclic Populations: C Field studies and reaction-diffusion models C C This is the AUTO code used to continue in alpha from one point at which C stability changes through an Eckhaus instability, to determine the stability C boundary for a range of alpha. C This code is essentially a modification of spectrum.f but with equations C added to calculate the second derivative of the real part of the eigenvalue C with respect to gamma. C Eckhaus instabilities occur at the point where this second derivative is zero 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). We first of all converted these equations C into coupled first order ODEs (to give four equations), and then we split them 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. The final 8 equations calculate C the first and second derivative of the real part of the eigenvalue with C respect to gamma. The second derivative (PAR(8)) is indicative of wave stability C in these predator-prey equations. C C SUBROUTINE STPNT defines the parameter values involved in continuation. C In addition to the parameters defined in spectrum.f we define the C first and second derivatives of the real part of \lambda with respect to C \gamma at \lambda=0, and also define a dummy variable (PAR(9)) which we use C in the first run to allow AUTO to find a solution. 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 c=par(1) !wave speed alpha=par(2) !ratio of prey to predator dispersal rate period=par(3) !wavelength gamma=par(4) !varies between 0 and 2pi, see text for details evalr=par(5) !real part of the eigenvalue (lambda) evali=par(6) !imaginary part of lambda drel=par(7) !first derivative of Re(lambda) with respect to gamma ddrel=par(8) !second derivative of Re(lambda) with respect to gamma p1=0.2d0 !kappa p2=0.15d0 !sigma p3=0.05d0 !mu 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) u1rgv=u(13) u2rgv=u(14) v1rgv=u(15) v2rgv=u(16) u1rek=u(17) u2rek=u(18) v1rek=u(19) v2rek=u(20) fkinet=(1-u1ptw)*u1ptw-(u1ptw*v1ptw)/(u1ptw+p1) gkinet=(p2*u1ptw*v1ptw)/(u1ptw+p1)-p3*v1ptw Du=alpha**0.5 Dv=(1/alpha)**0.5 a11=1-2*u1ptw-(p1*v1ptw)/((u1ptw+p1)**2) a12=(-u1ptw)/(u1ptw+p1) a21=(p2*p1*v1ptw)/((u1ptw+p1)**2) a22=((p2*u1ptw)/(u1ptw+p1))-p3 f(1)=period*u2ptw f(2)=period*((1/Du)*(-c*u2ptw-fkinet)) f(3)=period*v2ptw f(4)=period*((1/Dv)*(-c*v2ptw-gkinet)) 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)) f(13)=period*(u2rgv-u1r) f(14)=period*((1/Du)*(drel*u1r- & a11*u1rgv-a12*v1rgv-c*u2rgv)-u2r) f(15)=period*(v2rgv-v1r) f(16)=period*((1/Dv)*(drel*v1r- & a21*u1rgv-a22*v1rgv-c*v2rgv)-v2r) f(17)=period*(u2rek-2*u1rgv) f(18)=period*((1/Du)*(ddrel*u1r+2*drel*u1rgv- & a11*u1rek-a12*v1rek-c*u2rek)-2*u2rgv) f(19)=period*(v2rek-2*v1rgv) f(20)=period*((1/Dv)*(ddrel*v1r+2*drel*v1rgv- & a21*u1rek-a22*v1rek-c*v2rek)-2*v2rgv) 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 C The initial data corresponding to these parameters is C read in from the file eckhaus.dat par(1)=0.496355 !wave speed par(2)=1.0d0 !ratio of prey to predator dispersal rate par(3)=21.092494146 !wavelength par(4)=0.0d0 !gamma par(5)=0.0d0 !Re lambda par(6)=0.0d0 !Im lambda par(7)=0.0d0 !dlambda/dgamma par(8)=0.0d0 !d2lambda/dgamma2 par(9)=0.0d0 !unspecified dummy C 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) u1gv0=u0(13) u2gv0=u0(14) v1gv0=u0(15) v2gv0=u0(16) u1ek0=u0(17) u2ek0=u0(18) v1ek0=u0(19) v2ek0=u0(20) 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) u1gv1=u1(13) u2gv1=u1(14) v1gv1=u1(15) v2gv1=u1(16) u1ek1=u1(17) u2ek1=u1(18) v1ek1=u1(19) v2ek1=u1(20) C gamma=par(4) cg=cos(gamma) sg=sin(gamma) C fb(1)=u1ptw0-u1ptw1 fb(2)=u2ptw0-u2ptw1 fb(3)=v1ptw0-v1ptw1 fb(4)=v2ptw0-v2ptw1 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 fb(13)=u1gv0-u1gv1 fb(14)=u2gv0-u2gv1 fb(15)=v1gv0-v1gv1 fb(16)=v2gv0-v2gv1 fb(17)=u1ek0-u1ek1 fb(18)=u2ek0-u2ek1 fb(19)=v1ek0-v1ek1 fb(20)=v2ek0-v2ek1 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) u1cg=u(13) u2cg=u(14) v1cg=u(15) v2cg=u(16) u1ek=u(17) u2ek=u(18) v1ek=u(19) v2ek=u(20) 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 fi(4)=u1r*u1cg+u2r*u2cg fi(5)=u1r*u1ek+u2r*u2ek C 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----------------------------------------------------------------------