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