## Abstract

Many recent datasets on cyclic populations reveal spatial patterns with the form of periodic travelling waves (wavetrains). Mathematical modelling has identified a number of potential causes of this spatial organization, one of which is a hostile habitat boundary. In this paper, the author investigates the member of the periodic travelling wave family selected by such a boundary in models of reaction–diffusion type. Using a predator–prey model as a case study, the author presents numerical evidence that the wave generated by a hostile (zero-Dirichlet) boundary condition is the same as that generated by fixing the population densities at their coexistence steady-state levels. The author then presents analysis showing that the two waves are the same, in general, for oscillatory reaction–diffusion models with scalar diffusion close to Hopf bifurcation. This calculation yields a general formula for the amplitude, speed and wavelength of these waves. By combining this formula with established results on periodic travelling wave stability, the author presents a division of parameter space into regions in which a hostile boundary will generate periodic travelling waves, spatio-temporal disorder or a mixture of the two.

## 1. Introduction

Cycles of abundance in natural populations were one of the first predictions made by mathematical biology [1,2]. Over the past two decades, field data have shown that, in many cases, such population cycles are not synchronous across a habitat, comprising instead spatio-temporal patterns. The first reported instance of this was periodic travelling waves (ptws) in vole population numbers [3,4], meaning that there is a regular spatial pattern of abundance that moves across the domain at a constant speed. Other examples include ptws and patterns of regional synchrony in geometrid moths in Northern Fennoscandia [5–9], ptws and intermittent synchrony in red grouse [10,11], and ptws in larch budmoth in the European Alps [12–14]; see Sherratt & Smith [15] for further examples.

Modelling studies have demonstrated a number of potential causes for this spatial asynchrony, including invasions [16–21], heterogeneous habitats [13,22–24], migration between subpopulations [25], migration driven by pursuit and evasion [26] and hostile habitat boundaries [27–30]. This paper concerns the last of these. Hostile boundaries arise naturally at abrupt changes in landscape: for example, one of the best datasets on (cyclic) red grouse is from Kerloch Moor in northeastern Scotland, UK [11,31], which is bounded to the north by farms and woodland that is unsuitable habitat for red grouse, and which they do not cross [32]. I have shown previously that, for cyclic populations, a hostile boundary generates a ptw [27,28]. In this paper, I will investigate analytically the form of this ptw in reaction–diffusion models for which the local population dynamics have a low-amplitude limit cycle. I begin in §2 by discussing the generation of ptws by hostile boundaries in predator–prey systems; this acts as an illustrative and motivating example. In §3, I focus on low-amplitude cycles, and derive the leading-order form of the ptw generated by a hostile boundary when the local population dynamics are close to a Hopf bifurcation. In §4, I use this result to determine the stability of the ptw, and I consider the form of the solution when the ptw is unstable. I conclude in §5 by discussing the basis for hostile boundaries in behavioural ecology.

## 2. Boundary-generated periodic travelling waves in predator–prey systems

To provide a specific illustration of ptw generation by hostile boundaries, I consider the Rosenzweig–MacArthur model [33] for predator–prey interaction,
2.1aand
2.1bThis is a dimensionless formulation of the model in which *p* and *h* denote, respectively, predator and prey densities at space point *x* and time *t*; throughout the paper, I will consider only one-dimensional spatial domains. I will also assume throughout that the two populations have equal dispersal rates. This is appropriate for most aquatic systems, but many terrestrial predators disperse more rapidly than their prey. The prospects for extending the results in this paper to unequal dispersal rates will be discussed in §5.

In (2.1), the prey consumption rate per predator is an increasing saturating function of prey density with Holling type II form: *C*>0 reflects how quickly the function saturates. Parameters *A*>0 and *B*>0 are dimensionless combinations of the birth and death rates. Provided *A*>1+1/*C*, equations (2.1) have a unique homogeneous coexistence steady state (*h*_{s},*p*_{s}), where *h*_{s}=1/(*AC*−*C*) and *p*_{s}=(1−*h*_{s})(1+*Ch*_{s})/*C*. This is stable as a solution of the kinetic odes only if *C*<*C*_{Hopf}=(*A*+1)/(*A*−1), with population cycles for *C*>*C*_{Hopf}.

At a boundary that is hostile for the prey (say), the correct condition is in most cases of Robin type, ∂*h*/∂*x*=*Qh* [34,35]. Here, the parameter *Q* is large, reflecting the high mortality rate that would be experienced by individuals outside the habitat. It follows that the Robin condition can be approximated by *h*=0, and this approximation is in widespread use. In fact, I have shown previously [36] that, for ptw generation, the Robin condition is particularly well approximated by the Dirichlet condition, which I will use throughout this paper.

