Instructions on calculating the stability of waves along the wave family using AUTO

Part of the Electronic Supplementary Material for Sherratt, J.A. & Smith M.J. (2008) Periodic Travelling Waves in Cyclic Populations: Field studies and reaction-diffusion models. J. R. Soc. Interface 5, 483-505

The techniques in the previous example allowed us to determine that wave stability must change at least twice along the wave family in Fig. 2(a), when α=1, between c=0.5 and c=5.0. Furthermore, general theory implies that sufficiently fast waves are also stable. Therefore wave stability changes at least three three times along the wave family to equations (A4), when α=1. In this example we will demonstrate how to accurately calculate these stability transition points and track how the wave properties at these points varies with α. The files you need for this example are eckhaus.f, r.eckhaus, r.eckhaus.2, r.eckhaus.3, and eckhaus.dat.

The techniques used in this example were adapted from Rademacher et al. (2007)

For the specific model studied here we always find that the transition between stable and unstable waves is associated with the same type of change to the stability spectrum; an Eckhaus instability. This type of stability transition occurs when the stability spectrum close to the origin moves into our out of the right hand half of the complex plane. More precisely, the second derivative of the real part of λ with respect to γ at λ=0 changes sign. The technique we use can only detect stability transitions of this type. The alternative possibility of a Hopf bifurcation of the periodic travelling wave does occur for some other reaction-diffusion systems (Bordiougov and Engel, 2006). However, for this system, we find that the sign of the second derivative is always indicative of whether waves are stable or not. Our approach will therefore to be to redraw the travelling wave family by continuing in wave speed from a wave speed close to the Hopf bifurcation point (which is always unstable) and monitor wave stability throughout the continuation. We then get AUTO to label solutions where the stability transition occurs. From the previous examples we expect at least two points in which the stability changes.

PART 1: Initialise the solution

The first step is to initialise AUTO with a periodic solution. We have chosen one that is close to the Hopf bifurcation point. Using identical techniques to the previous section we generate an initial data file. This has the periodic solution to equations (A4) for a given wave speed (first 4 equations), then 8 equations corresponding to the real and imaginary parts of eigenfunction equations (A6), and then 8 more equations corresponding to the first and second derivatives of the real part of the eigenfunction equations with respect to γ (equations (A8)). In this special case of λ=0 the respective complex parts of these equations is zero (see Rademacher et al. 2007 for an explanation). These last 8 equations are all initialised with zero values. To obtain the precise initial values for these equations it is sufficient to perform an initial AUTO continuation in a dummy parameter and AUTO obtains their initial values through Newton's method (see Rademacher et al. 2007 for details). The necessary data is stored in the file eckhaus.dat and the equations file is eckhaus.f. To read in the solution type

@fc eckhaus

Then to perform the initial contnuation, type

@r eckhaus dat

Using a lowercase "r" is important here. This gives the output

  BR    PT  TY LAB    PAR(9)        L2-NORM       MAX U(1)      PAR(5)        PAR(6)        PAR(7)        PAR(8)        PAR(1)
   1     3  EP   2  3.897930E+01  1.046080E+00  2.042113E-01 -9.940254E-16  0.000000E+00  7.728872E-01 -3.896163E+01  4.963551E-01

which shows that AUTO continued in the dummy parameter (PAR(9)) for 3 iterations. Note that at the end of this continuation the value of λ (PAR(6)) is still zero but the value of the first (PAR(7)) and second (PAR(8)) derivatives of λ with respect to γ are not. Note also that we have started with a solution that has a wave speed close to the bifurcation point (PAR(1)=0.5) and, as expected, the second derivative of Re(λ) with respect to γ (PAR(8)) is negative, meaning that the stability spectrum at the zero eigenvalue extends into the right hand complex plane; this corresponds to an unstable periodic travelling wave solution. We need to save this starting point by typing

@sv eckhaus

PART 2: Continue in wave speed and monitor wave stability

We next perform a continuation in wave speed and monitor wave stability (PAR(8)) by typing

@R eckhaus 2

which gives the output

  BR    PT  TY LAB    PAR(1)        L2-NORM       PAR(5)        PAR(6)        PAR(7)        PAR(8)        PAR(3)        PAR(2)
   1    57  UZ   3  9.376881E-01  1.095253E+00  5.679339E-17  0.000000E+00  6.077955E-01  1.232157E-07  4.272838E+01  1.000000E+00
   1   112  UZ   4  1.467805E+00  1.123403E+00 -5.653356E-12  0.000000E+00  8.676132E-01 -1.127586E-07  9.640921E+01  1.000000E+00
   1   585  UZ   5  6.700426E+00  1.137477E+00  1.438813E-09  0.000000E+00  6.513091E+00  1.164886E-06  5.683096E+02  1.000000E+00
   1   786  EP   6  9.011210E+00  1.137812E+00  6.657638E-10  0.000000E+00  8.867451E+00  7.822464E-01  7.692810E+02  1.000000E+00

This shows that wave speed (PAR(1)) was increased from about 0.5 to 9, and AUTO detected a stability transition (PAR(8)=0) at speeds of 0.94, 1.47 and 6.70. You can save this output as a datafile by typing

@sv eckhaus.fig1

which saves the full dataset in p.eckhaus.fig1.

PART 3: Continue in α at the stability transition point

Extensive calculation of wave stability using the techniques in the previous example indicated that for all values of α with this predator-prey model, wave stability always changes through an Eckhaus instability. Rather than laboriously calculate wave stability along selected wave families for fixed α we instead perform continuations in α from the stability transition points (PAR(8)=0). In this example we will continue from labelled solution number 4 above, reducing α and monitoring how the wave speed (PAR(1)) and wavelength (PAR(3)) of those labelled solutions vary. To execute this run, first of all type

@ap eckhaus

which allows AUTO to detect the previous solution labels, then type

@R eckhaus 3

which gives the output

  BR    PT  TY LAB    PAR(2)        L2-NORM       PAR(1)        PAR(5)        PAR(6)        PAR(7)        PAR(3)        PAR(8)
   1   442  LP   7  2.638544E-01  1.097144E+00  9.190852E-01 -4.448141E-14  0.000000E+00  5.956973E-01  5.266620E+01 -1.127586E-07
   1  1861  EP   8  1.000091E+02  1.100652E+00  2.670358E+00  1.831696E-16  0.000000E+00  1.769976E+00  1.149577E+02 -1.127586E-07

This output shows that AUTO ended the continuation at α=100, a limit defined in the constants file. Since we started this continuation by reducing α, AUTO must have found some turning point at which α started to increase again. AUTO identified this point (a "fold") at iteration 442, which is the first output line above. This means that AUTO initially detected the stability boundary for reductions in α (PAR(2)) but then detected that the line of periodic solutions, turned, and continued for increasing α. If this output is saved using the command

@sv eckhaus.fig2

and wave speed (PAR(1)) is plotted against α (PAR(2)) from the file p.eckhaus.fig2 then this turning point is clearly visible.

PREVIOUS EXAMPLE,

INDEX