J. Phys. Chem. 1993,97, 1566-1 570
1566
Asymptotic Pseudospectral Method for Reaction-Diffusion Systems E. Badlot and J. Boissonrde' Centre de Recherche Paul Pascal, Universits Bordeaux I, Avenue Schweitzer, F- 336OO Pessac, France Received: July 20, 1992
We show that the asymptotic method of Young and Boris is appropriate to the time integration of the stiff differential equations generated in the Fourier space when collocation spectral methods are applied to large reaction-diffusion model systems. This results in an easily implementable algorithm suitable for studies of morphology and stability of nonquilibrium structures in extended systems. increase the complexity of the sparsity pattern of the linear system, so that these methods are still more involved. It seems more Numerical simulations have always been priviledged methods appropriate to resort to spectral methods that use a decomposition in the study of chemical dissipativestructures.' The spontaneous on Fourier modes and automatically include these periodic formation of stationary spatial patterns resulting from the sole boundary conditions.21J2 There are other benefits in these competition of reaction and diffusion in systems kept far from methods. The errors due to space discretization converge faster equilibrium by a permanent supply of fresh chemicals-the sowith the number of modes than any finite difference scheme for called Turing structures-has received a special att~ntion.~-~ any order do with the number of points. One gets advantage of Interest in these problems has been recently stimulated by this property as soon as a sufficient number of modes is used to experimental evidences of Turing patterns in the real ~ o r l d .A~ . ~ reach the "spectral convergence" regime, that is, when the growing number of numerical simulations in two or three wavelength of the highest wavenumber mode is smaller than any dimensions have been performed in order to find out the typical length of the solution. Obviously, this is not the case morphology and the selection of the different possible patterns, when rapid spatial variations-such as steep fronts-are present either from a biological point of view or as a support for more (this problem can be sometimes circumvented by the choice of general theories of pattern formation in nonlinear system^.^-^^ a different base of functions but the method becomes more From a numerical point of view, the problem comes down to the complex). On the opposite, when the spatial oscillations are solution of a system of parabolic partial differential equations roughly harmonic-an ubiquitous situation nearby transitions points and in critical regimes-the spatial approximationis much &/dt = F(c) DA$ (1) more accurate than with finite differences, even with a small number of modes. Furthermore, with finite differencesa spatial where c is a vector of concentrations, F(c) is a nonlinear function resolution which seemssufficientto properly describe the solution, (reaction terms), D is a diagonal matrix of diffusion coefficients, does not actually guarantee a correct description of thedynamical and Arcis the Laplacian operator (diffusion terms). In the models processes. The tendency of these methods to develop large phase used in pattern formation theory, stiffness do not generally occur errors is a serious drawback in computations on wavy or from the reactive terms, in contrast with problems involving steep fronts or stiff chemical kinetics like combustion problems. unstationary processes. A good example, taken from a wave Nevertheless, numerical stability remains limited due to the space equation, is given in the introductionof ref 22 and another example discretization. In one dimension, high-performance finite difwill be illustrated in section 3 in the context of reaction4iffusion ference algorithms can be easily implemented.'J9.*0 When the systems. This weakness is quite logical since phases convey a dimensionalityincreases,the job becomes much more demanding, long range and global information, whereas finite difference especially if one is interested in the slow critical dynamics found methods are local by nature. In spectral methods, the phase in the vicinity of transitions or in extended systems dynamics information is accurately and globally characterized by the real where, due to the presence of topological defects, slow motions and imaginary part of each mode. As a matter of fact, most persist in the asymptotic regime. Although the rates given by eq nonlinear theories are expressed in terms of Fourier modes, so 1 are small, the time step is severelylimited by numerical stability that the spectral methods provide a more natural description and considerations which leads to large computation times. The an immediate perception of the internal dynamics. Thus, although instabilities are primarily related to the intrinsic stiffness caused finite difference methods remain most appropriate to numerous by the necessity of keeping spatial step size small enough to stiff chemistry problem in the real life, spectral methods seems accurately reproduce the solutions. To circumvent this problem, more suited to theoreticians working of the nonlinear dynamics one must resort to implicit or semiimplicitalgorithms;l9.*0in these properties of reaction-diffusion equations such as bifurcation methods one must generally solve, at each time step, one or several and pattern formation, pattern selection, spatiotemporal chaos, huge sparse linear systems involving the jacobian of the second or chemical turbulence. Among these spectral methods, collomember of eq 1. Block alimination supplemented by technical cation methods (often called "pseudospectral") should be preferred tricks (two-dimensional problems) or iterative techniques (twoto Galerkin methods for they make easier the treatment of the and three-dimensional problems) is used. The stability of the nonlinear terms. Unfortunately, stability problems, similar to whole process is limited by the nonlinear terms and must be those met in finite difference problems, also limit the time step controlled by appropriate time-step monitoring. Hunding has during integration. Raort to implicit techniques to bypass the developed such elaborated algorithms for the computation of problem is still more difficult and elaborate than with finite Turing structures in various spheroidal geometries on a vector differencestechniques since the nonlinear terms and the Laplacian computer,%I I 13. I 5 operator are evaluated respectively in the real space and in the Fourier space. For these reasons, spectral methodsdid not become In theoretical studies of extended systems, periodic boundary as popular for the solution of reaction-diffusion, as they did in conditions are often chosen in order to minimize the creation and the field of fluid dynamics. the freezingof topological defectsat boundaries. These conditions
1. Introduction
+
0022-3654/93/2097- 1566$04.00/0
Q 1993 American Chemical Society
The Journal of Physical Chemistry, Vol. 97, No. 8, 1993 1567
Reaction-Diffusion Systems At the same time, the growing availability of powerful desk workstations makes more insistent the necessity of algorithms that present a reasonable level of efficiency (Le., much higher than classical explicit methods) but remain very flexible and easy to implement with basic library programs like fast Fourier transforms. Fifteen years ago, Young and Boris proposed an asymptotic method to solve some stiff ordinary differential equations.23 The numerical scheme was applied to the direct solution of reaction4iffusion equations and implemented in a package.24 In spite of its great simplicity,this method has received much less attention than those based on more sophisticated stiff ODE'S integrator^.^^ In this paper, we show that, in the context of a pseudospectral (collocation) method, this technique is most appropriate to integrate the stiff ODE'S that rule the evolution of the Fourier components. The validity of the algorithm is demonstrated on a sensitive instability problem and comparisons are made with a finite difference method. 2. Tbe Metbod
For simplicity and in view of practical implementation, we shall illustrate our presentation on a particular two species model in one dimension. Extension to other models and higher dimensionalities is straightforward. The Schnackenberg is one of the chemical scheme used in the study of Turing structures.4J7 The rateequations can be reduced, after appropriate scaling, to the following form: aulat = r ( a - u
+
+ D, a 2 U / a X 2
(2)
U~U)
aulat = r ( b - U*V)
+ D, a 2 q a X z
(3)
where u and u are the concentration of the two variable species, a and b are concentrations kept constant and uniform, y is a rate constant, and D, and D, are diffusion coefficients. The length of the system was scaled to L = 2 r with periodic boundary conditions. As a result of a Turing type symmetry-breaking instability, this model is known to developstationaryconcentration patterns for proper sets of parameters which have been extensively studied el~ewhere.~J~ To set the ideas, we fix the parameters at pertinent values: y = lo4, D, = 4r2,D, = 8 0 ~ The ~ . values of a, b, u, and u are of order unity. In a collocation scheme, the Fourier transforms of u and u are replaced by the discrete Fourier transforms: 1 N-l 1 N-l 'k -~up-2riJkfN 6, = - ~ v , e - 2 r ~ j k f N NI=O Nil0
-
N
N
..,,--I (4) 2 2 where the uj and uj are taken at the N collocation points given by xj = 2 r j / N . The equations for the Fourier components t(k and z)k ("spectral space") form a system of 2N ordinary differential equations: k=--,
+
rcl
dfik/dt = Y[abk,o (uzu)k] - (7
+ D,k2)ck
(5)
The variables u and u at the collocation points can be obtained from their Fourier components Gk and 8k by the inverse discrete Fourier transform:
j = 0, ...,N - 1 (7) It seems more advisable to perform the integration in the
spectral space (eqs 5 and 6)than in the real space ( q s 2 and 3), since the diffusion operator is in diagonal form and since the periodic boundary conditions are automatically included in the Fourier transform. Any explicit integrator for ordinary differential equations can be used. Each evaluation of the rcl functions (u2u)k is performed by successive application of the inverse transforms (7), evaluation of the nonlinear products N uj2uj, and application of the direct transform (u20)k = (I/N) The efficiency of the process is achieved by useof fast Fourier transform algorithms (FFT)for real functions, available in all basic mathematical libraries.27~28 To ensure numerical stability, the time step h must be kept smallerthan a critical value Atc. For common explicit algorithms, the order of magnitude of Atc is given by the typical relaxation time of the system. For a zero wave vector, this relaxation time of the model is T 1/y. Unfortunately, for Fourier modes with large wave vectors, the relaxation time becomes small. For instance, with 64 collocation points, our relaxation time for lkl = 32 is T = 1/D,k2 lod, Le., 2 orders of magnitude smaller than for k = 0. Thus, although the global dynamics are rather slow, the time step must be kept smaller than Atc lod, Increasing the number of modes in order to improve the solution increases stiffness and still degrades the stability conditions. Implicit methods should be used to circumvent this type of limitations, but they generally call for some approximation of the jacobian matrix of the second member of equations (5,6). Since N the values of the (uzu)k,stemming from the nonlinear terms, are obtained by numerical Fourier transforms, they are not given in an explicit form. Thus the required partial derivativescannot be computed easily and these methods are awkward to implement in practice. The method of Young and Boris allows for an elegant solution that recovers the stability condition Atc l/y. We first summarize the principles of this method, originally developed to solve stiff differential equations in the real space. More details can be found in refs 20,23, and 24. We assume that each evolution equation can be written in the form
-
-
-
-
dy/dt = Q ( t ) - y / T ( f ) F(t) (8) where Q(r) and ~ ( tcan ) depend on all the variables of the system. The indices of the variables have been dropped for clarity. The function T defines a relaxation time. Let us call t s the typical evolution time of the system. All the equations are integrated simultaneously with a second order in time predictor-corrector method, iterations being repeated until convergence is achieved, but two different sets of formulas are applied to the different equations according to relaxation time value. If T t s eq 8 is said to be normal and can be integrated with the classical second-order corrector formula:
-
+ h ) -YO) = (h/2)[Fn(t+ h) + F0)l
vn(t
(9)
where F,(t + h) and yn(t + h) are respectively the values of F(t + h) and y(t + h) obtained at iteration n. A predictor formula is obtained at first iteration by setting Fo(t + h) = F(t): Yl(t
+ h ) = r(t) + hF(t)
(10)
If T