Figure 1*a* illustrates a typical solution of (2.1) for *C*>*C*_{Hopf}, on a large domain with *h*=*p*=0 at the left-hand (hostile) boundary and the zero-flux condition ∂*p*/∂*x*=∂*h*/∂*x*=0 at the other boundary. After transients have decayed, the solution settles to a form that is independent of initial conditions (see [27]). There is a ptw in the bulk of the domain, moving in the positive *x* direction, with thin transition layers near the two boundaries. Detailed numerical investigation [27,28] shows that this behaviour is driven by the left-hand (hostile) boundary condition, with the zero-flux condition at the right-hand boundary playing no significant role. Intuitively, one can regard the hostile boundary condition as the factor that forces the solution away from spatially uniform oscillations.

For any given values of *A*, *B* and *C*>*C*_{Hopf}, (2.1) has a one-parameter family of ptw solutions; Kopell & Howard [37] showed that this is a general property of oscillatory reaction–diffusion systems with scalar diffusion, meaning that the diffusion coefficient is the same in each equation. Figure 2 illustrates this family for the parameter values used in figure 1, plotting the wavelength as a function of the prey amplitude. One particular wave in this family (indicated by the large dot) is selected by the condition at the hostile boundary. There is currently no analytical solution to this wave-selection problem for almost all oscillatory systems, including (2.1). However, Sherratt [41] obtained the solution for the system
2.2aand
2.2bwhere , with subscripts *x* and *t* denoting partial derivatives. This belongs to the ‘*λ*–*ω*’ class of equations [37,42]. Sherratt [41] considered (2.2) on a semi-infinite domain with *u*=*v*=0 at the boundary. I showed that, as transients decay, the solution away from the boundary approaches the ptw
2.3a
2.3b
and
2.3c

The system (2.2) is not a model for any particular physical or biological phenomenon. Rather, its significance is as the normal form of a reaction–diffusion system with scalar diffusion and oscillatory kinetics close to a standard supercritical Hopf bifurcation. For example, the reduction of the Rosenzweig–MacArthur model (2.1) to normal form is described in the appendix of Sherratt [18]. This gives equations of the form (2.2) with *ω*_{0} and *ω*_{1} depending on *A*, *B* and *C* as follows:
2.4aand
2.4bEquations (2.2) and (2.4) determine the leading-order behaviour of (2.1) as .

Under the reduction to normal form, it is the coexistence steady state (*p*,*h*)=(*p*_{s},*h*_{s}) that maps to (*u*,*v*)=(0,0). Therefore, the boundary condition *u*=*v*=0 in the *λ*–*ω* system corresponds to *p*=*p*_{s}, *h*=*h*_{s} in the predator–prey system. This boundary condition is relevant to some applications, for example when there is an abrupt change in habitat causing a transition from cyclic to non-cyclic dynamics [43], but such transitions are very much rarer than hostile boundaries. Figure 1*b* illustrates the generation of ptws by this boundary condition for (2.1). The solutions in figure 1*a* (*p*=*h*=0 at *x*=0) and figure 1*b* (*p*=*p*_{s}, *h*=*h*_{s} at *x*=0) are for the same values of the parameters *A*, *B* and *C*. Although there are many qualitative similarities between the two solutions, a fundamental difference is that the ptws are not the same in the two cases (figure 2).

