Research Article  Open Access
Philku Lee, George V. Popescu, Seongjai Kim, "A Nonoscillatory SecondOrder TimeStepping Procedure for ReactionDiffusion Equations", Complexity, vol. 2020, Article ID 5163704, 15 pages, 2020. https://doi.org/10.1155/2020/5163704
A Nonoscillatory SecondOrder TimeStepping Procedure for ReactionDiffusion Equations
Abstract
After a theory of morphogenesis in chemical cells was introduced in the 1950s, much attention had been devoted to the numerical solution of reactiondiffusion (RD) partial differential equations (PDEs). The Crank–Nicolson (CN) method has been a common secondorder timestepping 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 secondorder timestepping procedure for RD equations, called a variablemethod, 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 timestepping procedure is nonoscillatory and of a secondorder 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 reactiondiffusion (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], chloriteiodidemalonic 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 secondorder timestepping procedure. Timestepping 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 timestepping 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 timedependent multidimensional problems, either parabolic or hyperbolic [19, 20]. ADI reduces a multidimensional problem to multiple easytosolve onedimensional 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 variablemethod in which the timestepping 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 twocomponent 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 timestepping 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 twocomponent nonlinear RD equations. We adopt the successive overrelaxation (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 timestepping 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 spacetime 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 secondorder 5point 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 semiimplicit 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 timestepping procedure for the numerical solution of parabolic problems because it is stable and of secondorder 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 firstorder 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 secondorder 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 onedimensional (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 firstorder 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 timestepping 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 TwoComponent Nonlinear Problems
This section introduces an effective timestepping procedure for the numerical solution of twocomponent RD problems (1).
3.1. Linearization through Extrapolation
Once the spatial derivatives are approximated by secondorder 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 secondorder extrapolation schemes for . Then, the method for the twocomponent 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 variablemethod 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 secondorder 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 secondorder 5point 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 nearoptimal 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 E51620 3.60 GHz processor.
To solve the algebraic system at each time level, the SOR method is employed with the nearoptimal 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. OneComponent 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 predictorcorrector 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 secondorder accuracy, with its errors being smallest among all the methods for most cases. The new method is a hybrid timestepping procedure which assembles merits from the CN method (highaccuracy) 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 secondorder 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. TwoComponent Nonlinear RD Systems
Twocomponent RD systems enable to explain a much wider range of phenomena than their onecomponent counterparts. Many twocomponent models have been developed and numerically verified for dynamical patterning behaviors in biology and chemistry. In this section, we consider twocomponent 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 midpulse initial condition, andwhere (47b) is a leftpulse 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 (selfreplication 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 noflux boundary condition (47c).
(a)
(b)
To investigate bifurcation in the RD system and our method numerical accuracy, we present numerical solutions of the wavesplitting 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 oneorder 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 twocomponent 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 secondorder accuracy for the twocomponent Gray–Scott model in the 2D space. One should notice that the variable has also proved its secondorder accuracy, the same as the CN method. On the contrary, the implicit method involves considerable errors due to its firstorder 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 nearoptimal 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 steadystate 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 steadystate 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 secondorder timestepping 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 timestepping 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 timestepping procedure has proven nonoscillatory and having a secondorder 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 NSFMCB (1714157) awarded to George V. Popescu.
References
 A. M. Turing, “The chemical basis of morphogenesis,” Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences, vol. 237, no. 641, pp. 37–72, 1952. View at: Google Scholar
 A. Gierer and H. Meinhardt, “A theory of biological pattern formation,” Kybernetik, vol. 12, no. 1, pp. 30–39, 1972. View at: Publisher Site  Google Scholar
 J. Schnakenberg, “Simple chemical reaction systems with limit cycle behaviour,” Journal of Theoretical Biology, vol. 81, no. 3, pp. 389–400, 1979. View at: Publisher Site  Google Scholar
 P. Gray and S. K. Scott, “Autocatalytic reactions in the isothermal, continuous stirred tank reactor: isolas and other forms of multistability,” Chemical Engineering Science, vol. 38, no. 1, pp. 29–43, 1983. View at: Publisher Site  Google Scholar
 I. Prigogine and G. Nicolis, “Selforganisation in nonequilibrium systems: towards a dynamics of complexity,” in Bifurcation Analysis, pp. 3–12, Springer, Berlin, Germany, 1985. View at: Google Scholar
 I. Lengyel and I. R. Epstein, “Modeling of turing structures in the chlorite–iodide–malonic acid–starch reaction system,” Science, vol. 251, no. 4994, pp. 650–652, 1991. View at: Publisher Site  Google Scholar
 S. Kondo and T. Miura, “Reactiondiffusion model as a framework for understanding biological pattern formation,” Science, vol. 329, no. 5999, pp. 1616–1620, 2010. View at: Publisher Site  Google Scholar
 P. K. Maini and R. E. Baker, Developmental Biology: Mathematical Modelling of Development, Wiley, Hoboken, NJ, USA, 2001.
 K. U. Torii, “Twodimensional spatial patterning in developmental systems,” Trends in Cell Biology, vol. 22, no. 8, pp. 438–446, 2012. View at: Publisher Site  Google Scholar
 E. Crampin and P. Maini, “Reactiondiffusion models for biological pattern formation,” Methods and Applications of Analysis, vol. 8, no. 3, pp. 415–428, 2001. View at: Publisher Site  Google Scholar
 C. Chiu and N. Walkington, “An ADI method for hysteretic reactiondiffusion systems,” SIAM Journal on Numerical Analysis, vol. 34, no. 3, pp. 1185–1206, 1997. View at: Publisher Site  Google Scholar
 R. I. Fernandes, B. Bialecki, and G. Fairweather, “An ADI extrapolated CrankNicolson orthogonal spline collocation method for nonlinear reactiondiffusion systems on evolving domains,” Journal of Computational Physics, vol. 299, pp. 561–580, 2015. View at: Publisher Site  Google Scholar
 I. Sgura, B. Bozzini, and D. Lacitignola, “Numerical approximation of Turing patterns in electrodeposition by ADI methods,” Journal of Computational and Applied Mathematics, vol. 236, no. 16, pp. 4132–4147, 2012. View at: Publisher Site  Google Scholar
 F. Shakeri and M. Dehghan, “The finite volume spectral element method to solve Turing models in the biological pattern formation,” Computers & Mathematics with Applications, vol. 62, no. 12, pp. 4322–4336, 2011. View at: Publisher Site  Google Scholar
 R. I. Fernandes and G. Fairweather, “An ADI extrapolated CrankNicolson orthogonal spline collocation method for nonlinear reactiondiffusion systems,” Journal of Computational Physics, vol. 231, no. 19, pp. 6248–6267, 2012. View at: Publisher Site  Google Scholar
 J. Douglas Jr., “On the numerical integration of