Abstract
After a theory of morphogenesis in chemical cells was introduced in the 1950s, much attention had been devoted to the numerical solution of reaction-diffusion (RD) partial differential equations (PDEs). The Crank–Nicolson (CN) method has been a common second-order time-stepping procedure. However, the CN method may introduce spurious oscillations for nonsmooth data unless the time step size is sufficiently small. This article studies a nonoscillatory second-order time-stepping procedure for RD equations, called a variable-method, as a perturbation of the CN method. In each time level, the new method detects points of potential oscillations to implicitly resolve the solution there. The proposed time-stepping procedure is nonoscillatory and of a second-order temporal accuracy. Various examples are given to show effectiveness of the method. The article also performs a sensitivity analysis for the numerical solution of biological pattern forming models to conclude that the numerical solution is much more sensitive to the spatial mesh resolution than the temporal one. As the spatial resolution becomes higher for an improved accuracy, the CN method may produce spurious oscillations, while the proposed method results in stable solutions.
1. Introduction
As molecular imaging and single cell analysis is advancing our understanding of spatial processes shaping the cellular dynamics, new models of nonlinear dynamics are necessary. Originating in study of organism development, spatial pattern formation has received a large amount of research over the past decade. Among the most studied, the reaction-diffusion (RD) systems are generating patterns that have been shown to represent well morphogenesis. A theory of morphogenesis based on a RD model was initially proposed by Turing [1]. Gierer and Meinhardt [2] were the first to explore pattern formation in biological systems using the RD equation. Subsequently, several equations of RD type have been studied to understand patterning in developmental biology. Some were derived from phenomenological models (Gierer–Meinhardt) while other modeled simple reaction schemes (Schnackenberg trimolecular autocatalytic reactions model [3], Gray–Scott model [4], Brusselator model [5], chlorite-iodide-malonic acid, CIMA model [6]). Recent work on RD systems demonstrates that it can be used to understand biological patterns formation [7], while [8] and [9] reviewed RD systems can be used to investigate spatial patterning in developmental systems.
The RD model for biological pattern formation is defined as follows [10]. Let be a bounded domain in , , denote the boundary of , and for some :where , is the diffusion tensor, denotes the Laplace operator, is the outward normal derivative on the boundary , and is the reaction kinetics of the system given as
After Turing proposed a theory of morphogenesis in chemical cells in 1952 [1], much attention has been devoted to the numerical solution of RD problems of form (1); see [11–14] and references therein. Most of the numerical methods studied employed finite difference or finite element approximations for the spatial discretization, while some researchers use finite volume and collocation methods. Once the nonlinear reaction terms are treated (linearized or extrapolated), the Crank–Nicolson (CN) method can be applied as a second-order time-stepping procedure. Time-stepping procedures are required at each time step to solve a system of linear algebraic equations, which, although sparse, is compute intensive for multidimensional problems. In order to enhance efficiency of time-stepping procedures, one can adopt the alternating direction implicit (ADI) method as in [11–13, 15]. In particular, Fernandes et al. [12] introduce an ADI extrapolated CN orthogonal spline collocation method for RD problems.
ADI was invented as a perturbation of the CN method by Douglas, Peaceman, and Rachford in 1955 [16–18] and has been employed effectively for the calculation of numerical solution of various time-dependent multidimensional problems, either parabolic or hyperbolic [19, 20]. ADI reduces a multidimensional problem to multiple easy-to-solve one-dimensional problems, for an extra cost of a splitting error in , where is the time step size. However, the splitting error can be much larger than the sum of spatial and temporal discretization errors, unless the time step size is sufficiently small [21].
On the contrary, the CN method applied for nonsmooth data may introduce spurious oscillations to the numerical solution unless the time step size is sufficiently small to satisfy the maximum principle 12, which has been recognized in the original paper as well [22]. For this reason, whenever a larger time step or a higher spatial resolution is desirable/necessary, the (less accurate, first order) implicit method which is immune to oscillations has been used at least for several initial time steps with nonsmooth initial data [23]. The CN method and its perturbations (such as ADI) must be applied with care when the solution involves fast transitions or sharp edges; in particular, the time step size should be set very small, e.g., , where is the spatial grid size. In order to overcome the oscillation problem of the CN method applied for linear parabolic problems of nonsmooth data, the authors recently suggested a variable-method in which the time-stepping parameter of the conventional -method, , was determined based on local oscillatory characteristics of the solution and the data [24].
In this article, we apply the variable- method for the numerical solution of two-component nonlinear RD equations, as given in (1). The variable- method is a perturbation of the CN method which evolves the solution implicitly at points where the solution shows a certain portent of oscillations and maintains as a similar accuracy as the CN method with smooth data. The proposed method would be an adequate choice of time-stepping procedure for the numerical solution of RD partial differential equations (PDEs) when a larger time step or a higher spatial mesh resolution is desirable. We have performed a sensitivity analysis for the numerical solution of biological pattern forming models to the spatial and temporal grid sizes. It has been observed from various examples that accuracy of the numerical solution is much more sensitive to the spatial mesh resolution than the temporal one. When the spatial mesh resolution is set high for a higher accuracy, the method allows to keep the temporal resolution moderate or low. The suggested variable- method can result in a smooth/stable numerical solution by suppressing possible oscillations, unlike the CN method.
The article is organized as follows. Section 2 includes a brief review on the CN method and its spurious oscillatory behaviors, as preliminaries. In Section 3, a variable- method is presented for the numerical solution of two-component nonlinear RD equations. We adopt the successive over-relaxation (SOR) method to solve the resulting linear systems at each time level. Section 4 considers a heuristic technique for the choice of the optimal relaxation parameter for SOR. Section 5 gives numerical examples that show effectiveness of the variable- method applied to RD problems in 1D and 2D spaces. In Section 6, we summarize our experiments and present conclusions.
2. Preliminaries
In this section, we present a brief review of time-stepping procedures, for the numerical solution of linear parabolic equations of the formwhere is a diffusion coefficient and is a reaction/source term. We also consider difficulties arising when dealing with nonsmooth data (initial values, boundary conditions, and/or the source term).
2.1. The -Method: Difference Equation
Let be a rectangular domain in : . By partitioning , we obtain the space-time grid points:where , , and are prescribed positive integers and
Define the discrete domain, the set of the spatial grid points, byand denote the set of boundary grid points by and the set of interior grid points by .
Let for all functions defined in . Then, the second-order 5-point finite difference (FD) approximation of at readswhere the FD operators are defined as
For the temporal derivative , a convenient FD approximation can give
Expressing the spatial derivative by a weighted average of previous and current time values, we can formulate the -method for (3) aswhere is either or . A simple algebraic rearrangement of (10) in a vector form becomeswhere and , considered as column vectors. Popular choices of are 0, 1, and , which are, respectively, the explicit method (the forward Euler method), the implicit method (the backward Euler method), and the semi-implicit method (the Crank–Nicolson method).(i)Forward Euler method: when , algorithm (11) is stable when it satisfies
Although the explicit method is efficient for each time step, its stability condition in (12) enforces, choosing a small time step size ; it may become less efficient compared with other implicit methods. It is elementary in numerical analysis that when , the -method (10) is unconditionally stable.(ii)Crank–Nicolson method: when , (11) can be rewritten as
The CN method has been the most popular time-stepping procedure for the numerical solution of parabolic problems because it is stable and of second-order accuracy in both spatial and temporal directions. However, the CN method applied for nonsmooth data may introduce spurious oscillations to the numerical solution unless the algorithm parameters satisfy the maximum principle [22, 23]. As analyzed by the authors [24], the undesired oscillations are due to instability involved in the explicit half step of the CN method, the first term in the right side of (13). The variable -method proposed in [24] suppresses spurious oscillations, by evolving the solution implicitly at points where the solution shows a certain portent of oscillations or reduced smoothness, and maintains as a similar accuracy as the CN method with smooth data.(iii)Backward Euler method: when , algorithm (11) reads
Although the implicit method shows a first-order accuracy in the temporal direction, it never introduces spurious oscillations to its numerical solutions.
2.2. Numerical Oscillations of the CN Method
Although the CN method is unconditionally stable and of second-order accuracy in both spatial and temporal directions, it may introduce spurious oscillations into the numerical solution for nonsmooth data. For simplicity, consider a homogeneous diffusion equation with discontinuous initial values defined in the one-dimensional (1D) space:of which the exact solution is given byfor .
Figure 1 depicts the exact and numerical solutions evolved by the CN and the variable- method [24], while Figure 2 compares the numerical solutions at obtained by the three numerical methods. The numerical solutions are obtained with and . As one can see, spurious oscillations occur along with the step discontinuities in the numerical solution of the CN method, as shown in Figure 1(b). On the contrary, as in Figures 1(c) and 2, the variable- method results in an accurate numerical solution without any oscillations. One should notice from Figures 1(a) and 2 that the implicit method is immune to spurious oscillations, but its error is considerably large due to a first-order discretization error in the temporal direction. Although the CN method reveals spurious oscillations, it is quite accurate far from discontinuities. The variable- method results in numerical solutions which are smooth as for the implicit method and accurate as for the CN method associated with smooth data.
(a)
(b)
(c)
The variable- method is a hybrid time-stepping procedure that is based on the CN method and alternately using the implicit method at points, where the numerical solution shows a certain portent of oscillations or reduced smoothness (the wobble set).
3. A Variable- Method for Two-Component Nonlinear Problems
This section introduces an effective time-stepping procedure for the numerical solution of two-component RD problems (1).
3.1. Linearization through Extrapolation
Once the spatial derivatives are approximated by second-order finite difference schemes, as in Section 2.1, the semidiscrete problem for (1) is formulated as
Let numerical solutions be obtained up to the th time level, . For the numerical solution in the th level, we first extrapolate numerical solutions in the two previous levels to approximate the solution at the midpoint :
See [12], for details of second-order extrapolation schemes for . Then, the -method for the two-component RD problem reads:where , , and .
The linearized problem (19) can be resolved by solving for two separate components . Each component in (19) can be formulated as follows:where , , and denote, respectively, , , and , for or 2, and is a known source term. The -method (20) can be rewritten in a vector form as
We present here the main steps of variable- method for a nonoscillatory solution of (20); a complete study of the method for diffusion equation was published in [24].
3.2. The Variable- Method
The method begins with defining the wobble set, the set of wobble points, as a collection of the grid points where the solution has high fluctuations so that the implicit method should be applied for the numerical solution not to develop oscillations.
One can easily verify that numerical oscillations of the CN method occur when its explicit half step produces spurious oscillations. Such nonphysical oscillations may happen particularly when the time step size is larger than the stability limit of the explicit scheme. Thus, the wobble set may be formed to include points where the explicit half step of the CN method introduces undesired local extrema. It follows from (21) that the explicit half step of the CN method () reads
Let be an interior grid point and consider the four partial directions (made with eight vicinal points of ): four directions having , , , and from the positive -direction. When spurious oscillations are observed in at least one direction, we select the point as a wobble point.
Define an index function for local extrema (idxt) as
Let , , and be point indices and define
Then, the wobble set (for the computation of ) is defined as
Remark 1. The function iswb selects candidates for the wobble set from local extrema satisfying ; however, the conditionexcludes cases where a strict extremum in becomes a strict extremum in the same sense in . Thus, the wobble set (25) is the set of interior grid points where becomes a local extremum while is either a nonextreme or an extreme in the opposite sense, for at least one of four partial directions.
Having the wobble set, the parameter for the computation of can be assigned pointwise:Thus, the variable-method for (20) can be formulated asor, in a vector form after grouping variables:The variable- method is analyzed for its numerical stability and accuracy and verified for various examples [24]. It results in nonoscillatory numerical solutions of which the accuracy is almost second-order in time.
Remark 2. The ADI procedure was also applied to (19) for which the initial values show sharp transitions. It has been observed that ADI may introduce undesirable/discontinuous features to its solution unless the time step size is sufficiently small, i.e., . The main problem with ADI is that the diffusion becomes anisotropic, i.e., faster in the coordinate directions. When ADI is applied to the variable- formulation of (19), the anisotropic features are reduced significantly. However, it requires to set small enough for a reliable numerical solution.
The algebraic system in (29) will be solved by applying the SOR method, with its initial value at the time level being set asIn particular, SOR converges quite fast for an appropriate choice of the relaxation parameter .
In the following, we will consider how to tune the optimal relaxation parameter for SOR.
4. The Optimal SOR Parameter
In this section, we will try to find a relaxation parameter which is heuristically optimal. Let us begin with the 2D algebraic system of (21) with :where , , and is the dimension of the algebraic system. It is known that the optimal relaxation parameter for the SOR method can be determined as ([25], Section 4.3)where is the spectral radius of the Jacobi iteration matrix .
For simplicity, assume that the problem is defined on the unit square with a Dirichlet boundary condition. We further assume that the domain is partitioned into grids so that . Then, the eigenvalues of the second-order 5-point FD coefficient matrix read ([25], Section 6.5)and therefore the eigenvalues of can be formulated asfor . Note that the diagonal element of is
So, the eigenvalues of the Jacobi iteration matrix are given as
In order to find the maximum of , we first obtainfor some . Here, we have used and the approximation . Now the spectral radius of reads
Assuming that , we finally obtainfor some .
It follows from (32) and (39) that the optimal SOR parameter , corresponding to the spatial grid size and the time step size , can be written asfor some . The constant can be found experimentally from a selected set of (Δt, h), as summarized in the following:
Once is estimated as in (41), the parameter in (40) is near-optimal for various choices of .
5. Numerical Experiments
In this section, we present numerical experiments which show effectiveness of the variable- method. The algorithms are implemented in MATLAB and carried out on a desktop computer of Intel Xeon CPU E5-1620 3.60 GHz processor.
To solve the algebraic system at each time level, the SOR method is employed with the near-optimal parameter calculated as in (40), with being estimated with a small grid problem. The SOR iteration is stopped when the maximum difference of consecutive iterates becomes smaller than a tolerance :
The -error , measured at , is defined as follows:where is the exact solution.
5.1. One-Component RD System
To investigate accuracy of the variable- method, we consider the diffusion problem (15) studied earlier in Section 2.2. For a comparison purpose, we have implemented not only the -methods and the variable- method but also the implicit predictor-corrector -Padé (IPC-[0,2]) method [26] and Lawson and Morris (LM) local extrapolation method [23].
Table 1 presents the -error at when the five methods are applied for the numerical solution of (15), with various and . As one can see from the table, the CN method resolves its numerical solution poorly (due to oscillations), except for the case the method satisfies the maximum principle. On the contrary, the variable- method results in a second-order accuracy, with its errors being smallest among all the methods for most cases. The new method is a hybrid time-stepping procedure which assembles merits from the CN method (high-accuracy) and the implicit method (smoothness).
Now, consider a nonlinear RD problem of the formwith the boundary and initial values, as given in (15).
Figure 3 presents the numerical solutions evolved by the implicit method, the CN method, and the variable- method, when and (the mesh is the same as the one selected in Figure 1). Similar to the linear problem in Figure 1, spurious oscillations are introduced into the numerical solution of the nonlinear problem by the CN method only. It should be noticed that spurious oscillations of the CN method are damped out much faster for the nonlinear problem than the linear problem, which is due to the reaction kinetic term . For the nonlinear problem, It seems that the oscillations at early time steps do not affect the solution at later steps much. This observation explains a partial reason that the second-order CN method has been popular for the numerical solution of PDEs in mathematical biology. However, for other applications, such spurious oscillations at early moments may alter the numerical solution significantly so as for the CN method to be unstable; see Figure 4 below. It is important to develop an effective algorithm which can suppress spurious oscillations for convenient choices of algorithm parameters; the variable- method is effective and stable.
(a)
(b)
(c)
(a)
(b)
(c)
(d)
(e)
(f)
5.2. Two-Component Nonlinear RD Systems
Two-component RD systems enable to explain a much wider range of phenomena than their one-component counterparts. Many two-component models have been developed and numerically verified for dynamical patterning behaviors in biology and chemistry. In this section, we consider two-component models interested in the literature of biology and chemistry, to verify effectiveness of the variable- method.
5.2.1. Gray–Scott Model in 1D
We apply the numerical methods for the numerical solution of the Gray–Scott model [4, 27] defined as (1) associated with the following reaction kinetics:for any constants and . Let . We assign two sets of model constants and initial and boundary conditions as follows [28]:where (46b) is a mid-pulse initial condition, andwhere (47b) is a left-pulse initial condition.
In Figure 5, we present the propagation of the numerical solution of by the variable- method associated with (46a)–(46c) at the grid sizes and over . The initial midpulse splits in early moments to travel in both directions (self-replication of the pulse). As each of the pulses travels, it becomes thicker (bigger) up to a certain width and begins to replicate itself recursively.
(a)
(b)
Figure 6 depicts the propagation of the numerical solution of by the variable- method associated with (47a)–(47c) over , at the same resolution, as in Figure 5. One can clearly observe a traveling pulse which begins from the left edge point and reflects whenever it hits the boundary, due to the no-flux boundary condition (47c).
(a)
(b)
To investigate bifurcation in the RD system and our method numerical accuracy, we present numerical solutions of the wave-splitting problem (46a)–(46c) obtained with various spatial and temporal grid sizes, as shown in Figure 7. The image represents the numerical solution obtained with the mesh resolution . For example, the image is associated with the mesh resolution . One can easily point out from the images that the spatial resolution alters the numerical solution dramatically even with halved spatial grid sizes (compare the images horizontally), while the temporal resolution affects little the numerical solution even with one-order smaller temporal step sizes (compare them vertically).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
The main reason for such a sensitivity to the spatial resolution is that the RD does not have as much time before growing to reach the margins of the mesh in low spatial resolutions (of large ’s). When this happens, the RD pattern typically deteriorates and it does not travel in an appropriate speed nor reaches a condition to replicate itself on time, see ([29], Section 4.2) for similar observations. We summarize the experiments with the Gray–Scott model in 1D as follows:(i)Accuracy of the numerical solution is much more sensitive to the spatial mesh resolution than the temporal one.(ii)Thus, it is crucial to set a high spatial resolution (small ’s) for a desirable accuracy.(iii)As the spatial resolution becomes higher, the CN method may more likely produce spurious oscillations, while the variable- method results in stable solutions.
5.2.2. Gray–Scott Model in 2D
Note that the two-component Gray–Scott model is formulated as in (1) with the reaction kinetics given in (45). We choose problem coefficients as follows [15]:
For the purpose of error analysis, we select a smooth solution defined asand replace the reaction kinetics with :
Then, in (49) would be the exact solution of with the initial condition .
Table 2 summarizes the -error with and the elapsed time (CPU) for the implicit, CN, variable- methods for three different meshes refined by a factor of 2 in both spatial and temporal directions with . Since the solution including the initial condition is smooth over the entire time interval , the CN method introduces no spurious oscillations into its numerical solution and proves a second-order accuracy for the two-component Gray–Scott model in the 2D space. One should notice that the variable- has also proved its second-order accuracy, the same as the CN method. On the contrary, the implicit method involves considerable errors due to its first-order convergence in temporal direction. Figure 8 shows the numerical solution by the variable- method, and its error at for the Gray–Scott model, when and .
(a)
(b)
(c)
(d)
For all the three methods, the algebraic system is solved by the SOR method with its optimal relaxation parameter being calibrated from the lowest resolution, . That is, the constant in (41) is evaluated using the experimentally optimal with and then (40) is utilized to estimate for other grid sizes . With the near-optimal parameter and an effective initialization scheme in (30), for both the CN and variable- methods, the SOR method has converged in (59) iterations in average for solving the two algebraic systems (for and ) in a time level. SOR is comparable with ADI in efficiency when the parameter is set optimal and the initialization is carried out accurately; SOR has proven its efficiency for the numerical solution of elliptic obstacle problems [30]. For the Gray–Scott model in 2D, the variable- method becomes about a third most expensive computationally than the CN method, due to the wobble set processing.
5.2.3. Gierer–Meinhardt Model
The Gierer–Meinhardt model [2] is (1) defined in with the following reaction kinetics and parameters:for which various numerical methods have been developed [12, 31, 32]. We cast the experiment employing coefficients and the initial condition used in [32]:
The initial values are depicted in Figure 9. In this section, we restrict our attention to the dynamics of of the model.
(a)
(b)
In order to investigate effectiveness of the variable- method and oscillatory behaviors of the CN method as well, we have carried out numerical experiments for the Gierer–Meinhardt model with a relatively low spatial resolution. Figure 4 presents numerical solutions at two different times for of the Gierer–Meinhardt model with the spatial resolution . When the time step size is set , the variable- method evolves the numerical solution as shown in Figures 4(a) and 4(d), for which the final steady-state pattern is the same as that in [32]. On the contrary, with the same , the CN method has produced a quite different pattern as in Figures 4(b) and 4(e), due to the nonsmooth initial values (Figure 9). However, the CN method can recover the correct steady-state pattern when it runs with , as depicted in Figures 4(c) and 4(f). For a similar accuracy, the variable- method (taking 170 s) is about 7 times more efficient than the CN method (taking 1242 s).
We summarize our experiments with the Gray–Scott and Gierer–Meinhardt models in 2D as follows:(i)The variable- method shows the same accuracy as the CN method for problems of smooth data(ii)For nonsmooth data, the variable- method evolves a smooth solution for all choices of , while the CN method introduces spurious oscillations to alter the solution unless the time step size is sufficiently small(iii)When a large time step size is desirable, the suggested method is a few times more efficient than the CN method for a similar accuracy
Remark 3. Although the variable- method can employ larger time step sizes than the CN method to get stable numerical solutions for problems of nonsmooth data, one may not set the time step size too large, due to an accuracy issue rather than the stability issue. Furthermore, for nonlinear problems, the overall stability of the numerical algorithm can be determined by not only grid sizes but also numerical schemes including methods of dealing with the nonlinear terms.
Figure 10 presents numerical solutions and their errors of for the Gray–Scott model (48)–(50) at by the variable- method, varying with fixed . Compared with Figure 8, the solutions show stability and a good accuracy, although the time step size is as large as . As shown in the bottom line in Figure 10, all three spatially different cases show the same level of errors since the entire errors are dominated by temporal direction errors. We conclude from this example that grid sizes in both temporal and spatial directions would not significantly affect the stability of the proposed method when the initial condition is smooth and the nonlinearity is not severe.
As an example of nonsmooth data and severe nonlinearity, we select the Gierer–Meinhardt model (51) and (52) to simulate with large temporal step sizes. When , the proposed algorithm introduced a rapid decay of solution values independently of the spatial grid size, so that the pattern is not formed appropriately. We believe that it is due to the error incorporated with the reaction term (19) when is approximated by the extrapolation scheme (18). However, when , our method produces stable solutions for all choices of spatial grid sizes. Figure 11 presents numerical solutions of for the Gierer–Meinhardt model at by the variable- method with fixed and various . Note that for Gierer–Meinhardt model, the pattern forming is slow down as the spatial grid size becomes smaller, as shown in and of Figure 11; this tendency has been observed for all other choices of . This is another example that accuracy of the numerical solution is much more sensitive to the spatial mesh resolution than the temporal one.
(a)
(b)
(c)
(d)
(e)
(f)
(a)
(b)
(c)
(d)
(e)
(f)
6. Conclusions
The Crank–Nicolson (CN) method has been a popular second-order time-stepping procedure for the numerical solution of systems of nonlinear RD equations. However, the CN method may introduce spurious oscillations for nonsmooth data unless the time step size is sufficiently small. We have studied a nonoscillatory time-stepping procedure for RD equations, called a variable- method, as a perturbation of the CN method. In each time level, the new method detects points of potential oscillations and resolves the solution applying the implicit method locally at those points. The proposed time-stepping procedure has proven nonoscillatory and having a second-order temporal accuracy, although the initial conditions are nonsmooth. Various examples have been considered to show effectiveness of the method. We also have performed a sensitivity analysis for the numerical solution of biological pattern forming models to conclude that the numerical solution is much more sensitive to the spatial mesh resolution than the temporal one.
Data Availability
The experiment data and figures used to support the findings of this study are included within the article, and the MATLAB codes for the experiments are available from the corresponding author upon request.
Conflicts of Interest
The authors declare no conflicts of interest.
Acknowledgments
This research was supported by NSF-MCB (1714157) awarded to George V. Popescu.