I performed a detailed numerical study of the members of the ptw family selected by the boundary conditions (*p*,*h*)=(0,0) and (*p*,*h*)=(*p*_{s},*h*_{s}). I calculated the wavelength and prey amplitude of these ptws for a range of different *C* values, with *A* and *B* fixed. A quantitative study such as this demands careful investigation of numerical accuracy, and I present full details of my numerical method and convergence tests in appendix A. Figure 3 illustrates my calculations, with results for (*p*,*h*)=(0,0) and (*p*,*h*)=(*p*_{s},*h*_{s}) plotted in grey and black, respectively. Although the ptws generated by the two boundary conditions are in general different, they appear to become the same as *C* approaches *C*_{Hopf}. This suggests that, although the formula (2.3) was derived for a boundary condition equivalent to (*p*,*h*)=(*p*_{s},*h*_{s}), it might also be relevant for the much more widespread condition (*p*,*h*)=(0,0). In §3, I will present an analytical study showing that (2.3) does indeed give the ptw selected by the hostile boundary condition (*p*,*h*)=(0,0), to leading order as .

## 3. Periodic travelling wave generation for *C* near *C*_{Hopf}

### (a) Analytical investigation

The Hopf bifurcation in the kinetics of (2.1) is of standard supercritical type [44], example 3.1, so that (2.2) is the appropriate normal form. More precisely, there are invertible coordinate and parameter changes and a time reparametrization that together transform (2.1) into
3.1aand
3.1bwith correction terms that are *O*(*R*^{4}) [[44], §3.5; [45], pp. 150–152]. Here, *R*=(*U*^{2}+*V* ^{2})^{1/2} and *μ*=*C*−*C*_{Hopf}. The constants *Λ*_{0}, *Λ*_{1}, *Ω*_{0}, *Ω*_{1} and *Ω*_{2} depend on the parameters *A* and *B*.

From the viewpoint of the normal form, the coordinate change involved in transforming (2.1) to (3.1) is relevant only in the vicinity of the coexistence steady state (*p*_{s},*h*_{s}). However, it can be applied globally to (2.1), and I denote by (*U*_{bdy},*V* _{bdy}) the image of *p*=*h*=0 under this coordinate change. I now rescale (3.1) using
which gives (2.2) with
3.2It follows that, for small *μ*, the leading-order behaviour of (2.1) subject to the hostile boundary condition *p*=*h*=0 is given by the solution of (2.2), (3.2) subject to the boundary condition
3.3Here, I am assuming a semi-infinite spatial domain, which I will take as *x*>0. I will present formal arguments showing that, away from the boundary, this leading-order behaviour consists of the ptw (2.3). It is important to emphasize that my arguments are not rigorous, and a rigorous proof is a natural challenge for future work.

Following standard practice [37,42], I rewrite (2.2) in terms of the amplitude *r*=(*u*^{2}+*v*^{2})^{1/2} and phase . This gives
3.4aand
3.4bon *x*>0, with the hostile boundary condition being
3.5aand
3.5bat *x*=0. A typical numerical solution of (3.4) subject to (3.5) is shown in figure 4. Numerical solutions must of course be done on a finite domain, and I impose the zero-flux condition *r*_{x}=*θ*_{x}=0 at *x*=*L* for some suitably large *L*, but this only affects the solution in the immediate vicinity of that boundary. Otherwise, the solution at large times consists of a thin region near *x*=0 in which *r* and *θ*_{x} vary spatially and oscillate in time, with the majority of the domain having constant *r* and *θ*_{x}, corresponding to a ptw in *u* and *v*.

The spatial variations in *r* and *θ*_{x} become increasingly localized near *x*=0 as *μ* (>0) is decreased, with the frequency of the temporal oscillations also increasing. Therefore, I consider separately the ‘outer’ solution in the bulk of the domain, and a boundary layer near *x*=0. Figure 5 presents detail of the solution near *x*=0 for a small value of *μ*. This shows that the temporal oscillations in *r* and *θ*_{x} are more localized to *x*=0 than the spatial variation. Therefore, I look for outer solutions in which *r* and *θ*_{x} are independent of time, with the form
3.6as . Here, *K* is a constant, independent of *μ*; the *μ*^{−1}*αt* term arises as a consequence of the corresponding term in (3.4*b*). Substituting (3.6) into (3.4) gives
3.7aand
3.7bwhere .

For the boundary layer, I substitute
into (3.4). Neglecting terms that are *O*(1) as , this gives
3.8aand
3.8b

