Details on the computation of absolute spectra for periodic travelling waves in the Lambda-Omega equations using MATLAB and AUTO.

Part of the Supplementary Online Material for Smith, M.J., Rademacher, J.D.M. and Sherratt, J.A. Absolute stability of wavetrains can explain spatiotemporal dynamics in reaction-diffusion systems of lambda-omega type, in SIAM Journal of Applied Dynamical Systems, 8(3), 1136-1159, 27 August 2009

This online example illustrates how to calculate the absolute spectrum for periodic travelling wave solutions to the Lambda-Omega equations. The approach given is mostly based on the techniques outlined in Rademacher et al. (2007), and was partly implemented by modifying the AUTO code developed by Jens Rademacher for that paper. We use two different software packages in this analysis, MATLAB and AUTO. The MATLAB script can be run on any operating system with MATLAB and the "symbolic maths toolbox" installed. The instructions relating to AUTO will only work on UNIX/LINUX operating systems with AUTO installed, although modification of the methodology for other operating systems should be straightforward.

AUTO can be obtained from http://indy.cs.concordia.ca/auto/. We used AUTO 97 to run these examples, which required us to modify the AUTO problem size restrictions. It appears that newer versions of AUTO allow for larger problem sizes by default however we have not tested this. It is straightforward to recompile AUTO for larger problem sizes should you need to. Instructions to do this are in the AUTO manual (available at http://indy.cs.concordia.ca/auto/ ) under the heading "Restrictions on problem size". For your information we detail our own AUTO problem size limits in the file auto.h.

This online example is also an extension to the previous example we developed to illustrate the calculation of travelling wave families and essential stability using AUTO. We refer you to that page for more details.

Before carrying out these methods we recommend reading section 4 of the manuscript, which gives an overview of the procedure.

1. Obtain spatial and temporal eigenvalues at the branch points

In this first step we use MATLAB to calculate the values of λ and ν associated with the four branch points in the generalised absolute spectrum for a given travelling wave solution. In our study we found that these points alone can be used to predict the emergent spatiotemporal dynamics in simulations of the lambda-omega equations. Furthermore, these solutions are convenient intial starting values for AUTO continuation to calculate the generalised absolute spectrum. The MATLAB script can be viewed and downloaded by clicking Calcbranchpoints.m.

In this section we step through the mathematical rationale and procedure behind the MATLAB script. As a reminder, the lambda-omega equations are given by

(A1a) Eq1a

(A1b)Eq1b

where the variables u and v are the abundances of the interacting dynamical components, t is time and x is the spatial coordinate. For the purpose of studying the absolute spectrum we work with equations (A1) in amplitude-phase form as given by

(A2a) Eq2a

(A2b)Eq2b

where amplitude and phase . Kopell and Howard (1973) calculated an exact formula for the essential stability of periodic travelling wave solutions to equations (A2) that is given by

(A3) Eq3

Moreover, Sherratt(2003) showed that when equations (A1) are simulated on a sufficiently large semi-infinite domain with zero-Dirichlet boundary conditions, with a value of ω_1 such that condition (A3) is satisfied, then periodic travelling waves result with amplitude given by

 (A4)Eq4

Our study performs stability spectral analysis on waves selected by Dirichlet boundary conditions (given by equation (A4)) to predict the spatiotemporal dynamics emerging in PDE simulations.

The main manuscript outlines the linear stability analysis of solutions to equations (A2) that leads to the dispersion relation

(A5)

For each value of λ, this equation has four solutions of ν, which it is convenient to denote as

(A6)Eq6

The index, k, of ν_k(λ), is called the Morse index. For absolute stability, we must consider values of ν, such that Re(ν_i) = Re(ν_{i+1}) for some i = 1, 2 or 3. The set of such λ's is known as the “generalised absolute spectrum”. The absolute spectrum is a subset of the generalised absolute spectrum and is defined as the set of λ for which Re(ν_2) = Re(ν_3).

Rademacher et al. (2007) showed that the curves of generalised absolute spectrum are connected at singularities, where two of the four values of ν are equal. Such points are known as “branch points”. At such points the derivative of equation (A5) with respect to ν equals zero, which implies

(A7)Eq7

Given a value of ω_1 we can therefore use the wave selection formula (A4) and then equations (A5) and (A7) to calculate the ν and λ values at the branch points. This turns out to be a convenient way to calculate the absolute stability of travelling wave solutions because we find that the absolute spectrum attains its largest real part at branch points. It is also a convenient method for obtaining initial values for points on the generalised absolute spectrum, for the purposes of continuation.

