GDNL#3/97 Revised version arXiv:patt-sol/9712007v1 23 Dec 1997 Validation and Calibration of Models for Reaction-Diffusion Systems Rui Dilão1 and Joaquim Sainhas2 Grupo de Dinâmica Não-Linear, Departamento de Fı́sica Instituto Superior Técnico, Av. Rovisco Pais 1096 Lisboa Codex, Portugal Abstract Space and time scales are not independent in diffusion. In fact, numerical simulations show that different patterns are obtained when space and time steps (∆x and ∆t) are varied independently. On the other hand, anisotropy effects due to the symmetries of the discretization lattice prevent the quantitative calibration of models. We introduce a new class of explicit difference methods for numerical integration of diffusion and reaction-diffusion equations, where the dependence on space and time scales occurs naturally. Numerical solutions approach the exact solution of the continuous diffusion equation for finite ∆x and ∆t, if the parameter γN = D∆t/(∆x)2 assumes a fixed constant value, where N is an odd positive integer parametrizing the alghorithm. The error between the solutions of the discrete and the continuous equations goes to zero as (∆x)2(N+2) and the values of γN are dimension independent. With these new integration methods, anisotropy effects resulting from the finite differences are minimized, defining a standard for validation and calibration of numerical solutions of diffusion and reaction-diffusion equations. Comparison between numerical and analytical solutions of reaction-diffusion equations give global discretization errors of the order of 10−6 in the sup norm. Circular patterns of travelling waves have a maximum relative random deviation from the spherical symmetry of the order of 0.2%, and the standard deviation of the fluctuations around the mean circular wave front is of the order of 10−3 . 1 2 rui@sd.ist.utl.pt jsainhas@sd.ist.utl.pt Phone: (351)-(1)-8417617, FAX: (351)-(1)-8419123 1 1 - Introduction After the important work of Turing [1952] on the chemical basis of morphogenesis, theoretical and experimental studies have shown that solutions of reaction-diffusion (RD) partial differential equations present self-organizing properties, common to biological and chemical extended systems. These properties appear as emerging coherent patterns or structures in spatially extended media. Due to the nonlinear nature of local kinetic mechanisms in RD systems, local fluctuations can grow and propagate by diffusion to the surrounding media. For this reason, extended media supporting local nonlinear kinetic mechanisms are called excitable media [Winfree, 1990; Keener & Tyson, 1992]. The prototype experimental system for the study of pattern formation in extended media are the Belousov-Zhabotinsky reaction [Zaikin & Zhabotinsky, 1970; Zhabotinsky & Zaikin, 1973] and, in biological systems, the aggregation patterns of colonies of Dictyostelium discoideum, [Tyson et al., 1989; Steinbock et al., 1993]. Due to the similarity of dynamic behaviors, these experimental systems provide insights into the mechanisms of pattern formation (morphogenesis). A large collection of experimental evidence of spatiotemporal patterns in nature can be found in Field, Körös & Noyes [1972], Meinhardt [1982], Harrison [1993], Cross & Hohenberg [1993], Kock & Meinhardt [1994] and Epstein & Showalter [1996]. Reaction-diffusion (partial differential) equations are the simplest mathematical models for the study of pattern formation in excitable media. The coherent patterns appearing in numerical simulations of RD equations are undamped travelling waves, spiral travelling waves, stable strips and spotty structures (Turing patterns), [Castets et al.,1990]. From the qualitative point of view, the same structures can be obtained with cellular automata models, [Markus & Hess, 1990; Durrett & Griffeath, 1993]. However, cellular automata models cannot be calibrated and validated from real parameters (kinetic reaction rates and diffusion coefficients), as the physical and chemical mechanisms that characterize the observed phenomena are generally continuous models. These parameters are specially important as different pattern topologies are obtained when they are varied (bifurcation phenomena). The possibility to validate a physical or chemical mechanism for the emergence of a given coherent pattern lies on the existence of reliable numerical alghorithms. If simulation models show the same qualitative behavior as real systems, model parameters can be calibrated with experimental data. If the calibrated model has the same behavior as the real system, the model is validated, and the association of specific parameter values to the real system characterizes its properties, [Oreskes et al., 1994]. A classical example of this methodology is obtained in the study of diffusion. In this case, the diffusive properties of a medium are quantified by the diffusion coefficient, measured by comparing experimental concentration profiles with the analytical solutions of the diffusion equation. Therefore, 2 diffusion, originated by the random impacts of the solvent molecules with suspended particles, is quantifiable by the diffusion coefficient, measuring the mean square displacement of the suspended particles [Chandrasekhar, 1943]. As experimental RD systems are in general nonlinear, the calibration and validation of models should be based on the comparison between the data obtained with numerical solutions of RD equations and data acquired in experiments. This rises the problem of comparison between numerical and analytical solutions of nonlinear equations — benchmarking. 1.1 - Problems with numerically obtained patterns with RD equations The prototype mathematical model for the study of pattern formation in excitable media is the two dimensional system of RD equations   2 ∂ 2ϕ ~ ∂ ϕ ~ ∂ϕ ~ T ~ ~ (1.1) = X(~ ϕ) + D · + ∂t ∂x2 ∂y 2 where ϕ ~ = (ϕ1 , ϕ2 ) are the concentrations of some chemical species or morphogens [Tur~ T = (D1 , D2 ) are diffusion coefficients and X ~ T = (X1 , X2 ) is a vector field ing, 1952], D describing the local kinetics of the system. Systems of partial differential equations of this type describe, for example, the space-time pattern formation in the Belousov-Zhabotinsky reaction and the interaction between cells separated by permeable membranes, [Turing, 1952; Zhabotinsky & Zaikin, 1973; Murray, 1993]. In general, for non-linear vector fields, there are no general analytical techniques to obtain solutions of (1.1) and we must rely on numerical methods. Stability and consistency conditions for finite difference integration methods insure that the truncation errors of numerical solutions vanish when the discrete space and time steps go simultaneously to zero [Smith, 1985; Sewell, 1988]. However, for large integration times the accumulated error can grow indefinitely [John, 1971], distorting concentration patterns. The simplest finite difference explicit method used in numerical integration of eq. (1.1) is the system of difference equations
D1 ∆t t t t t t ϕ + ϕ + ϕ + ϕ − 4 ϕ 1,i−1,j 1,i+1,j 1,i,j−1 1,i,j+1 1,i,j (∆x)2
D2 ∆t t t t t t t t ϕt+∆t ϕ + ϕ + ϕ + ϕ − 4 ϕ 2,i−1,j 2,i+1,j 2,i,j−1 2,i,j+1 2,i,j 2,i,j = ∆tX2 (ϕ1,i,j , ϕ2,i,j ) + (∆x)2 (1.2) where the subscripts (i, j) refer to the coordinates (i∆x, j∆x) of the vertices of a squared symmetric lattice in the plane. The solutions at time t, ϕt1,i,j and ϕt2,i,j , are obtained iteratively through given initial functions ϕ1 (x, y, 0) = f1 (x, y) and ϕ2 (x, y, 0) = f2 (x, y), ϕ01,i,j = f1 (i∆x, j∆x) and ϕ02,i,j = f2 (i∆x, j∆x). The centered difference scheme (1.2) is t t ϕt+∆t 1,i,j = ∆tX1 (ϕ1,i,j , ϕ2,i,j ) + 3 stable if γ = max{D1 , D2 }∆t/(∆x)2 ≤ 1/4, and the spatial truncation error is of the order of (∆x)2 , [Smith, 1985; Sewell, 1988]. In the literature on numerical methods, it is generally accepted that for the integration of two-dimensional diffusion and RD equations, the parameter γ = D∆t/(∆x)2 must be as close as possible to the value γ = 1/4, [Carslaw & Jaeger, 1959; Crank, 1975; Smith, 1985; Sewell, 1988]. This is justified by the simple fact that computation time is faster and the choice of γ is not very important if the purpose is to determine equilibrium solutions [Press ~ = ~0, equilibrium solutions of the difference scheme (1.2) coincide et al., 1989]. In fact, if X with the equilibrium solution of (1.1). However, if the goal is to obtain information about propagation velocities of diffusion fronts, or to study the pattern formation in reactiondiffusion processes, this approach presents serious problems, as we show now. The first observation is that numerically computed concentration profiles obtained with RD equations depend strongly on the magnitude of γ = D∆t/(∆x)2 , within the stability condition γ ≤ 1/4. To show this, we take the Brusselator model, [Prigogine & Lefevre, 1968], in an extended two dimensional media, defined by the equations  2  ∂ ϕ1 ∂ 2 ϕ1 ∂ϕ1 2 = k1 A − (k2 B + k4 )ϕ1 + k3 ϕ1 ϕ2 + D1 + ∂t ∂x2 ∂y 2 (1.3)  2  ∂ ϕ2 ∂ 2 ϕ2 ∂ϕ2 2 = k2 Bϕ1 − k3 ϕ1 ϕ2 + D2 + ∂t ∂x2 ∂y 2 where the kinetic terms correspond to the system of autocatalitic reactions A →k1 ϕ1 , B + ϕ1 →k2 ϕ2 + D, 2ϕ1 + ϕ2 →k3 3ϕ1 and ϕ1 →k4 E, the ki ’s are the reaction rates, and the concentrations A and B are kept constant. System (1.3) has an equilibrium unstable solution for ϕ1 = Ak1 /k4 and ϕ2 = Bk2 k4 /(Ak1 k3 ). In Fig. 1, different evolved patterns are shown, obtained with numerical integration method (1.2), for the same initial condition on a lattice of 200×200 cells and several values of γ = max{D1 , D2 }∆t/(∆x)2 . The initial condition has been chosen to produce spiral patterns, [Klevecz et al., 1992]. As γ increases within the stability condition γ ≤ 1/4, patterns change shape and, after some threshold, the spiral disappears leading to concentric waves. For small values of γ, the actualization process induces a squared spurious symmetry in the concentration profiles, showing spatial anisotropy. These simulations show that different nonequilibrium patterns are obtained without changing either the physical parameters (diffusion and kinetic coefficients), or the initial conditions. On the other hand, p as ∆x = D∆t/γ, different lattice lengths are obtained as we vary γ, for the same number of lattice sites. Therefore, the propagation velocity and physical dimensions of patterns obtained numerically, as well as its symmetry properties, are strongly dependent on the magnitude of γ, leading to deformations in the patterns, together with quantitative changes in the propagation velocities of wave fronts. 4 On the other hand, the ad hoc scaling introduced in Fig. 1, show that space and time variables are not independent. This is a well known symmetry associated with the diffusion equation, which is invariant under the linear substitutions x′ = ax and t′ = a2 t [John, 1971]. So, we are faced with the problem of knowing which of the simulations of Fig. 1, better approximates the solution of the continuous RD equation, for the prescribed initial data. To avoid the problem of numerical diffusion anisotropy in patterns of RD systems, some authors use implicit integration schemes as the Crank-Nicholson or ADI methods, [Press et al., 1989], which are computer time consuming. With explicit integration discrete methods, Barkley [1990] proposed an explicit nine point formula difference approximation for the Laplace operator, Pearson [1995] adds random noise, and Weimar et al. [1992] analyse several “masks” for the discrete approximation to the diffusion equation. In the context of ecological modelling, Dejak et al. [1987] makes an error analysis as a function of the scaling parameter γ. For cellular automata models, Markus & Hesse [1990] use a random algorithm. All these strategies to overcome the problem of diffusion anisotropy are heuristic and estimates of global numerical errors are lacking. On the other hand, these methods do not show the role of the heuristic scaling relation introduced in depicting the simulations of Fig. 1, which otherwise change the physical dimensions of the system. For a general discussion on “plausible-looking results which are qualitatively incorrect” see Ruuth [1995]. Another problem when simulating diffusion processes with explicit methods is related to the maximum propagation velocity of “suspended particles”. For the discrete method (1.2), in one time step, local concentration changes are due to next neighbors contributions and the maximum propagation velocity is vmax = 2 ∆x ∆t . But, for the diffusion equation, vmax = +∞, as is well know. The solution of the discrete explicit method (1.2) converges to the solution of Eq. (1.1), if lim∆t→0 lim∆x→0 γ =constant [John, 1971]. But as γ = D∆t/(∆x)2 , this last condition implies that lim∆t→0 lim∆x→0 vmax = +∞, but for finite ∆x and ∆t, vmax remains finite. Therefore, the best choice of γ remained an open question. According to the microscopic interpretation of diffusion by the theory of Brownian motion, the value of γ should be, γ = 1/2, γ = 1/4 and γ = 1/6 in one, two and three dimensions, respectively [Murray, 1993]. So, it seems natural to choose one of the above values of γ in simulations. However, these values correspond to the upper limit of stability of explicit methods in one, two and three dimensions, respectively, and introduce numerical instabilities due to roundoff effects. 5 1.2 - An overview of results and organization of the paper In this paper, we derive a new class of explicit difference methods for the numerical integration of diffusion and reaction-diffusion equations and a rigorous criterion for the choice of the space-time scaling parameter γ = D∆t/(∆x)2 in order that: a) Explicit difference methods for the numerical integration of RD equations give accurate results for the transient solutions of diffusion and reaction-diffusion processes, enabling to calibrate model parameters with experimental data. b) The relation between space and time scales, lim∆t→0 lim∆x→0 γ =constant, a property of diffusion equation, is correctly interpreted. c) The spatial distribution of concentrations do not depend on the symmetries of the discretization lattice. In particular, diffusion anisotropy in dimensions two and three are minimized, without additional hypotheses or data modification. The calibration of the space and time scales results from a choice of a constant value for γ, dependent on the order of the global discretization error: γ ≡ γN , for an error of the order of (∆x)2(N+2) , where N ≥ 1 is an odd integer. The main consequences of the approach developed here are: i) There exists an optimum value of γ, γ = γN , for which the new class of explicit difference methods shows the smallest global error, when its solutions are compared with the solutions of the continuous diffusion and RD equation — Fig. 2. The value of γN is dimension independent, defining a space-time scaling relation. ii) The algorithm is explicit and the accuracy of the simulations can be arbitrarily increased, for finite ∆x and ∆t — Eqs. (2.6), (2.34) and (2.47). iii) Numerical lattice anisotropy in diffusion and reaction-diffusion equations are minimized — Fig. 5. iv) In the context of the probabilistic interpretation of diffusion through Brownian motion, the transition probabilities from one lattice cell to the adjacent cells is exactly calculable, being dependent on the number of neighborhood cells taken into account in the discretization. The values of the transition probabilities determine the optimum Laplacian mask operator [Weimar et al., 1992], as a function of neighborhood order. v) Solutions of RD equations with cylindrical or spherical symmetries numerically coincide, independently of the space dimension of simulations — Fig. 8. In the next section we introduce the main technique for the discretization of the diffusion equation in dimensions 1, 2 and 3, in order to account for lattice anisotropies. The main results of the paper are stated in Theorems A, B and C, and the exact solutions of a linear RD equation are compared with the numerically calculated solution — benchmarking. In Sec. 3 we summarize the main results of the paper from the applied point of view. 6 2 - Main Results 2.1 - One space dimension We begin with the diffusion equation in one space dimension ∂2φ ∂φ =D 2 ∂t ∂x (2.1) and we take a lattice Σ in the (x, t) slab with x ∈ R and 0 ≤ t ≤ T . With ∆x and ∆t as the increments in the variables x and t, the lattice Σ has coordinates (n∆x, k∆t), with n ∈ Z and k ∈ N. If φ(x, t) is a solution of the parabolic equation (2.1), satisfying the initial condition φ(x, 0) = f (x), this solution, restricted to Σ, will be denoted by φkn = φ(n∆x, k∆t). As it is well known, for the initial condition f (x), the solution of the diffusion equation (2.1) is Z +∞ (x−y)2 1 (2.2) f (y) √ φ(x, t) = e− 4Dt dy 4πDt −∞ If f (x) is a discontinuous and integrable, the function φ(x, t) is of class C ∞ (R), for every t > 0. So, without loss of generality, we consider that φ(x, t) is of class C ∞ (R). On the other hand, if f (x) has compact support, φ(x, t) 6= 0 for every t > 0 and finite x. This implies that diffusive propagation has infinite velocity. In order to discretize the diffusion equation (2.1) on the lattice Σ, we introduce the usual finite difference approximation to the space and time derivatives, φ(x, t + ∆t) − φ(x, t) ∂φ ≃ ∂t ∆t ∂ 2φ φ(x − i∆x, t) + φ(x + i∆x, t) − 2φ(x, t) ≃ ∂x2 i2 (∆x)2 (2.3) where i is some positive integer. To account for large propagation velocities in the discretization of the diffusion equation, we rewrite (2.1) in the form N X ∂φ ∂ 2φ =D αi 2 ∂t ∂x i=1 (2.4) where the αi are non negative constants such that N X αi = 1 (2.5) i=1 Introducing the finite difference approximation (2.3) into (2.4), we obtain the discrete equation N X
αi k k k+1 k vn−i + vn+i − 2vnk (2.6) vn = vn + γN 2 i i=1 7 where γN = D ∆t (∆x)2 (2.7) The integer N is the space width of the finite difference equation (2.6) and the constants αi measure the connectivity strength between each cell and its neighbors, up to width N . With vn0 = φ0n = f (n∆x), all the vnk , for k ≥ 0, are determined recursively. For N = 1, the relation between the recursive solutions of (2.6) and the solutions of the diffusion equation for the same initial function f (x) = φ(x, 0) is easily obtained through the maximum principle. As a matter of fact, if γN is a sufficiently small constant, then |||v(x, t)||| ≤ ||f (x)|| (2.8) where ||f (x)|| = supx |f (x)| and |||v(x, t)||| = supx,t |v(x, t)| is the sup norm taken for x ∈ R and 0 ≤ t ≤ T . Under these conditions, if ∆x, ∆t → 0, vnk converges uniformly to the solution of the diffusion equation. In the case N = 1, the uniform convergence is achieved under the stability condition γ1 = D∆t/(∆x)2 ≤ 1/2, (Lemma II, §7.2 of [John, 1971]). For the general case of N ≥ 1, we apply Schwartz inequality to (2.6): ||vnk+1 || So, if 0 < γN PN αn n=1 n2 ≤ ≤ 1 2 ||vnk ||.|1 − 2γN N X αi i=1 i2 |+ ||vnk ||.|2γN N X αi i=1 i2 | and the constants αi are non negative, we have ||vnk+1 || ≤ ||vnk || (2.9) With ||vn0 || ≤ ||f (x)|| and iterating (2.9), we obtain inequality (2.8). Therefore, the finite difference method (2.6) is stable and obeys the maximum principle if γN = D 1 ∆t ∗ ≤ γN = PN 2 (∆x) 2 n=1 αn n2 (2.10) and αi ≥ 0, for i = 1, . . . , N . By the Lax Equivalence Theorem, stability of (2.6) is a necessary and sufficient condition for the initial-value problem of the diffusion equation to be properly posed [Richtmeyer & Morton, 1967]. However, from the numerical point of view there is no control on the discrete increments ∆x and ∆t, or on the values of the parameter γN . The next theorems solve these problems showing that convergence of the solutions of finite difference diffusion equation to the solutions of the diffusion equation is obtained in the limit N → ∞, for finite ∆x and ∆t. To be more precise, we denote by v N (x, t) the solution of the difference equation (2.6) on the lattice Σ, for some choice of the width N . Our main objective is to establish the conditions under which limN→∞ v N (x, t) = φ(x, t), for finite ∆x and ∆t. 8 Theorem A. If N is an odd positive integer, then there exists a positive constant γN , and real constants α1 , . . . , αN , solutions of the system of equations N i−1 γN 2 X αn n2(i−1) = 0 , − i! (2i)! n=1 i = 1, . . . , N N N X 2 γN αn n2N = 0 − (N + 1)! (2(N + 1))! n=1 (2.11a) (2.11b) ∗ ∗ Moreover, if α1 , . . . , αN are non negative and γN ≤ γN , with γN given by (2.10), the solution of the difference equation (2.6) approaches the solution of the continuous equation (2.1) with an error of the order of (∆x)2(N+2) . The values of αi as a function of γN are obtained solving the system of linear inhomogeneous equations (2.11a) and γN is largest real root of the polynomial (2.11b). Remarks on Theorem A: As a matter of fact, we have tested the positivity of the ∗ constants αi , as well as the stability condition γN < γN , up to N = 23, for odd N . Numerical tests indicate that the constants αi are positive for N odd, only if γN is the largest real root of (2.11b). This suggests the stronger assertion that, limn→∞ v 2n+1 (x, t) = φ(x, t). For N even the theorem is false. Proof: Let φ(x, t) be a solution of the diffusion equation (2.1). We suppose that φ(x, t) is of class C ∞ (R) in both variables x and t. We define the difference operator N X αn Λφ = φ(x, t + ∆t) − φ(x, t) − γN (φ(x + n∆x, t) + φ(x − n∆x, t) − 2φ(x, t)) (2.12) 2 n n=1 PN where n=1 αn = 1 and γN = D∆t/(∆x)2 . Applying the Taylor theorem to the operator Λφ, we obtain, ! N N i−1 X γN ∂ 2i φ 2 X 2i 2(i−1) Λφ = (∆x) γN − αn n ∂x2i i! (2i)! n=1 i=1 ! N   N X 2 ∂ 2(N+1) φ γ 2(N+2) 2N 2(N+1) N + O (∆x) α n − + (∆x) γ n N (N + 1)! (2(N + 1))! n=1 ∂x2(N+1) (2.13) where we have used the relation, 2i ∂ iφ i∂ φ = D ∂ti ∂x2i Choosing γN and α1 , . . . , αn such that N i−1 γN 2 X αn n2(i−1) = 0 , − i! (2i)! n=1 i = 1, . . . , N N N X γN 2 αn n2N = 0 − (N + 1)! (2(N + 1))! n=1 9 (2.14a) (2.14b) 
we have, for large N and sufficient small ∆x, |Λφ| = O (∆x)2(N+2) . Let vik be a solution of the finite difference equation (2.6), and let ΛN v be the operator defined by N X
αn k k+1 k k N vi+n + vi−n − 2vik (2.15) Λ v = vi − vi − γN 2 n n=1 As, ΛN v = 0, we have   ||Λφ − ΛN v|| = O (∆x)2(N+2) (2.16) With, eki = φ(i∆x, k∆t) − vik and subtracting (2.15) from (2.12), Λφ − ΛN v = ek+1 − eki − γN i N X
αn k ei+n + eki−n − 2eki 2 n n=1 = Λφ − ΛN v + eki + γN ek+1 i N X
αn k ei+n + eki−n − 2eki 2 n n=1 and Therefore, by the Schwartz inequality, || = ||φ(x, t + ∆t) − v N (x, t + ∆t)|| ≤ ||Λφ − ΛN v|| + ||φ(x, t) − v N (x, t)|| ||ek+1 i (2.17) ∗ on the lattice Σ, provided αi ≥ 0, for i = 1, . . . , N , and 0 < γN ≤ γN . As, ||φ(x, 0) − N v (x, 0)|| = 0 in Σ, and iterating (2.17), |||φ(x, t) − v N (x, t)||| ≤    T  O (∆x)2(N+2) = O (∆x)2(N+2) ∆t (2.18) for all (x, t) ∈ Σ. To finish the proof, we investigate the conditions under which equations (2.14) have real and positive solutions. The system of N equations (2.14a) can be written in the form ~ A~ α=L (2.19) ~ have = (j 2 )i−1 , and the vectors α ~ and L where A is a Vandermond matrix with aij = λi−1 j Q i−1 components αi and (2i)!γN /2(i!), respectively. So, as det A = i<j (λj − λi ) 6= 0, the equation (2.19) has solutions ai ≡ ai (γN ), for every N > 0, where ai (γN ) are polynomials in γN . Introducing the polynomials ai (γN ) into (2.14b), and as the αi (γN ) have degree N − 1 in γN , (2.14b) has a real root if N is odd. We now prove that (2.14b) has a positive real root. With (2.14a) for i = 2, . . . , N and (2.14b), we construct the linear equation ~′ B~ α=L 10 (2.20) ~ ′ has components where B is a Vandermond matrix with elements bij = λij = (j 2 )i , and L i (2i + 2)!γN /2((i + 1)!). As det B 6= 0, the solutions of equation (2.20) are polynomials in γn of degree N . We denote these polynomial by α′i (γN ). Introducing these solutions into (2.14a), for i = 1, we have N X a′i (γn ) − 1 = 0 (2.21) p(γn ) := i=1 As det B 6= 0, solving (2.20) for γn = 0, we have a′i (0) = 0, and p(0) = −1. Therefore, if N is odd, p(γN ) has a positive real root. Thus, Theorem A is proved. In Tab. 1 we present the values of γN and αi , for N = 1, 3 and 5, determined by ∗ Theorem A, as well as the value of the parameter γN , defining the stability limit of the difference equation (2.6). Theorem A has several important consequences for the numerical simulation of diffusion equations, namely: 1) The constant γN = D∆t/(∆x)2 introduces a coupling between the space and time scales, showing that in the finite differences approximation to the diffusion equation the space and time scales are not independent, for N odd and finite ∆x and ∆t. As γN is a constant, it defines a scaling relation for the discrete diffusion processes. This scaling relation is the best one in the sense of error minimization in the sup norm. 2) The global error between the numerical and exact solutions of the diffusion equations obtained by the finite difference method (2.6) depend only on the space increment ∆x and is time independent. Therefore, choosing N , the accuracy of the numerical solutions can be increased by decreasing ∆x, and as γN is a constant, this leads to the decrease of the time step according to the relation ∆t = γN (∆x)2 /D. This property will be analysed in more detailed below. On the other hand, the increase of accuracy can be obtained, for fixed ∆x and ∆t, increasing the number of neighborhood cells in (2.6). 3) If the αi are positive constants, we can interpret the connectivity coefficients αi in (2.6) as being proportional to the transition probability of a particle to jump from a cell to its neighborhood cells. More precisely, if Pi−>j is the transition probability of a particle initially at [i∆x − ∆x/2, i∆x + ∆x/2] to be at [j∆x − ∆x/2, j∆x + ∆x/2], after a time ∆t, by (2.6), we have Pi−>i = 1 − 2γN ∗ γN , N X αi i=1 i2 =1− γN ∗ , γN Pi−>j = γN α|i−j| |i − j|2 with |i − j| ≤ N (2.22) if γN ≤ Tab. 2. The above numerical method can be easily extended for systems of nonlinear RD partial differential equations. To be more specific we take the RD equation ∂2φ ∂φ = D 2 + f (φ) ∂t ∂x 11 (2.23) where f (φ) is a continuous function. The finite differences approximation to (2.23) can be written as N X
αi k k k vnk+1 = vnk + γN v + v − 2v + ∆tf (vnk ) (2.24) n−i n+i n 2 i i=1 P where, as before, αi = 1. By Theorem A, as γN is a scaling constant, ∆t = γN (∆x)2 /D and the error in the sup norm between the exact and approximate solutions of (2.23) and (2.24) is of the order of (∆x)2 . In applications, the finite difference approximation (2.24) is sufficient to obtain qualitative simulation results of RD equations. On the other hand, due to the characteristic finite scales in biological systems, some authors consider (2.24) as the basic RD equation [Turing, 1952]. However, in chemical systems we encounter stiff local kinetics. In some of this cases, the Euler type approximation (2.24) remains convergent. In order to compare the solutions of the diffusion and RD equations (2.1) and (2.23) with the solutions of the difference equations (2.6) and (2.24), we take the simplest integrable reaction-diffusion system obtained by coupling the pseudo-first order reaction A + X →β 2X with diffusion, [Tilden, 1974], where X and A represent two chemical species and β is the reaction rate. Representing the chemical species and its concentration by the same symbol, the RD equation for the pseudo-first order reaction in infinitely extended media is ∂X ∂2X = D 2 + βAX (2.25) ∂t ∂x where A and β are kept constants. We consider now the one dimensional lattice with coordinates (i∆x). The initial distribution of X, localized at x = 0, is represented by the function X(x, t = 0) = ∆xX0 δ(x), whose mean value in the interval [−∆x/2, ∆x/2] is X0 . So, the solution of (2.25) with this initial distribution is eβAt −x2 /4Dt e X(x, t) = ∆xX0 √ 4πDt (2.26) We now approximate the reaction-diffusion equation (2.25) by the discrete evolution equation N X
αi k k+1 k k k vn = vn + γN v + v − 2v + ∆tβAvnk (2.27) n−i n+i n 2 i i=1 where γN = D∆t/(∆x)2 and the constants αi are given by Theorem A, Tab. 1. Under these conditions, γN is a constant and ∆t is determined by the relation ∆t = γN (∆x)2 /D. The initial conditions corresponding to the localized distribution of X are then v00 = X0 , vi0 = 0 for i 6= 0 Note that if X(x) is an initial distribution of matter and we take the class of funcR R tions YX = {Y (x) : i∆x±∆x/2 Y (x)dx = i∆x±∆x/2 X(x)dx, for all i ∈ Z}, the numerical 12 solution of the diffusion equation at time t is the same for any Y ∈ YX . Therefore, we define ǫ = supx∈[i∆x−∆x/2,i∆x+∆x/2] |X(x) − X̄i |, where X̄i is the mean value of X(x) in the interval [i∆x − ∆x/2, i∆x + ∆x/2], and ε is i independent. The constant ε is the incertitude in the initial function X(x) associated to the discretization step ∆x. In Fig. 2 we show the global error, calculated with the norm of the supremum, between the solution (2.26) and the solution of (2.27), as a function of γN , for N = 1 and P N = 3. In the case N = 3, and in order to obey the relation αi = 1, the values of αi , i = 1, 2, 3, are function of γ3 , obtained by solving (2.11a): a1 = 3/2 − 13γ3 /4 + 5γ32 /2, a2 = −3/5 + 4γ3 − 4γ32 and a3 = 1/10 − 3γ3 /4 + 3γ32 /2. We now estimate the computing time needed to calculate the solution of the discrete diffusion equation as a function of the number of neighbors taken into account into (2.6) and of the fixed time step ∆t, for a finite lattice with zero flux boundary conditions. If L is the number of integer sites of the lattice, the number of memory calls to calculate the concentration at time t is, disregarding boundary conditions, TNdif ≃ L(2N + 1)Integer Dt t = L(2N + 1)Integer ∆t γN (∆x)2 (2.28) To decrease the global error four orders of magnitude, the increase in memory calls, or relative excess computing time, is dif ≃ T(N+2)/N dif TN+2 = TNdif 2N + 3 γN 2N + 1 γN+2 (2.29) dif dif dif = 0.88. Therefore, = 0.84 and T7/5 = 0.79, T5/3 With the values of Tab. 1, we have, T3/1 the computing time is faster with increasing N , with better accuracy in the global error. However, the incertitude ε in the initial function X(x), introduced by the discretization step ∆x, increases. This is a consequence of the scaling relating γN =constant. However, for RD equations the global error can be larger for larger N . In fact, the error associated to (2.24) is of the order of ∆t = γN (∆x)2 /D and γN+2 > γN . To decrease the global error in the numerical integration of RD equations with increasing N , we must therefore ask that, ∆N+2 t ≤ ∆N t, where ∆N t is the time step taken in (2.24) with N neighbors. This last condition implies that γN+2 (∆N+2 x)2 ≤ γN (∆N x)2 . With L(N+2) (∆N+2 x) = L(N) (∆N x), where L(N) is the number of integer sites of the lattice for fixed N , the relative excess computing time to decrease the global error four orders in the discrete RD equation is, by (2.28), r 2N + 3 γN+2 (L(N+2) )3 (2N + 3)γN RD (2.30) ≥ T(N+2)/N = (N) 3 2N + 1 γN (L ) (2N + 1)γN+2 and ∆ (N+2) x≤∆ (N) 13 x r γN γN+2 (2.31) So, as γN+2 > γN , the global error in RD systems decreases with increasing N if the lattice space increment decreases according to (2.31). With the values in Tab. 1, we have, ∆(3) x ≤ 0.69∆(1) x, ∆(5) x ≤ 0.81∆(3) x, ∆(7) x ≤ 0.86∆(5) x. In this case, computing time RD RD RD is slower and we have, T3/1 = 2.41, T5/3 = 1.59 and T7/5 = 1.37. Therefore, in order to decrease the global numerical error in the integration of RD equations, it is sufficient to increase the number of neighborhood connectivities and simultaneously decrease the space lattice increment ∆x according to the relation (2.31). For example, in the simulations of Fig. 2, the number of iterations k to reach the time t = 0.002 varies with γN , k = Integer(t/∆t) = Integer(Dt/γN (∆x)2 ). For N = 1 and N = 3 we have, k = 120 and k = 95, respectively. According to (2.29), the relative excess dif computing time is T3/1 = 0.79, and therefore, T3dif < T1dif . For N = 3, to decrease the global error, we should have chosen ∆x = 0.0069, according to (2.31). But, in this case, RD T3/1 = 2.41, and k = 289. 2.2 - Two space dimensions To extend the previous results for the two-dimensional diffusion and RD equation, we consider a squared lattice with spatial coordinates (i∆x, j∆x) with i, j ∈ Z. In order to enumerate lattice points in the neighborhood of the generic point (i∆x, j∆x), we associate to each integer n ≥ 1 the set of integer coordinates Jn = {(r, s) : r 2 + s2 = dn , r, s ∈ Z , dn ∈ N}, where {dn } is the sequence of positive integers such that r 2 + s2 = dn has integer solutions, {dn } = {1, 2, 4, 5, 8, 9, 10, 13, . . .}. For example, J1 = {(−1, 0), (1, 0), (0, −1), (0, 1)} J2 = {(−1, −1), (1, −1), (−1, 1), (1, 1)} (2.32) J3 = {(0, 2), (2, 0), (0, −2), (−2, 0)} J4 = {(2, 1), (1, 2), (−2, 1), (−1, 2), (2, −1), (1, −2), (−2, −1), (−1, −2)} Let hn be the number of elements in each set Jn , hn = #Jn . Writing the diffusion equation in two space dimensions as  2  M X ∂ φ ∂ 2φ ∂φ αi =D + ∂t ∂x2 ∂y 2 i=1 (2.33) PM where αi are non negative constants with n=1 αn = 1, the finite difference approximation to (2.33) is M X X
αn k+1 k k k vi+r,j+s − vi,j (2.34) vi,j = vi,j + 4γN d h n=1 n n (r,s)∈Jn where the index n refers to the neighborhood order relative to the central cell with lattice coordinates (i, j), and γN = D∆t/(∆x)2 , Fig. 3. The integer M is the space width of 14 the finite difference equation (2.34). In order to distinguish the error order index N and the neighbor order M we have introduced two different indices. The role of this indexes is specified in Theorem B below. The stability condition for (2.34) is easily derived (Sec. 2.1). The two-dimensional finite difference equation (2.34) is stable and obeys the maximum principle if γN = D 1 ∆t ∗ ≤ γN = PM 2 (∆x) 4 n=1 αn dn (2.35) and the constants αn are non negative. The solutions of the finite difference equation (2.34) approach the solutions of the two-dimensional diffusion equations (2.33) under the following conditions: Theorem B. Let R1 be the rank of the matrix of the linear system of inhomogeneous equations  X M m−1   X γN αn 4 m 2m − r 2m−2k s2k = 0 (2.36) k m! (2m)! 2k n=1 dn hn (r,s)∈Jn in the variables αn , with, 1 ≤ m ≤ N + 1, 0 ≤ k ≤ m/2, M = 1 + N + (N + 1)2 /4, and let γN be as in Theorem A, for odd N . Let R2 be the rank of the matrix of the linear system of inhomogeneous equations (2.36), with M = R1 . Then, there exist real constants α1 , . . . , αM , solutions of (2.36), with M ≥ R2 and R2 ≤ 1 + N + (N + 1)2 /4. Moreover, if M > R2 , the linear system (2.36) has an infinite number of solutions. If, for some M , ∗ ∗ α1 , . . . , αM can be chosen non negative and γN ≤ γN , where γN is given by (2.35), then the solution of the difference equation (2.34) approaches the solution of the continuous equation (2.33) with an error of the order of (∆x)2(N+2) . Remarks on Theorem B: It can happen that an error of the order of (∆x)2(N+2) occurs for an infinite set of solutions of (2.36), and to obtain the non negativity of the αn , we must increase M . For example, if N = 3, R1 = 7 and M = 7, and an > 0 for n = 1, . . . , 7. For N = 5, R1 = 14 and non negative solutions for the αn ’s is obtained only if M = 15. In this case, it is possible to show that an > 0, for n = 1, . . . , 15 and α15 ∈ [3.7404 × 10−6 , 4.747 × 10−6 ]. We have tested numerically the choice of positive constants αn up to N = 7. Proof: Let φ(x, y, t) be a solution of the two-dimensional diffusion equation (2.33). We suppose that φ(x, y, t) is analytic in x, y and t. Let us define the difference operator M X αn Λφ = φ(x, y, t+∆t)−φ(x, y, t)−4γN d h n=1 n n 15 X (r,s)∈Jn (φ(x + r∆x, y + s∆x, t) − φ(x, y, t)) (2.37) where PM n=1 αn = 1 and γN = D∆t/(∆x)2 . By the Taylor expansion,   +∞ X m X ∂ m φ(x, y, t) 1 m m−k k ǫ1 ǫ2 φ(x + ǫ1 , y + ǫ2 , t) = φ(x, y, t) + m−k ∂y k m! k ∂x m=1 (2.38) k=0 Now suppose that (r ′ , s′ ) ∈ Jn . By definition of Jn we have that (−r ′ , s′ ) ∈ Jn and P (r ′ , −s′ ) ∈ Jn . So, with ǫ1 = r∆x and ǫ2 = s∆x, if m is odd, (r,s)∈Jn r m−k sk = 0 and we have X (φ(x + r∆x, y + s∆x, t) − φ(x, y, t)) (r,s)∈Jn     2m +∞ X 2m 2m X X 2m ∂ φ(x, y, t)  (∆x) r 2m−k sk  = 2m−k k k (2m)! ∂x ∂y m=1 k=0 (r,s)∈Jn Analogously, if k is odd, r,s∈Jn r 2m−k sk = 0, and X (φ(x + r∆x, y + s∆x, t) − φ(x, y, t)) P (r,s)∈Jn = +∞ X m X m=1 k=0     2m 2m ∂ φ(x, y, t)  X 2m−2k 2k  (∆x) r s (2m)! 2k ∂x2m−2k ∂y 2k 2m (2.39) (r,s)∈Jn Taking time derivatives of the two-dimensional diffusion equation, m   X ∂ 2m φ m ∂ mφ m = D k ∂x2m−2k ∂y 2k ∂tm (2.40) k=0 introducing (2.39) and (2.40) into (2.37), and truncating the expansion for m = N + 1, we obtain for the operator Λφ,      X M N+1 m m 2m X X X αn ∂ φ γN m − 4γN 2m r 2m−2k s2k  Λφ = (∆x)2m 2m−2k 2k k 2k ∂x ∂y m! (2m)! d h n=1 n n (r,s)∈Jn m=1 k=0   + O (∆x)2(N+2) (2.41) Introducing the operator M Λ v= k+1 vi,j − k vi,j M X αn − 4γN d h n=1 n n X (r,s)∈Jn k k vi+r,j+s − vi,j

