Instructions on calculating the stability spectrum of selected waves 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

Below we give instructions and AUTO code for calculating the stability of periodic travelling wave solutions. Those wishing to modify this code to suit their own system of interest should be aware that it requires more background work than is required for the previous section. However, the files supplied should allow the user to calculate the stability of three members of the wave familydrawn in the previous example. The files you need for this example are spectrum.f, r.spectrum, r.spectrum.2, r.spectrum.3, spectrum.dat.1, spectrum.dat.2, and spectrum.dat.3.

For illustrative purposes we selected three travelling wave solutions from the wave family corresponding to α=1 in Fig.2(b). These had wave speeds of 0.5, 1.0 and 5.0. From the figure it can be seen that the solution with c=1.0 is stable and the others are not. The instructions illustrate how this is determined.

The equations file required for these calculations is spectrum.f. To avoid problems associated with numercal inaccuracies we get AUTO to calculate the periodic travelling wave solutions at the same time as the periodic eigenfunctions. Therefore the first four equations are the travelling wave ODEs associated with equations (2) (equations A4). They are similar to those used in the preceeding example except we have used z=x-ct as the travelling wave coordinate for convenience. To represent the eigenfunction equations (A6) we need to convert them into a system of first order ODEs in the standard way. Since the eigenfunctions and eigenvalues can be complex we also need to calculate the equations seperately for the real and imaginary parts. Hence equations (A6) become a system of 8 first-order ODEs. Therefore in total we solve a system of 12 first-order ODEs.

PART 1: Calculate the initial solutions

To calculate the periodic eigenfunctions (A6) requires initial periodic solutions. These can be calculated in several ways. In general we utilise a standard technique to obtain eigenvalues to equation (A6) when γ=0. Mathematically speaking, the technique involves discretising equations (A6) using finite differences and using programs such as MATLAB or LAPACK to calculate some of the eigenvalues and corresponding eigenfunctions. Since this is a relatively standard technique, we omit details; however MATLAB code to do this is available from the authors on request. For the purposes of this example we just consider continuation in γ starting at one of the eignvalues for γ=0, namely λ=0. This is always an eigenvalue due to translation invariance. We have calculated the eigenfunctions when λ=0 for c=0.5, 1.0 and 5.0, and these are saved as spectrum.dat.1, spectrum.dat.2, and spectrum.dat.3, respectively. These data files contain one column with the normalised spatial coordinate, followed by four columns with the values corresponding to the periodic travelling wave solution to (A4) and finally eight columns with the values corresponding to the eigenfunctions for when to λ=0 in the system of equations (A6), split into real and imaginary parts.

PART 2: Continue in γ to draw the eigenvalue spectrum

The first step is to get AUTO to use the data files as initial data. First copy the data file of your choice to be named "spectrum.dat". Then edit the spectrum.f file to ensure that the wave speed (par(1)) and wavelength (par(3)) declared in STPT are those corresponding to the data file, as detailed in spectrum.f. You can uncomment these lines by removing the "C" at the start of the relevant line and comment them by adding "C". Once this is done, type

@fc spectrum

and this should create the file q.dat from spectrum.dat, in the format required by AUTO.

We next draw the stability spectrum by getting AUTO to do a continuation in γ (par(4)) while tracking the real (par(5)) and imaginary (par(6)) parts of λ. We first do two pre-continuation steps to ensure that the periodic solutions used are accurate (using finite differences to calculate the initial data can introduce small errors). First type

@r spectrum dat

which gives the output

  BR    PT  TY LAB    PAR(4)        L2-NORM       MAX U(1)      MAX U(2)      PAR(5)        PAR(6)        PAR(1)        PAR(2)
   1    15  EP   2  1.108419E-01  1.046496E+00  2.088819E-01  2.797320E-02  4.711503E-04  3.951714E-03  4.999823E-01  1.000000E+00

This AUTO run continues in γ (PAR(4) from the initial data (when γ=0) for a few iterations. We next continue in γ in the opposite direction and get AUTO to label the solution when γ=0. This gives a more accurate solution starting point than the one in the initial data. First of all save the previous solution by typing

@sv spectrum

then type

@R spectrum 2

which gives the output

  BR    PT  TY LAB    PAR(4)        L2-NORM       MAX U(1)      MAX U(2)      PAR(5)        PAR(6)        PAR(1)        PAR(2)
   1    14  UZ   3 -3.491091E-06  1.046496E+00  2.088819E-01  2.797321E-02  4.824896E-13 -1.262518E-07  4.999823E-01  1.000000E+00
   1    20  EP   4 -1.203557E-01  1.046496E+00  2.088819E-01  2.797319E-02  5.524371E-04 -4.280686E-03  4.999823E-01  1.000000E+00

Iteration 14 was determined by AUTO to have a γ (PAR(4)) value sufficiently close to zero and we use that solution as our starting point for the final continuation step. We first append the previous solution to those already calculated by typing

@ap spectrum

We then continue the solution in γ from γ=0 to γ=2π by typing

@R spectrum 3

which gives the output

  BR    PT  TY LAB    PAR(4)        L2-NORM       MAX U(1)      MAX U(2)      PAR(5)        PAR(6)        PAR(1)        PAR(2)
   1   226  UZ   5  6.283186E+00  1.046496E+00  2.088814E-01  2.797306E-02  8.254225E-02  1.390798E-01  4.999823E-01  1.000000E+00
   1   234  EP   6  6.518038E+00  1.046496E+00  2.088819E-01  2.797330E-02  8.185430E-02  1.443033E-01  4.999823E-01  1.000000E+00

We then save this output by typing

@sv spectrum.fig1

which results in the file p.spectrum.fig1. The part of the spectrum that has been calculated can be visualised by plotting the real (PAR(5)) against the imaginary (PAR(6)) parts of the eigenvalues from this file. Note that in this example, when c=0.5, the real part of the stability spectrum extends into the right hand complex plane (positive real part), meaning that this periodic travelling wave solution is unstable. If the above procedure is repeated for c=5.0 then the spectrum again extends into the right hand complex plane, and so the periodic solution is also unstable. However, for c=1.0, the part of the spectrum calculated by this procedure lies in the left hand complex plane. This is also true starting from any of the eigenvalues when γ=0. The periodic travelling wave is therefore stable in this case. In fact, we have found that in all cases for these equations, stability changes occur via an Eckhaus instability (see the next example), so that the portion of the spectrum near λ=0 is all that is required to detect stability.

NEXT EXAMPLE,

PREVIOUS EXAMPLE,

INDEX