The supplied MATLAB code therefore calculates a PTW solution for a user-defined value of ω_1 using equation (A4) and then uses equations (A5) and (A7) to calculate the λ and ν values associated with the branch points of the generalised absolute spectrum. Specifically, the code substitutes equation (A7) into (A5) to give a fourth order polynomial in ν, zeros of which give one value of ν for each of the four branch points, corresponding to the repeated roots. It then calculates the four λ values associated with these repeated roots by substituting them individually back into (A7). The code next substitutes these λ values individually into (A5) and solves for ν again to give all four ν_k values associated with each branch point.

The MATLAB code then reports whether the selected wave is absolutely stable or not; if any of the branch points have ν values with equal real part and Morse index 2 AND a value of λ with a positive real part then the selected wave is absolutely unstable. Secondly, the code prepares datafiles that act as starting values for AUTO continuation of the generalised absolute spectrum.

2. AUTO continuation of the generalised absolute spectrum

This section gives the methodology and example code for the continuations used to construct the generalised absolute spectra given in Fig. 4 of the main manuscript. Our general methodology is to perform continuations of solutions to the dispersion relation equation (A5), with the condition that the real parts of two of the ν values must be equal. This gives the generalised absolute spectrum, of which the absolute spectrum is the subset defined by the ν values with equal real part having Morse index 2.

Our example here steps through the continuation of the generalised absolute spectrum from each of the branch points when ω_1=1.0. This parameter value results in an essentially and absolutely stable wave being selected by zero Dirichlet boundary conditions. The four initial parameter value files for these continuations can be obtained by running the MATLAB script with ω_1=1.0 or can be downloaded by clicking LO_BP1.txt, LO_BP2.txtLO_BP3.txt or LO_BP4.txt. Note that LO_BP1.txt and LO_BP2.txt have ν values with equal real part whose Morse index is 1 and 3 and LO_BP3.txt and LO_BP4.txt have ν values with equal real part whose Morse index is 2.

The order of the ν_k values in these data files is important and it should be checked, if using this MATLAB code, that the first two ν values and their eigenvectors (explained below) are equal. If this is not true then the continuation will not work. This modification can be made by copying and pasting within the text files if the ordering of the ν values is not correct. To use the AUTO code provided here the data files need to follow the precise format of the datafiles produced by the MATLAB code.

Line of miscellaneous text
ω_0			ω_1			R
Re(λ)			Im(λ)
Re(ν_double1)	Im(ν_double1)
Re(_double1)	Re(_double1)	Im(_double1)	Im(_double1)		
Re(ν_double2)	Im(ν_double2)
Re(_double2)	Re(_double2)	Im(_double2)	Im(_double2)	
Re(ν_other1)	Im(ν_other1)
Re(_other1)	Re(_other1)	Im(_other1)	Im(_other1)	
Re(ν_other2)	Im(ν_other2)
Re(_other2)	Re(_other2)	Im(_other2)	Im(_other2) 

where and are eigenvectors associated with specific values of ν.

For each continuation we use the dispersion relation in matrix form, as given by equation (10) in the main manuscript, rewritten here in extended form as (A8a) (A8b)

where k=1,2,3,4, and and are the eigenvectors associated with specific values of ν_k. We also impose the necessary constraint that two of the values of ν_k have the same real part, which ensures that continuations proceed along the generalised absolute spectrum.

It is convenient to perform a simultaneous continuation for all four values of ν_k to allow their Morse index to be monitored. For our system this necessitates the simultaneous contination of sixteen equations in total (4 ν_k values in 2 equations for separate real and imaginary parts). In addition we need initial values for all of the eigenvectors and . These eigenvectors are also calculated by the MATLAB code. The procedure is to substitute ν_k, λ, ω_1 and R into equation (A8) to give ; we then take and .

A further complication, as noted in the manuscript, is that a modification to equations (A8) is necessary to allow continuations to begin from the branch points because two of the ν_k values are equal and AUTO cannot distinguish them at the initial continuation step. The necessary modification is a desingularisation procedure outlined in Rademacher et al. (2007). This allows the first two equal ν_k values in the datafiles to be separated in the initial continuation step (their imaginary parts become different but their real parts remain equal). For this desingularisation step we replace the four equations for one of the repeated ν_k values with the desingularised equations (A9a) (A9b)