in order to make |Λφ − ΛM v| = O (∆x)2(N+2) , we must have  X M m−1   αn 4 m 2m γN − k m! (2m)! 2k n=1 dn hn 16 X (r,s)∈Jn r 2m−2k s2k = 0 (2.42) with m = 1, . . . , N + 1 and k = 0, . . . , m. For fixed N , we have ((N + 2)2 + (N + 2))/2 − 1 equations and M unknown constants α1 , . . . , αM . However, the equations in (2.42) are not independent. Dependent equations in (2.42) occur for values of k, say k ′ and k ′′ ,


2m
P m 2m 2m−2k′ 2k′ s = such that k ′ = m − k ′′ , because m (r,s)∈Jn r k′ = k′′ , 2k′ = 2k′′ and P 2m−2k′′ 2k′′ s . Therefore, we can restrict the range of the indix k in (2.42) to the (r,s)∈Jn r integers k ≤ m/2, obtaining (2.36).
Let bm be the cardinality of the set { m k : 0 ≤ k ≤ m/2}. By the Pascal triangle rule, it is straightforward to show by induction that, bm+2 = bm + 2, with b1 = 1 and b2 = b3 = 2. Therefore, for each m, the number of equations in (2.42), with k ≤ m/2, is bm = 3/4 + (−1)m /4 + m/2. Now let us denote by cm the number of equations in (2.42) up to order m and with k ≤ m/2. So, we have the recurrence relation cm+1 = cm + bm+1 , with c1 = 1. The solution of this recurrence is cm = −1/8 + (−1)m /8 + m + m2 /4, and the number of equations in (2.42), with k ≤ m/2, is 1 1 1 M1 = − + (−1)N+1 + N + 1 + (N + 1)2 8 8 4 Let R1 be the rank of the matrix of the linear inhomogeneous system (2.36), for N odd and M = M1 = N + 1 + (N + 1)2 /4. Let γN be as in Theorem A. If R1 = M1 , then M = M1 = N + 1 + (N + 1)2 /4, and the linear inhomogeneous system (2.36) has a solution. If, R1 < M1 , let R2 be the rank of the matrix of the linear inhomogeneous system (2.36), for N odd and M = R1 . Solving a linearly independent subsystem of (2.36) of dimension R2 × R2 , for α1 , . . . , αR2 , we obtain a solution that annulates R2 equations in (2.36). If we let M > R2 this linearly independent subsystem of (2.36) has an infinite number of solutions, dependent of M − R2 variables αn . But, in this case (R1 < M1 ), these solutions also annulate the remaining dependent equations in (2.36). Otherwise, as ′ the αn ’s do not depend on ϕ, |Λϕ| = O(∆x)2(N +2) , with N ′ < N , and this contradicts Theorem A for the choice ϕ(x, y, t) = ϕ(x, t). Finally, if αi ≥ 0, for 1 ≤ i ≤ M and some M ≥ R2 , by the maximum principle, the solution of the difference equation (2.34) approaches the solution of the continuous equation (2.33) with an error of the order of (∆x)2(N+2) . Therefore, Theorem B is proved. We show in Tab. 3 the values of the constants αi and γN for N = 1 and N = 3 (M = 3 and M = 7), for the discrete diffusion equation (2.34), as well as the values of stability limit (2.35). To obtain comparable results in dimensions 1 and 2, from the point of view of global errors, we must take the corresponding difference methods for the same value of N . (Note that in dimension 1, M = N ). The scaling constant is the same for the same value of N , but the connectivities with neighborhood sites must follow the order defined by Theorems A and B. 17 The connectivity relations determined by the above theorem have a simple geometric meaning. As a matter of fact, for each N , the local connectivities with neighborhood sites, with αi > 0 in (2.34), follow a spherically symmetric front, Fig. 3, suggesting the possibility of elimination of lattice symmetries in numerical solutions of RD equations. To test the symmetry properties of the solutions of the discrete methods (2.34), we take as prototype model the Brusselator, (1.3). By (2.34), the difference RD equation for the Brusselator is ϕk+1 1,i,j =ϕk1,i,j M X αn + 4γN d h n=1 n n X (r,s)∈Jn ϕk1,i+r,j+s − ϕk1,i,j
+ ∆t(k1 A − (k2 B + k4 )ϕk1,i,j + k3 (ϕk1,i,j )2 ϕk2,i,j ) k ϕk+1 2,i,j =ϕ2,i,j + M D2 X αn + 4γN D1 n=1 dn hn ∆t(k2 Bϕk1,i,j − X (r,s)∈Jn ϕk2,i+r,j+s − ϕk2,i,j k3 (ϕk1,i,j )2 ϕk2,i,j ) (2.43)
where γN = Dmax ∆t/(∆x)2 and Dmax = max{D1 , D2 }. In general, for systems of RD equations, the integration algorithm is a set of equations of type (2.43), one equation for each diffusive variable ϕi . The consistency with the scaling relation γn =constant, given by Theorem B, shows that we can only have D2 = D1 or D2 = 0, in agreement, respectively, with the simulation and model choices of Keener & Tyson [1992] and Zhabotinsky & Zaikin [1973]. To maintain compatibility with these limits, the diffusion terms are weighted by the non dimensional factors Di /Dmax , where p Dmax = max{D1 , D2 }. The space scale is now determined by ∆x = Dmax ∆t/γN , as otherwise we would be simulating diffusion below the mean free path of one of the diffusive variables. Under these conditions and for the same parameter values as in Fig. 1, we choose N = 1 (M = 3) in (2.43) and we have followed the time evolution of ϕ1 and ϕ2 with (2.43). In Fig. 4a) we represent the spatial distribution of ϕ1 at the time t = 435, for ∆t = 1/6, in a square lattice of 525 × 525 sites, zero flux boundary conditions, and initial conditions ϕ01,i,j = Ak1 /k4 , ϕ02,i,j = Bk2 k4 /(Ak1 k3 ), ϕ01,262,262 = 1.1(Ak1 /k4 ), for (i, j) 6= (262, 262) ϕ02,262,262 = 1.1(Bk2 k4 /(Ak1 k3 )) (2.44) The system develops undamped spherically symmetric travelling waves, originated at the perturbed lattice point (i, j) = (262, 262). To analyse the radial symmetry of wave fronts, we measured the distance from a point in the wave front to the point where the perturbation has been initiated, as a function of the polar angle θ. We denote this distance by rθ . In this case, the reference point at the wave front corresponds to a radial local maxima of 18 ϕ1 . The normalized distance from the wave front to the center, rθ /r̄, as a function of the polar angle θ, is represented in Fig. 5. From this numerical simulation we conclude that the maximum relative deviation from spherical symmetry of the wave fronts oscillates randomly in θ with a maximum relative error of the order of 0.2%. The relative deviation as a function of θ oscillates randomly, and the standard deviation of the fluctuations around the mean circular wave front is 1.2 × 10−3 . Numerical analysis for linear RD equations show the same type of behavior for the global error as in Fig. 2. Other patterns with spiral symmetry that appear in experimental systems are plotted in Figs. 2.3b-d, [Agladze & Krinsky, 1982; Ross et al., 1988]. 2.3 - Three space dimensions For the three-dimensional diffusion equation, we consider a cubic lattice with spatial coordinates (i∆x, j∆x, k∆x) with i, j, k ∈ Z. We associate to each integer n ≥ 1 the set of integer coordinates Jn = {(q, r, s) : q 2 + r 2 + s2 = dn , q, r, s ∈ Z , dn ∈ N}, where {dn } is the sequence of positive integers such that q 2 + r 2 + s2 = dn has integer solutions, {dn } = {1, 2, 3, 4, 5, 6, 8, 9, 10, 11, . . .}. For example, J1 = {(−1, 0, 0), (0, −1, 0), (0, 0, −1), (0, 0, 1), (0, 1, 0), (1, 0, 0)} J2 = {(−1, −1, 0), (1, −1, 0), (−1, 1, 0), (1, 1, 0), (1, 0, 1), (−1, 0, 1), (2.45) (1, 0, −1), (−1, 0, −1), (0, 1, 1), (0, −1, 1), (0, 1, −1), (0, −1, −1)} Writing the diffusion equation in three space dimensions as M X ∂φ αi =D ∂t i=1  ∂ 2φ ∂ 2φ ∂ 2φ + + 2 ∂x2 ∂y 2 ∂z  (2.46) PM where n=1 αn = 1 and M is the space width of the connectivity constants αi , the finite difference approximation to (2.46) is m+1 vi,j,k = m vi,j,k M X αn + 6γN d h n=1 n n X (q,r,s)∈Jn m m vi+q,j+r,k+s − vi,j,k
(2.47) where, hn = #Jn , the index n refers to the order of neighborhood of each cell and γN = D∆t/(∆x)2 , Fig. 6. The three-dimensional finite difference equation (2.47) is stable and obeys the maximum principle if γN = D 1 ∆t ∗ ≤ γ = P N M (∆x)2 6 n=1 αn dn (2.48) The solutions of the finite difference equation (2.47) approache the solutions of the three-dimensional diffusion equations (2.46) under the following conditions: 19 Theorem C. Let R1 be the rank of the matrix of the linear system of inhomogeneous equations  M     m−1  X γN 6 m1 2m1 X αn m 2m − q 2m−2m1 r 2m1 −2m2 s2m2 = 0 m2 2m2 n=1 dn hn m! m1 (2m)! 2m1 (q,r,s)∈Jn (2.49) in the variables αn , with, 1 ≤ m ≤ N + 1, 0 ≤ m1 ≤ m, 0 ≤ m2 ≤ m1 , M = 3 + 13N/3 + 3N 2 /2 + N 3 /6, and let γN be as in Theorem A, for odd N . Let R2 be the rank of the matrix of the linear system of inhomogeneous equations (2.49), with M = R1 . Then, there exist real constants α1 , . . . , αM , solutions of (2.49), with M ≥ R2 and R2 ≤ 3 + 13N/3 + 3N 2 /2 + N 3 /6. Moreover, if M > R2 , the linear system (2.49) has an infinite ∗ number of solutions. If, for some M , α1 , . . . , αM can be chosen non negative and γN ≤ γN , ∗ where γN is given by (2.48), then the solution of the difference equation (2.47) approaches the solution of the continuous equation (2.46) with an error of the order of (∆x)2(N+2) . The proof of Theorem C is analogous to the one of Theorem B. In Tab. 4 we show the values of the constants αi ’s, for N = 1 and N = 3. In the case of N = 3, the solution of (2.49) for positive αi can not be obtained for M = R2 = 9, and we have introduced the additional constant α10 . So, system (2.49) has an infinite set solutions, as a function of α10 . Numerical analysis shows that it is possible to choose αi > 0, with i = 1, . . . , 10, if α10 ∈ [0.01278, 0.01785], and γ3 < γ3∗ . In this case, the error between the solutions of (2.46) and (2.47) is of the order of (∆x)10 , for any choive of α10 in the interval [0.01278, 0.01785]. To test the symmetry properties of spherical wave fronts of the explicit method (2.47), for N = 1, we again consider the spatially extended Brusselator model with zero flux boundary conditions. We take as initial conditions a central perturbation to the steady state. In order to maintain the computations within the range of a personal computer, the three dimensional lattice has been chosen with physical dimensions of 175 × 175 × 175 cells, and the time increment was ∆t = 0.1. In Fig. 7 we present a three dimensional view of the concentration of the chemical species ϕ1 . In this case, as we choose the time step as independent increment, the dimension of the simulation corresponds to a cube of side p L = 175∆x = 175 D1 ∆t/γ1 ≃ 136, where D1 = 1.0 and D2 = 0. To analyse the symmetry properties of the numerically generated pattern in Fig. 7 and make the error analysis, we compare one-dimensional concentration profile extracted from one, two and three-dimensional simulations. If ϕ(x, t) is a solution of the one-dimensional RD equation, it is also a solution of the two and three-dimensional RD equations, with adapted initial conditions. Therefore, one-dimensional concentration profiles calculated numerically in the three cases should coincide. In Fig. 8 we present the three concentration profiles calculated with the finite difference methods (2.6), (2.34) and (2.47). 20 The matching of the concentration profiles shown in Fig. 8 show that the global errors between exact and numerical solutions of three-dimensional RD equations, as well as its symmetry properties, are the same as in the one and two dimensional cases analysed in Sec. 2.1 and Sec. 2.2. This fact enables to extrapolate properties of spherical and cylindrical solutions in two and three space dimensions from simpler and less computer time consuming one-dimensional simulations. 3 - Conclusions We have presented a general class of finite difference approximations to diffusion and RD partial differential equations whose solutions approach the solution of the continuous system, within a global error of any prescribed order. We have shown that the difference equations obey a scaling relation, relating space and time integration steps. This scaling relation is reminiscent from the well known invariance of the diffusion equation for transformations of the form x′ → ax and t′ → a2 t. The spurious lattice symmetries that appear in numerical simulations of RD equations have been controlled and minimized below reasonable errors. Numerical experiments, Fig. 1, show that specific patterns in RD systems are related with the temporal and spatial scales arising naturally in these systems. Therefore, any RD pattern obtained numerically without the proper calibration of the space and time scales can be qualitatively incorrect when compared with the exact solution of the continuous system of RD equations. 21 For computer calculation in single precision, the case N = 1 is in general sufficient to avoid lattice spurious symmetries. To be more specific, for the case N = 1 and global errors of the order of (∆x)6 , the numerical alghorithms for the integration and RD equations are: dimension 1: 1 k k vik+1 =vik + (vi−1 + vi+1 − 2vik ) + ∆tf (vik ) 6 dimension 2: 1 k k+1 k k k k k + vi+1,j + vi,j−1 + vi,j+1 − 4vi,j ) =vi,j + (vi−1,j vi,j 9 1 k k k k k k + (vi−1,j−1 + vi+1,j+1 + vi+1,j−1 + vi−1,j+1 − 4vi,j ) + ∆tf (vi,j ) 36 (4.1a) dimension 3: 1 k k+1 k k k k k + vi+1,j,m + vi,j−1,m + vi,j+1,m + vi,j,m−1 =vi,j,m + (vi−1,j,m vi,j,m 18 1 k k k k k + vi,j,m+1 − 6vi,j,m ) + (vi−1,j−1,m + vi+1,j−1,m + vi−1,j+1,m 36 k k k k k + vi+1,j+1,m + vi+1,j,m+1 + vi−1,j,m+1 + vi+1,j,m−1 + vi−1,j,m−1 k k k k k + vi,j+1,m+1 + vi,j−1,m+1 + vi,j+1,m−1 + vi,j−1,m−1 − 12vi,j,m ) k + ∆tf (vi,j,m ) and the space and time steps are calibrated according the scaling relation 1 D∆t = 2 (∆x) 6 (4.1b) The finite difference methods (4.1) enable to measure relevant physical quantities in numerical simulations, in order to calibrate system parameters. The strategy of calibration and validation of numerical methods leads to the possibility of analysing computer simulations as experiments with the same type of data analysis and conclusions. This is particularly important for modeling nonlinear RD systems, where analytical solution of transient processes are difficult or impossible to find. Several authors have developed different strategies to calibrate the numerical solutions of diffusion and RD equations. In particular, Dejak et al. [1987], based on the comparison of numerical and known solution of diffusion equations found the heuristic condition γ = 1/6. In the context of cellular automata simulations Weymar et al. [1992] compared the preservation of symmetries of several Laplacian mask operators as a function of neighborhood connectivities. The results presented here give the precise foundation for these heuristic techniques. This is particularly important for the simulation and bench-marking of three-dimensional systems where complex topological structures are expected to appear. 22 On the other hand, in the two-dimensional case, the coefficients in (4.1a) are similar to those obtained in the numerical analysis of elliptic problems [Boisvert, 1981] which evolves the discretization of the Laplacian operator. Acknowledgement: This work has been partially supported by the PRAXIS XXI Project PRAXIS/PCEX/P/FIS/26/96 (Portugal). J. S. is supported by the Junta Nacional de Investigação Cientı́fica grant number BD 922. 23 References Agladze, K. I. & Krinsky, V. I. [1982] “Multi-armed vortices in an active chemical medium,” Nature 296, 424-426. Barkley, D., Kness, M, & Tuckerman, L. S. [1990] “Spiral-wave dynamics in a simple model of excitable media: The transition from simple to compound rotation,” Phys. Rev. A 42, 2489-2492. Boisvert, R. F. [1981] “Families of high order accurate discretization of some elliptic problems,” SIAM J. Sci. Stat. Comput. 2, 268-284. Crank, J. [1975] The Mathematics of Diffusion (Oxford Uni. Press, Oxford). Carslaw, H. S., & Jaeger, J. C. [1959] Conduction of Heat in Solids (Oxford Uni. Press, Oxford). Castets, V., Dulos, E., Boissonade, J. & De Kepper, P. [1990] “Experimental evidence of a sustained standing Turing-type nonequilibrium chemical pattern,” Phys. Rev. Lett. 64, 2953-2956. Chandrasekhar, S. [1943] “Stochastic Problems in Physics and Astronomy,” Rev. Mod. Phys. 15, 1-89. Cross, M. C. & Hohenberg, P. C. [1993] “Pattern formation outside of equilibrium,” Rev. Mod. Phys. 65, 851-1112. Dejak, C., Lalatta, I. M., Messina, E. & Pecenik, G. [1987] “A two-dimensional diffusion model of the Venice lagoon and relative open boundary conditions,” Ecol. Model. 37, 21-45. Epstein, I. R. & Showalter, K. [1996] “Nonlinear chemical dynamics: Oscillations, patterns, and chaos,” J. Phys. Chem. 100, 13132-13147. Durrett, R. & Griffeath, D.[1993] “Asymptotic behavior of excitable cellular automata,” Experimental Mathematics 2, 183-208. Field, R. J., Körös, E. & Noyes, R. M. [1972] “Oscillations in chemical systems. II. Thorough analysis of temporal oscillation in the bromate-cerium-malonic acid system,” J. Amer. Chem. Soc. 94, 3649-3665. Harrison, L. G. [1993] Kinetic Theory of Living Patterns (Cambridge Uni. Press, Cambridge). John, F. [1971] Partial Differential Equations (Springer, Berlin). Keener, P. & Tyson, J. [1992] “The dynamics of scroll waves in excitable media,” SIAM Review 34, 1-39. 24 Klevecz, R. R., Bolen, J. and Durn, O. [1992[ “Self-organization in biological tissues: Analysis of asynchronous and synchronous periodicity, turbulence and synchronous chaos emergent in coupled chaotic arrays,” Int. J. Bifurcation and Chaos 2, 941-953. Kock, A. J. & Meinhardt, H. [1994] “Biological pattern formation: From basic mechanisms to complex structures,” Rev. Mod. Phys. 66, 1481-1507. Meinhardt, H. [1982] Models of Biological Pattern Formation (Academic Press, London). Markus, M. & Hess, B. [1990] “Isotropic cellular automaton for modelling excitable media,” Nature 347, 56-58. Murray, J. D. [1993] Mathematical Biology (Springer, Berlin). Oreskes, N., Shrader-Frechette, K. & Belitz, K. [1994] “Verification, validation, and confirmation of numerical models in the earth sciences,” Science 263, 641-646. Pearson, J. E. [1993] “Complex patterns in a simple system,” Science 261, 189-192. Press, W. H., Flannery, B. P., Teukolski, S. A. & Vetterling, W. T. [1989] Numerical Recipes (Cambridge Uni. Press, Cambridge). Prigogine, I. & Lefever, R. [1968] “Symmetry breaking instabilities in dissipative systems. II,” J. Chem. Phys. 48, 1695-1700. Richtmeyer, R. D. & Morton, K. W. [1967] Difference Methods for Initial-Value Problems (Wiley, New York). Ross, J., Müller, S. C. & Vidal, C. [1988] “Chemical waves,” Science 240, 460-465. Ruuth, S. J. [1995] “Implicit-explicit methods for reaction-diffusion problems in pattern formation,” J. Math. Biol. 34, 148-176. Sewell, G. [1988] The Numerical Solutions of Ordinary and Partial Differential Equations (Academic Press, London). Smith, G.D. [1985] Numerical Solution of Partial Differential Equations (Oxford Uni. Press, Oxford). Steinbock, O., Siegert, F., Müller, S. & Weijer, C. J. [1993] “Three-dimensional waves of excitation during Dictyostelium morphogenesis,” Proc. Natl. Acad. Sci. USA 90, 7332-7335. Tilden, J, [1974] “On the Velocity of Spatial Wave Propagation in the Belousov Reaction,” J. Chem. Phys. 60, 3349-3350. Tyson, J. J., Alexander, K. A., Manoranjan, V.S. & Murray, J. D. [1989] “Spiral waves of cyclic AMP in a model of slime mold aggregation,” Physica D 34, 193-207. Turing, A. M. [1952] “The chemical basis of morphogenesis,” Philo. Trans. Roy. Soc. Lond. Ser. B237, 5-72. 25 Winfree, A. T. [1990] “Stable particle-like solutions to the nonlinear wave equations of three-dimensional excitable media,” SIAM Review 32, 1-53. Weimar, J. R., Tyson, J. J. & Watson, L.T. [1992] “Diffusion and propagation in cellular automaton models of excitable media,” Physica D 55, 309-327. Zaikin, A. N. & Zhabotinsky, A. M. [1970] “Concentration wave propagation in twodimensional liquid-phase self-oscillating system,” Nature 225, 535-537. Zhabotinsky, A. M. & Zaikin, A. N. [1973] “Autowave processes in a distributed chemical system,” J. Theor. Biol. 40, 45-61. 26 Table Captions Table 1: Parameters for the discrete diffusion equation (2.6), for N = 1, 3 and 5, according to Theorem A. We have tested numerically the positiveness of the constants αi up to N = 23, for odd N . In this case, for example, a23 ≃ 10−20 . Other values of γN are: γ7 = 0.713841, γ9 = 0.896295 and γ11 = 1.07875. Table 2: Transition probabilities to neighborhood cells in the Brownian motion interpretation of diffusion. These transition probabilities depend on the scaling relation γN =constant and of the non negativity of the connectivity coefficients αi . Table 3: Parameters for the discrete diffusion equation (2.34), according to Theorem B, for N = 1 and N = 3. Table 4: Parameters for the discrete diffusion equation (2.47), for N = 1 and N = 3, under the conditions of Theorem C. In the case N = 3, equations (2.49) have an infinite set of solutions. The constants αn , with n = 1, . . . , 10, are positive if α10 ∈ [0.01278, 0.01785]. For N = 3, the parameters were obtained with the choice of α10 in the middle of the interval [0.01278, 0.01785]. 27 Figure Captions Figure 1: Temporal evolution of non-equilibrium travelling waves obtained with the Brusselator model (1.3), for several values of the parameter γ ≤ 1/4. The initial condition has been chosen to produce spiral patterns at γ = 0.01. The simulation parameters are D1 = 0.008 , D2 = 0.0016, ki = 1.0, A = 1.0 and B = 2.3. We represent the concentration profiles for the chemical species ϕ1 . Time evolution has been calculated on a lattice of 200 × 200 cells, with zero flux boundary conditions, the same initial condition and a time step ∆t = 0.1. To maintain the same radial velocity, the lattices have been scaled p according to the relation ∆x = (∆tDmax /γ), where Dmax = max{D1 , D2 }. Otherwise, propagation velocities would be different, and dependent on the integration space and time steps. Figure 2: Error between the solution (2.26) and the solution of (2.27), as a function of γN , for N = 1 and N = 3, with the parameter values of Tab. 1. The error function is supi |X(i∆x, k∆t) − vik |, calculated at the time t = 0.002. The simulations have been performed with β = 1.0, A = 1.0, D = 1, ∆x = 0.01, ∆t = γN (∆x)2 /D and the initial value X0 = 0.01. For N = 1, the error is in fact minimized when γ1 = 1/6 (Theorem A) and the global error is (0.01)6 = 10−12 . For N = 3, the global error is of the order of 10−20 , which is below the limit of precision of computer calculation in single precision (≃ 10−6 ). Figure 3: Space width and connectivity constants αn for the discrete version of the two dimensional diffusion equation (2.34). If follows from Theorem B that local connectivities with neighborhood sites, with αi > 0 in (2.34), follow a spherically symmetric front, suggesting the possibility of elimination of lattice symmetries in numerical solutions of RD equations. Figure 4: a) Two dimensional patterns generated by the discrete difference RD equation (2.43) with N = 1, for the Brusselator, with initial conditions (2.44) and the parameters of Fig. 1. The depicted figures have been calculated in a lattice of 525 × 525 sites, with zero flux boundary conditions. The anisotropy analysis of the wave fronts in a) is presented in Fig. 5. From c) to d) we present other patterns generated with different initial conditions, qualitatively similar to patterns found in experimental systems. Figure 5: a) Anisotropy analysis of the pattern of Fig. 4a). Numerically obtained wave fronts deviate from spherical symmetry, within a maximum relative deviation of the order of 0.2%. The standard deviation of the fluctuations around the mean circular wave front is 1.2 × 10−3 . 28 Figure 6: Space width M or neighborhood order for the discretization of the three dimensional diffusion equation (2.47), for N = 1 and N = 3. Figure 7: Three dimensional spherical symmetric travelling waves obtained the finite difference approximation (2.47) for the Brusselator model in three-dimensional extended media. The kinetic parameters in this simulation are the same as in Fig. 1, and the diffusion coefficients are D1 = 1.0 and D2 = 0. Figure 8: Two and three dimensional plane wave solutions, and one-dimensional concentration profiles of the variable ϕ1 , at t = 125, for the Brusselator model in extended media, calculated with one, two and three-dimensional RD equations. The parameters are the same as in Fig. 7, with N = 1. 29 N error 1 (∆x)6 3 (∆x)10 5 (∆x)14 ∗ γN 1/2 0.666869 0.866792 ∗ γN (< γN ) 1/6 0.348977 0.531397 α1 1 0.670287 0.454996 α2 0.308768 0.443529 α3 0.020945 0.095103 α4 0.006218 α5 0.000154 Table 1 30 N Pi−>i Pi−>i±1 1 2/3 1/6 3 0.4767 0.2339 5 Pi−>i±2 Pi−>i±3 0.0269 0.0008 Pi−>i±4 Pi−>i±5 0.386938 0.241783 0.058923 0.005615 0.000206 0.000003 Table 2 31 N M R1 R2 error 1 3 3 3 (∆x)6 3 7 7 7 (∆x)10 di hi = #Ji ∗ γN 3/10 0.455079 ∗ γN (< γN ) 1/6 0.348977 α1 2/3 0.307488 1 4 α2 1/3 0.330497 2 4 α3 0 0.162856 4 4 α4 0.153393 5 8 α5 0.023198 8 4 α6 0.006326 9 4 α7 0.016243 10 8 Table 3 32 N M R1 R2 error 1 2 3 2 (∆x)6 3 10 10 9 (∆x)10 di hi #Ji ∗ γN 1/4 0.407538 ∗ γN (< γN ) 1/6 0.348977 α1 1/3 0.098725 1 6 α2 2/3 0.401411 2 12 α3 0.107430 3 8 α4 0.115899 4 6 α5 0.078248 5 24 α6 0.125812 6 24 α7 0.031315 8 12 α8 0.021208 9 30 α9 0.004633 10 24 α10 0.015319 11 24 Table 4 33 This figure "fig1.gif" is available in "gif" format from: http://arxiv.org/ps/patt-sol/9712007v1 This figure "fig2.gif" is available in "gif" format from: http://arxiv.org/ps/patt-sol/9712007v1 This figure "fig3.gif" is available in "gif" format from: http://arxiv.org/ps/patt-sol/9712007v1 This figure "fig4.gif" is available in "gif" format from: http://arxiv.org/ps/patt-sol/9712007v1 This figure "fig5.gif" is available in "gif" format from: http://arxiv.org/ps/patt-sol/9712007v1 This figure "fig6.gif" is available in "gif" format from: http://arxiv.org/ps/patt-sol/9712007v1 This figure "fig7.gif" is available in "gif" format from: http://arxiv.org/ps/patt-sol/9712007v1 This figure "fig8.gif" is available in "gif" format from: http://arxiv.org/ps/patt-sol/9712007v1