I will not consider the full solution of the boundary layer equations (3.8). Rather I focus entirely on the behaviour for large . This behaviour must match the outer solution, which is independent of time. It follows that the oscillations in the boundary layer must decay at a rate that is beyond all algebraic orders (e.g. exponential) at large . Therefore, for matching between the boundary layer and outer solutions, it is sufficient to consider time-independent solutions of (3.8). Moreover, any non-zero value of at large could not be matched by the outer solution. It follows that as
3.9with as . Here, ‘t.s.t.’ denotes ‘transcendentally small terms’. The constant of integration *K* must be the same as in the outer solution (3.6) to enable matching. Substituting (3.9) into (3.8) and equating leading-order terms gives
3.10aand
3.10bwhere .

To remove the singularity in (3.10) at I make the change of variables
3.11Note that *ξ* is an appropriate coordinate because is a (rescaled) solution amplitude, and must therefore be non-negative. Substituting (3.11) into (3.10) gives
3.12a
3.12b
and
3.12c
The steady states of (3.12) are where
3.13I require a solution of (3.12) with as . However, this does not necessarily imply that *p* and *q* tend towards steady-state levels: a closed *p*–*q* loop in the plane would also be a suitable asymptotic solution. I have not proved that there cannot be such a loop, but numerical calculations of the phase plane for (3.12*a*,*b*) show that the end points of trajectories are either infinity or one of the two steady states (figure 6). Moreover, the steady state with *p*>0 is not a suitable asymptote because , which means that cannot be bounded away from zero and positive. Therefore, the required solution of (3.12) must approach as . Calculation of the eigenvalues at these steady states shows that there is a unique eigenvalue −*Γ* with negative real part.

This discussion implies that, for a solution of (3.12) to have a form suitable for matching to the outer solution, it must have for some constant *A*>0 as . Also . Using the definition of *ξ* in (3.11), these translate to
3.14as .

I turn now to the outer solution. The leading-order behaviour of this as must match (3.14). This requires
3.15I have not attempted a systematic investigation of the solution of (3.7). Rather, noting that as , I looked for a solution in which . Substituting this constraint into (3.7) yields the exact solution
3.16a
3.16b
and
3.16cBecause as , this solution satisfies (3.15). I make no claims of uniqueness for (3.16), and there may be solutions satisfying (3.15) for either the same or different values of *K*. However, numerical results (discussed in §3*b*) are consistent with (3.16) being the appropriate solution.

The key implication of (3.16) is that as the solution approaches a ptw with amplitude and phase gradient opposite in sign to *ω*_{1}. The formula (3.13) for *Γ* shows that this is exactly the same ptw as that given in (2.3). Hence, to leading order for small *μ*, the ptw solutions of (2.1) generated by the boundary conditions *p*=*p*_{s}, *h*=*h*_{s} and *p*=*h*=0 are the same. In fact, my calculations do not depend on the values of *u*_{bdy} and *v*_{bdy}, and thus any Dirichlet boundary conditions select the same ptw to leading order for small *μ*.

### (b) Numerical verification

The non-rigorous nature of my calculations means that it is important to test them numerically. Conceptually, this is straightforward: one solves (2.2), (3.2) on a relatively large domain with (3.3) at one boundary and a zero-flux condition at the other. After a solution time that is sufficiently long to allow transients to decay, one can then read off the amplitude at the centre of the domain. I used a semi-implicit finite difference scheme, and, to improve accuracy, I solved using three different time steps and applied a convergence acceleration formula. Details of this are given in appendix A, which also discusses the various numerical convergence tests that I performed.

The basic difficulty with this procedure is that, as *μ* decreases, *ω*_{0} increases. In many situations, the value of *ω*_{0} is irrelevant to the solution of (2.2) because *ω*_{0} does not feature in the phase-amplitude equations (3.4), and in the physics literature most authors remove *ω*_{0} via a ‘gauge transformation’ [42]. However, (3.5) is not invariant under a gauge transformation, and the value of *ω*_{0} does affect the solution I am studying. Moreover, my numerical convergence tests show that preservation of accuracy requires the numerical time step to decrease in proportion to *μ*^{2}. Consequently, numerical solutions rapidly become computationally expensive as *μ* is decreased, and it is unfeasible to solve for very small values of *μ*.