where and are eigenvectors that are estimated during the initial AUTO continuation step, and initially η=0.

The AUTO equations file used in this example is LO.f. This code reads in the initial data from the file LO.dat when run for the first time. It has a parameter "it" to denote whether to perform the continuation using the desingularised equations ("it=0") or the normal equations ("it=1").

For our problem we have 16 equations and require 16 boundary conditions and 8 integral constraints. A detailed discussion on why we need these conditions, and the form they take, is outwith the scope of this tutorial but see Rademacher et al. (2007) for details.

Our initial AUTO run from the branch point performs a continuation in a dummy parameter (PAR(14)) to allow AUTO to obtain an appropriate eigenvector for the desingularised equations. The LO.f file obtains its initial values from LO.dat. For this first example rename the file LO_BP1.txt to LO.dat. The parameters file for this first continuation is the file r.LO-1a. Running LO.f with this parameters file gives the output (click to view pdf). This shows that constant values were obtained for all of the parameters (PAR(14) is the dummy parameter).

The next step is to perform a continuation from the branch point by increasing or decreasing η from η=0. The parameters file to do this for this branch point is r.LO-1b. This continuation gives the output.

Which shows that 50 continuation steps were performed, increasing η (PAR(9)) from η=0 to η=1.2... and obtaining corresponding values for the other parameters. In this case we can see that at this final continuation step the real parts of ν_3(PAR(10)) and ν_4(PAR(12)) are equal and are less than the real parts of ν_1=ν_2(PAR(6)). This part of the generalised absolute spectrum therefore has two pairs of ν values with equal real parts, whose Morse indices are 1 and 3.

The final continuation step is to continue using the normal equations (not desingularised) in Re(λ)(PAR(4)) for a greater number of iterations. This is really an optional step and continuations could also be performed in η if desired. The parameters file to do this is r.LO-1c and you can also modify the parameter "it" in LO.f to be "it=1". Example output from such a continuation is given by the output.

In this case the continuation ended because Re(λ)(PAR(4))=30, a termination point defined in the parameters file. Analysis of the full data from this continuation shows that λ remains on the real line and with ν values with equal real part having Morse index 1 and 3. An identical continuation procedure can be performed for the other branch point with λ real. And the procedure for the two branch points with complex conjugate λs is also similar. However, these latter cases raise a further complication to do with the changing of the Morse index at "triple points". The next section describes how to deal with this.

3. Dealing with triple points in AUTO continuation.

There are several possible ways to deal with the occurrence of triple points in the generalised absolute spectrum of which we detail the one we used in this study. We begin by starting with one of the complex conjugate branch points, with the initial conditions file LO_BP3.txt (copied to LO.dat). We first use the parameters file r.LO-3a with LO.f to obtain starting values. The next continuation used the parameters file r.LO-3b to continue in η. This gives the output

Note that an additional user-defined output ocurred at label 7. This occurred because the continuation went through a triple point. In this case Re(ν_1)=Re(ν_2)=Re(ν_4) (PAR(6)=PAR(6)=PAR(12)). (Note that generally a triple point only requires that the real parts of three ν values are the same). If we continue the procedure in the normal way, using the parameters file r.LO-3c, then we locate a second triple point in the output

Examination of the datafiles for these continuations reveals that the ν values with equal real part have Morse index 2 at the branch point, this changes to Morse index = 3 after the first branch point is crossed and then changes back to Morse index = 2 beyond the second triple point (label 9).

We wish to identify the missing part of the absolute spectrum that occurs between these triple points. To do this we save the previous continuation as a separate file, e.g. using the command ("@sv LO-TP"). We then copy one of the saved files to another file using the command "cp q.LO-TP LO2.dat". We then use LO2.dat as an initial parameter file for a second AUTO equations file LO2.f. This AUTO file is almost identical to LO.f but has been modified to read LO2.dat and swap the first and the fourth values of ν. In this case, this procedure is sufficient to enable AUTO to continue along a different branch of the generalised absolute spectrum; implementation with the parameters file r.LO2-3a gives the output

Examination of the data file reveals that the continuation proceeds in Re(λ) along the real line with the morse index of the ν values with equal real part flipping between 2 and 3 at the branch points. Moreover, the continuation automatically turns at the two branch points that are real in λ with ν values with equal real part that have Morse indices 1 and 3.