Investigation of higher-order terms in the asymptotic expansions of *r* and *θ* for small *μ* shows that the leading-order correction to (2.3) is *O*(*μ*^{1/2}). Therefore, I expect an approximately linear relationship between ptw amplitude and *μ*^{1/2} when *μ* is small. I used this to circumvent the inability to determine the ptw amplitude numerically for very small *μ*. I determined the values of the amplitude for a number of relatively small *μ* values and then performed a linear regression analysis for amplitude as a function of *μ*^{1/2}. This enabled an effective estimation of the amplitude when *μ*=0. Figure 7 illustrates this procedure, showing that it confirms that (2.3*b*) gives the limiting amplitude as .

## 4. Periodic travelling wave stability

Figure 2 illustrates the family of ptws that exist as solutions of (2.1) for a given set of parameters. An important subsidiary issue is the stability of these solutions. I am not aware of any results on the stability of solutions, such as those illustrated in figures 1 and 4, that satisfy a Dirichlet condition at one boundary of a semi-infinite domain, and tend asymptotically to a ptw. However, stability of ptws on infinite domains has been well studied. In particular, for oscillatory reaction–diffusion equations with scalar diffusion, it is a general result that ptws of sufficiently low amplitude are unstable, whereas those of high amplitude are stable [37,48]). In figures 1 and 3, I have chosen parameter values for which the particular member of the ptw family selected by the hostile boundary condition is stable as a solution of (2.1). However, this is not always the case. Figure 8 shows two solutions of (2.1) with a hostile boundary for parameter values giving an unstable ptw: as a result, there are disordered spatio-temporal oscillations.

Unstable solutions of a spatio-temporal system are either ‘convectively’ or ‘absolutely’ unstable [49]. In the former case, all unstable linear modes propagate, whereas an absolutely unstable solution has stationary unstable linear modes. These two types of instability lead to qualitatively different solution forms in a domain with a hostile boundary. As well as generating a ptw, this boundary applies a potentially destabilizing perturbation to that ptw. If the ptw is convectively unstable, then the unstable components of this perturbation all travel away from the boundary as they grow. Therefore, the ptw remains intact close to the boundary, with a transition to spatio-temporal disorder occurring in the interior of the domain. An example of this is shown in figure 8*a*. By contrast, if the ptw is absolutely unstable, then growing linear modes remain fixed in the vicinity of the boundary, and spatio-temporal disorder occurs throughout the domain (figure 8*b*).

For a general oscillatory system such as (2.1), even numerical determination of whether a ptw is convectively or absolutely unstable is a major challenge [50,51]. However, for the *λ*–*ω* system (2.2), detailed calculations can be carried out using the phase-amplitude equations (3.4), for which ptws are homogeneous solutions. Using this approach, it has been shown that the ptw solution (2.3) of (2.2) is absolutely unstable if |*ω*_{1}|>1.576465 and convectively unstable if 1.110468<|*ω*_{1}|<1.576465 [52]; the condition for stability is |*ω*_{1}|<1.110468 [37,41]. The combination of these results with (2.4) makes it straightforward to predict the stability of the ptws generated by a hostile boundary for the model (2.1), when *C* is close to *C*_{Hopf}. Figure 9 shows the division of the *A*–*B* plane into the three cases of stable, convectively unstable and absolutely unstable ptws. The corresponding solutions have the same qualitative form as figures 1, 8*a* and 8*b*, respectively. This ability to make a detailed prediction of the qualitative solution form is an important consequence of the formula (2.3) for the ptw generated by a hostile boundary for parameters close to Hopf bifurcation.

## 5. Discussion

Historically, most field studies on cyclic populations monitored only temporal changes in abundance. However, the past two decades have seen an increasing number of spatio-temporal datasets [15], together with advances in spatio-temporal statistics that permit their analysis [53–55]. In many cases, this has revealed ptws. Hypotheses on the causes of this spatial organization are common, but I am not aware of any example in which the cause has been established empirically. Therefore, mathematical modelling has played a key role, establishing various factors as viable causes of ptw formation. One of these is a hostile habitat boundary. Simulation-based studies have shown that this is an effective generator of ptws in cyclic populations [27,29,30,41]. However, there remains a pronounced lack of analytical results on this mechanism. The key contribution of this paper is a formula for the amplitude (and hence wavelength and frequency) of the ptw generated by a hostile habitat boundary, valid to leading order as the parameters approach a Hopf bifurcation in the local dynamics. This formula enables, for the first time, quantitative predictions of the way in which ecological parameters affect the spatio-temporal oscillations generated by hostile boundaries. Throughout the paper, I have illustrated my results using the Rosenzweig–MacArthur model for predator–prey interactions. However, my formula is valid for any reaction–diffusion model of a cyclic population with scalar diffusion.

The significance of the restriction to scalar diffusion is that the appropriate normal form close to Hopf bifurcation is the *λ*–*ω* system (2.2). More generally, the normal form contains a linear dispersion (i.e. cross-diffusion) term, making it a complex Ginzburg–Landau equation [42]. Some relevant results are known for this case. In particular, an exact formula for the ptw generated by a zero-Dirichlet boundary condition is known: it is the asymptotic wave for a stationary Nozaki–Bekki hole [56,57]. This boundary condition corresponds to fixing the populations at their coexistence steady state in the cyclic population model, and it would provide a starting point for studying the case of non-scalar diffusion. Numerical investigation shows that new phenomena can occur in this case, including folds in the ptw solution branch and non-trivial stationary waves [58]. A key step in my calculation is the exact solution (3.16). Bekki & Nozaki [56] calculated their exact hole solutions using Hirota's bilinear method [59], and it is possible that this method could be used to derive the generalization of (3.16).

An example of a cyclic population for which ptws occur, and for which a hostile boundary is a plausible cause, is field voles in Kielder Forest (northern England, UK). This population cycles with a period of about 4 years, and extensive spatio-temporal fieldwork has shown that the cycles are organized into ptws, moving across the forest at about 15–20 km yr^{−1} [4,60]. Kielder Forest contains a large (10 km^{2}) reservoir—one of the largest man-made bodies of water in Europe. Intuitively, one expects that the large open space of the reservoir will facilitate the hunting of voles by birds. These are a major predator: across the UK as a whole, birds are responsible for about 30 per cent of vole predation [61], and data for Kielder Forest itself shows that tawny owls alone remove 10–15% of the vole population per annum [62]. To the best of my knowledge, there are no data on the spatial distribution of avian predation, but its expected localization near the reservoir edge would make this habitat boundary very hostile for voles, suggesting that this might be the cause of the observed ptw in vole abundance.

The response of natural populations to habitat boundaries is an active area of research within behavioural ecology. The best-studied population type is forest-dwelling birds. Analyses of movement trajectories in response to call playback and during homing over long distances both show a strong tendency to use long within-forest paths in order to avoid crossing gaps [63–65]. Butterflies have also been studied extensively, and show a strong preference for moving in habitat corridors rather than crossing open areas [66–68]. Similarly, some fish travel significant distances within seagrass or reef habitats to avoid crossing gaps [69,70]. For mammals, most of the relevant data concern the behavioural response to roads. Hedgehogs [71], voles [72], mice [72], chipmunks [73] and moose [74] all exhibit significant avoidance of large roads. There are also a few studies of the response of mammals to natural habitat boundaries. For example, the forest-dependent red squirrel [75] and the grassland-dwelling Franklin's ground squirrel [76] both show a preference for detours over crossing habitat gaps. Detailed study of these various behaviours would require a much more detailed model than linear diffusion for animal movement [77]. However, within the context of a reaction–diffusion model, a zero-Dirichlet condition is an appropriate representation of avoidance of a habitat boundary. Therefore, these various behavioural studies argue strongly for the importance of understanding in detail the implications of hostile boundary conditions for overall population dynamics.

## Appendix A. Numerical methods

I solved (2.1) and (2.2) using a semi-implicit finite difference (Crank–Nicolson) scheme with a uniform space step and a constant time step, which I denote by *δx* and *δt*, respectively.

**A.1. Numerical solution of (2.1)**

My simulations were all for *A*=3 and *B*=4, and all used a Dirichlet condition at the left-hand boundary *x*=0 (either *p*=*h*=0 or *p*=*p*_{s}, *h*=*h*_{s}) and zero flux (∂*p*/∂*x*=∂*h*/∂*x*=0) at the right-hand boundary *x*=*L*, say. I solved for 0<*t*<*T*; the values of *L* and *T* were chosen so that the effect of further increase on the wavelength and prey amplitude of the ptws was less than that of errors owing to numerical discretization. For my final choices of space and time steps (discussed below), *L*=800 and *T*=24 000 were appropriate provided *C*≥2.2, but larger values of *T* were required for *C* closer to *C*_{Hopf}=2.

At the end of the simulation, I calculated the wavelength as the average distance between points at which *h*=*h*_{s} with ∂*h*/∂*x*>0, within the region . For the amplitude of *h*, I calculated the mesh points with the largest and smallest values of *h*, again in the region . For both of these points, I fitted a quadratic through the point and its immediate neighbours; the turning points of these quadratics gave the maximum and minimum of *h*, with the amplitude being their difference.

Convergence tables of the wavelength and prey amplitude for the two Dirichlet boundary conditions are given in the electronic supplementary material. The error is *O*(*δt*+*δx*^{2}). Based on these convergence tables, I used *δx*=0.25 and *δt*=0.003 for figure 1, which give errors of a little under 0.1 per cent. For figure 3, greater accuracy was required, and I achieved this using convergence acceleration (see [78], §3.9 for review). Because convergence is linear in *δt*, I applied the Aitken acceleration formula using the values given by *δt*=0.012, 0.006 and 0.003. I calculated these Aitken values for two spatial discretizations, *δx*=0.5 and 0.25, and then obtained a final numerical estimate by applying the Richardson extrapolation formula. This is more appropriate than the Aitken formula because of the quadratic convergence in *δx*. The final values given by this procedure are accurate to within 0.003 per cent for my test case *A*=3, *B*=4, *C*=5.

**A.2. Numerical solution of (2.2)**

In comparison with the predator–prey model (2.1), estimation of numerical accuracy for (2.2), (3.2) subject to (3.3) is significantly easier. This is because for *u*_{bdy}=*v*_{bdy}=0 the ptw amplitude is known to be given by (2.3*b*). I used the calculated value of this amplitude to assess numerical convergence for the test case *α*=1.5, *β*=0 and *ω*_{1}=0.8; these are the parameter values used in figures 4, 5 and 7. A convergence table of the ptw amplitude for these parameters with *μ*=1 is given in the electronic supplementary material. This also tabulates numerical convergence with the time step *δt* (for fixed *δx*) for , and . This shows that, as *μ* is decreased, the time step must also be decreased to maintain accuracy, with *δt*∝*μ*^{2} giving an approximately constant accuracy level.

Based on these convergence tables, I used *δx*=0.04 and *δt*=0.0025 for figure 4, which has ; this gives an error of about 0.2 per cent in the ptw amplitude. For figure 7, greater accuracy is required, and the small values of *μ* used in some of the simulations lead to potentially very long run times. Therefore, I ran for *δt*=0.04*μ*^{2}, 0.02*μ*^{2} and 0.01*μ*^{2} (with fixed *δx*=0.2), after which I estimated the ptw amplitude using the Aitken acceleration formula [78], §3.9; this is appropriate because of the linear convergence in *δt*. The resulting error is about 0.02 per cent.

Finally, it is important to comment that the case *u*_{bdy}=*v*_{bdy}=0 used for my convergence tests does have the special feature that there are no spatio-temporal oscillations near the left-hand boundary. Such oscillations might, in principle, affect the numerical accuracy of the ptw amplitude, especially for small *μ* when they become more localized and of higher frequency. As a test of this I compared the percentage errors in the amplitude of the ptw generated by *u*_{bdy}=*v*_{bdy}=0 and *u*_{bdy}=*v*_{bdy}=1; in the latter case, I calculated errors relative to a high-accuracy estimate of the amplitude. This showed that the errors were approximately the same in the two cases (see the electronic supplementary material for details).

- Received December 31, 2012.
- Accepted March 13, 2013.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.