# A Nonoscillatory Second-Order Time-Stepping Procedure for Reaction-Diffusion Equations

Authors: Philku Lee; George V. Popescu; Seongjai Kim
Journal: Complexity
 (2020)

---

## Abstract

After a theory of morphogenesis in chemical cells was introduced in the 1950s, much attention had been devoted to the numerical solution ofreaction-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.


---

## Body

## 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, thereaction-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 ℝd, d=1,2,3, Γ=∂Ω denote the boundary of Ω, and J=0,T for some T>0:(1)∂u∂t−DΔu=fu,Ω×J,∂u∂ν=0,Γ×J,u=u0,Ω×t=0,where u=u1,u2T, D=diagD1,D2T is the diffusion tensor, Δ denotes the Laplace operator, ∂/∂ν is the outward normal derivative on the boundary Γ, and fu is the reaction kinetics of the system given as(2)fu=f1u1,u2f2u1,u2.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 OΔt2, where Δt 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 principle12, 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., Δt=OΔx2, where Δx 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, θ∈0,1, 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. Section2 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 oflinear parabolic equations of the form(3)∂u∂t−Δu=f,Ω×J,∂u∂ν=0,Γ×J,u=u0,Ω×t=0,where D>0 is a diffusion coefficient and f 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 ℝ2: Ω=ax,bx×ay,by. By partitioning Ω×J, we obtain the space-time grid points:(4)xij,tn≔xi,yj,tn,i=0,1,…,nx,j=0,1,…,ny,n=0,1,…,nt,where nx, ny, and nt are prescribed positive integers and(5)xi=ax+i⋅Δx,yj=ay+j⋅Δy,tn=n⋅Δt,Δx=bx−axnx,Δy=by−ayny,t=Tnt.Define the discrete domain, the set of the spatial grid points, by(6)Ωd=xi,yj:0≤i≤nx,0≤j≤ny,and denote the set of boundary grid points by Γd=Ωd∩Γ and the set of interior grid points by Ωd0=Ωd\Γd.Letgijn=gxij,tn for all functions g defined in x,t. Then, the second-order 5-point finite difference (FD) approximation A of −Δ at t=tn reads(7)Auijn=A1uijn+A2uijn,A1uijn≔−δxxuijn,A2uijn≔−δyyuijn,where the FD operators are defined as(8)δxxuijn=ui−1,jn−2uijn+ui+1,jnΔx2,δyyuijn=ui,j−1n−2uijn+ui,j+1nΔy2.For the temporal derivative∂u/∂t, a convenient FD approximation can give(9)uijn+1−uijnΔt.Expressing the spatial derivative by a weighted averageθ∈0,1 of previous and current time values, we can formulate the θ-method for (3) as(10)uijn+1−uijnΔt+Aθuijn+1+1−θuijn=fijn+θ,where fijn+θ is either fxij,tn+θ or θfijn+1+1−θfijn. A simple algebraic rearrangement of (10) in a vector form becomes(11)I+θΔtAun+1=I−1−θΔtAun+Δtfn+θ,where un=uijn0≤i≤nx,0≤j≤ny and fn+θ=fijn+θ0≤i≤nx,0≤j≤ny, considered as column vectors. Popular choices of θ∈0,1 are 0, 1, and 1/2, 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θ=0, algorithm (11) is stable when it satisfies(12)μx+μy≤12,μx=ΔtΔx2,μy=ΔtΔy2.Although the explicit method is efficient for each time step, its stability condition in (12) enforces, choosing a small time step size Δt; it may become less efficient compared with other implicit methods. It is elementary in numerical analysis that when θ≥1/2, the θ-method (10) is unconditionally stable.(ii)
Crank–Nicolson method: whenθ=1/2, (11) can be rewritten as(13)I+Δt2Aun+1=I−Δt2Aun+Δtfn+1/2.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 θij=1 at points xij 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θ=1, algorithm (11) reads(14)I+ΔtAun+1=un+Δtfn+1.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:(15)∂u∂t−uxx=0,x,t∈0,1×0,T,u0,t=u1,t=0,t∈0,T,ux,0=0if0<x<14,1if14≤x<34,0if34≤x<1,of which the exact solution is given by(16)ux,t=∑k=1∞4kπsinkπ2sinkπ4sinkπxe−k2π2t,for x,t∈0,1×0,T.Figure1 depicts the exact and numerical solutions evolved by the CN and the variable-θ method [24], while Figure 2 compares the numerical solutions at T=1.0 obtained by the three numerical methods. The numerical solutions are obtained with Δt=0.01 and Δx=Δy=0.025. 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.Figure 1
Propagation of the exact and numerical solutions for (15): (a) the exact solution and the numerical solution by (b) the CN method and (c) the variable-θ method [24], for 0≤t≤T=1.0, when Δt=0.01 and Δx=0.025.
(a)
(b)
(c)Figure 2
The numerical solutions of (15) at T=1.0, compared with the exact solution.The variable-θ method is a hybrid time-stepping procedure that is based on the CN method θ=1/2 and alternately using the implicit method θ=1 at points, where the numerical solution shows a certain portent of oscillations or reduced smoothness (the wobble set).


## 2.1. Theθ-Method: Difference Equation

LetΩ be a rectangular domain in ℝ2: Ω=ax,bx×ay,by. By partitioning Ω×J, we obtain the space-time grid points:(4)xij,tn≔xi,yj,tn,i=0,1,…,nx,j=0,1,…,ny,n=0,1,…,nt,where nx, ny, and nt are prescribed positive integers and(5)xi=ax+i⋅Δx,yj=ay+j⋅Δy,tn=n⋅Δt,Δx=bx−axnx,Δy=by−ayny,t=Tnt.Define the discrete domain, the set of the spatial grid points, by(6)Ωd=xi,yj:0≤i≤nx,0≤j≤ny,and denote the set of boundary grid points by Γd=Ωd∩Γ and the set of interior grid points by Ωd0=Ωd\Γd.Letgijn=gxij,tn for all functions g defined in x,t. Then, the second-order 5-point finite difference (FD) approximation A of −Δ at t=tn reads(7)Auijn=A1uijn+A2uijn,A1uijn≔−δxxuijn,A2uijn≔−δyyuijn,where the FD operators are defined as(8)δxxuijn=ui−1,jn−2uijn+ui+1,jnΔx2,δyyuijn=ui,j−1n−2uijn+ui,j+1nΔy2.For the temporal derivative∂u/∂t, a convenient FD approximation can give(9)uijn+1−uijnΔt.Expressing the spatial derivative by a weighted averageθ∈0,1 of previous and current time values, we can formulate the θ-method for (3) as(10)uijn+1−uijnΔt+Aθuijn+1+1−θuijn=fijn+θ,where fijn+θ is either fxij,tn+θ or θfijn+1+1−θfijn. A simple algebraic rearrangement of (10) in a vector form becomes(11)I+θΔtAun+1=I−1−θΔtAun+Δtfn+θ,where un=uijn0≤i≤nx,0≤j≤ny and fn+θ=fijn+θ0≤i≤nx,0≤j≤ny, considered as column vectors. Popular choices of θ∈0,1 are 0, 1, and 1/2, 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θ=0, algorithm (11) is stable when it satisfies(12)μx+μy≤12,μx=ΔtΔx2,μy=ΔtΔy2.Although the explicit method is efficient for each time step, its stability condition in (12) enforces, choosing a small time step size Δt; it may become less efficient compared with other implicit methods. It is elementary in numerical analysis that when θ≥1/2, the θ-method (10) is unconditionally stable.(ii)
Crank–Nicolson method: whenθ=1/2, (11) can be rewritten as(13)I+Δt2Aun+1=I−Δt2Aun+Δtfn+1/2.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 θij=1 at points xij 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θ=1, algorithm (11) reads(14)I+ΔtAun+1=un+Δtfn+1.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:(15)∂u∂t−uxx=0,x,t∈0,1×0,T,u0,t=u1,t=0,t∈0,T,ux,0=0if0<x<14,1if14≤x<34,0if34≤x<1,of which the exact solution is given by(16)ux,t=∑k=1∞4kπsinkπ2sinkπ4sinkπxe−k2π2t,for x,t∈0,1×0,T.Figure1 depicts the exact and numerical solutions evolved by the CN and the variable-θ method [24], while Figure 2 compares the numerical solutions at T=1.0 obtained by the three numerical methods. The numerical solutions are obtained with Δt=0.01 and Δx=Δy=0.025. 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.Figure 1
Propagation of the exact and numerical solutions for (15): (a) the exact solution and the numerical solution by (b) the CN method and (c) the variable-θ method [24], for 0≤t≤T=1.0, when Δt=0.01 and Δx=0.025.
(a)
(b)
(c)Figure 2
The numerical solutions of (15) at T=1.0, compared with the exact solution.The variable-θ method is a hybrid time-stepping procedure that is based on the CN method θ=1/2 and alternately using the implicit method θ=1 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 Section2.1, the semidiscrete problem for (1) is formulated as(17)∂u∂t+DAu=fu,t∈0,T,u=u0,t=0.Let numerical solutions be obtained up to thenth time level, n>0. For the numerical solution in the n+1th level, we first extrapolate numerical solutions in the two previous levels to approximate the solution at the midpoint tn+1/2:(18)u˜n+1/2≔32un−12un−1.See [12], for details of second-order extrapolation schemes for n≥0. Then, the θ-method for the two-component RD problem reads:(19)uijn+1−uijnΔt+DAθuijn+1+I−θuijn=fu˜n+1/2ij,where u=u1,u2T, D=diagD1,D2T, and θ=diagθ1,θ2T.The linearized problem (19) can be resolved by solving for two separate components un+1=u1n+1,u2n+1T. Each component in (19) can be formulated as follows:(20)uijn+1−uijnΔt+DAθuijn+1+I−θuijn=fijn+1/2,where u, D, and θ denote, respectively, uk, Dk, and θk, for k=1 or 2, and fn+1/2 is a known source term. The θ-method (20) can be rewritten in a vector form as(21)I+θΔtDAun+1=I−1−θΔtDAun+Δtfn+1/2.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 thewobble set, the set of wobble points, as a collection of the grid points where the solution has high fluctuations so that the implicit method θ=1 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Δt 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 (θ=1/2) reads(22)un,∗≔I−Δt2DAun.Letxij be an interior grid point and consider the four partial directions (made with eight vicinal points of xij): four directions having 0°, 45°, 90°, and 135° from the positive x-direction. When spurious oscillations are observed in at least one direction, we select the point xij as a wobble point.Define an index function for local extrema (idxt) as(23)idxta,b,c=0,ifmina,c<b<maxa,c,1,ifb=maxa,c,−1,ifb=mina,c,2,ifmaxa,c<b,−2,ifb<mina,c.LetP, Q, and R be point indices and define(24)iswbP,Q,R,n=1,if idxtuPn,∗,uQn,∗,uRn,∗≠0andidxtuPn,∗,uQn,∗,uRn,∗+idxtuPn,uQn,uRn<4,0,otherwise.Then, the wobble set (for the computation ofun+1=uijn+1) is defined as(25)Wn=xij∈Ωd0∣iswbi−1,j,i,j,i+1,j,n=1or iswbi−1,j−1,i,j,i+1,j+1,n=1or iswbi,j−1,i,j,i,j+1,n=1or iswbi+1,j−1,i,j,i−1,j+1,n=1.Remark 1.
The functioniswb selects candidates for the wobble set from local extrema satisfying idxtuPn,∗,uQn,∗,uRn,∗≠0; however, the condition(26)idxtuPn,∗,uQn,∗,uRn,∗+idxtuPn,uQn,uRn<4,excludes cases where a strict extremum in un becomes a strict extremum in the same sense in un,∗. Thus, the wobble set (25) is the set of interior grid points xij where uijn,∗ becomes a local extremum while uijn 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 un+1 can be assigned pointwise:(27)θijn+1≔θxij,tn+1=1,ifxij∈Wn,1/2,otherwise.
Thus, thevariable-θmethod for (20) can be formulated as(28)uijn+1−uijnΔt+DAθijn+1uijn+1+1−θijn+1uijn=fijn+1/2,or, in a vector form after grouping variables:(29)I+θn+1ΔtDAun+1=I−1−θn+1ΔtDAun+Δtfn+1/2.
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., Δt=OΔx2,Δy2. 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 Δt 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 tn+1 being set as(30)un+1,0=2un−un−1.
In 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.


## 3.1. Linearization through Extrapolation

Once the spatial derivatives are approximated by second-order finite difference schemes, as in Section2.1, the semidiscrete problem for (1) is formulated as(17)∂u∂t+DAu=fu,t∈0,T,u=u0,t=0.Let numerical solutions be obtained up to thenth time level, n>0. For the numerical solution in the n+1th level, we first extrapolate numerical solutions in the two previous levels to approximate the solution at the midpoint tn+1/2:(18)u˜n+1/2≔32un−12un−1.See [12], for details of second-order extrapolation schemes for n≥0. Then, the θ-method for the two-component RD problem reads:(19)uijn+1−uijnΔt+DAθuijn+1+I−θuijn=fu˜n+1/2ij,where u=u1,u2T, D=diagD1,D2T, and θ=diagθ1,θ2T.The linearized problem (19) can be resolved by solving for two separate components un+1=u1n+1,u2n+1T. Each component in (19) can be formulated as follows:(20)uijn+1−uijnΔt+DAθuijn+1+I−θuijn=fijn+1/2,where u, D, and θ denote, respectively, uk, Dk, and θk, for k=1 or 2, and fn+1/2 is a known source term. The θ-method (20) can be rewritten in a vector form as(21)I+θΔtDAun+1=I−1−θΔtDAun+Δtfn+1/2.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 thewobble set, the set of wobble points, as a collection of the grid points where the solution has high fluctuations so that the implicit method θ=1 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Δt 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 (θ=1/2) reads(22)un,∗≔I−Δt2DAun.Letxij be an interior grid point and consider the four partial directions (made with eight vicinal points of xij): four directions having 0°, 45°, 90°, and 135° from the positive x-direction. When spurious oscillations are observed in at least one direction, we select the point xij as a wobble point.Define an index function for local extrema (idxt) as(23)idxta,b,c=0,ifmina,c<b<maxa,c,1,ifb=maxa,c,−1,ifb=mina,c,2,ifmaxa,c<b,−2,ifb<mina,c.LetP, Q, and R be point indices and define(24)iswbP,Q,R,n=1,if idxtuPn,∗,uQn,∗,uRn,∗≠0andidxtuPn,∗,uQn,∗,uRn,∗+idxtuPn,uQn,uRn<4,0,otherwise.Then, the wobble set (for the computation ofun+1=uijn+1) is defined as(25)Wn=xij∈Ωd0∣iswbi−1,j,i,j,i+1,j,n=1or iswbi−1,j−1,i,j,i+1,j+1,n=1or iswbi,j−1,i,j,i,j+1,n=1or iswbi+1,j−1,i,j,i−1,j+1,n=1.Remark 1.
The functioniswb selects candidates for the wobble set from local extrema satisfying idxtuPn,∗,uQn,∗,uRn,∗≠0; however, the condition(26)idxtuPn,∗,uQn,∗,uRn,∗+idxtuPn,uQn,uRn<4,excludes cases where a strict extremum in un becomes a strict extremum in the same sense in un,∗. Thus, the wobble set (25) is the set of interior grid points xij where uijn,∗ becomes a local extremum while uijn 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 un+1 can be assigned pointwise:(27)θijn+1≔θxij,tn+1=1,ifxij∈Wn,1/2,otherwise.
Thus, thevariable-θmethod for (20) can be formulated as(28)uijn+1−uijnΔt+DAθijn+1uijn+1+1−θijn+1uijn=fijn+1/2,or, in a vector form after grouping variables:(29)I+θn+1ΔtDAun+1=I−1−θn+1ΔtDAun+Δtfn+1/2.
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., Δt=OΔx2,Δy2. 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 Δt 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 tn+1 being set as(30)un+1,0=2un−un−1.
In 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 isheuristically optimal. Let us begin with the 2D algebraic system of (21) with θ=1/2:(31)ℒu=r∈ℝm×m,where ℒ=I+Δt/2DA, r=I−Δt/2DAu+Δtfn+1/2, and m>0 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)(32)ω^=21+1−ρTJ2,where ρTJ is the spectral radius of the Jacobi iteration matrix TJ.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 intoN×N grids so that h=Δx=Δy=1/N. Then, the eigenvalues of the second-order 5-point FD coefficient matrix A read ([25], Section 6.5)(33)λk,ℓA=1h24−2coskπN−2cosℓπN,1≤k,ℓ≤N−1,and therefore the eigenvalues of ℒ can be formulated as(34)λk,ℓℒ=1+Δt2Dλk,ℓA=1+DΔth22−coskπN−cosℓπN,for 1≤k,ℓ≤N−1. Note that the diagonal element of ℒ is(35)p≔1+Δt2D4h2=1+2DΔth2.So, the eigenvalues of the Jacobi iteration matrixTJ are given as(36)λk,ℓTJ=p−λk,ℓℒp=DΔth2coskπ/N+cosℓπ/N1+2DΔt/h2.In order to find the maximum ofλk,ℓTJ, we first obtain(37)maxkcoskπN=cosπN=1−c1h2,for some c1>0. Here, we have used h=1/N and the approximation cosx≈1−x2/2. Now the spectral radius of TJ reads(38)ρTJ=maxk,ℓλk,ℓTJ=1−c1h21+12h2DΔt−1.Assuming thatc1h2<h2/2DΔt<1, we finally obtain(39)ρTJ=1−c2h2DΔt,for some c2>0.It follows from (32) and (39) that the optimal SOR parameter ω^Δt,h, corresponding to the spatial grid size h and the time step size Δt, can be written as(40)ω^Δt,h=21+c0h/DΔt,for some c0>0. The constant c0 can be found experimentally from a selected set of (Δt, h), as summarized in the following:(41)aDetermineω^Δt0,h0for prescribed grid sizesΔt0,h0,heuristically.bSolve40forc0:c0=DΔt0h02ω^Δt0,h0−1.Oncec0 is estimated as in (41), the parameter ω^Δt,h in (40) is near-optimal for various choices of Δt,h.


## 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 c0 being estimated with a small grid problem. The SOR iteration is stopped when the maximum difference of consecutive iterates becomes smaller than a tolerance ε=10−6:(42)un,k−un,k−1∞<ε.TheL∞-error E∞tn, measured at t=tn, is defined as follows:(43)E∞tn≔un−u^tn∞,where u^ 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 0,2-Padé (IPC-[0,2]) method [26] and Lawson and Morris (LM) local extrapolation method [23].Table1 presents the L∞-error E∞T at T=1.0 when the five methods are applied for the numerical solution of (15), with various Δt and Δx. 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).Table 1
L
∞-error E∞T for (15) at T=1.0.
Δ
t
,
Δ
x
Implicit
CN
IPC-0,2
LM
Variable-θ
0.1
,
0.05
9.0
⋅
10
−
4
2.4
⋅
10
−
1
6.0
⋅
10
−
5
8.4
⋅
10
−
5
8.4
⋅
10
−
5
0.01
,
0.05
2.8
⋅
10
−
5
4.8
⋅
10
−
7
1.6
⋅
10
−
6
3.2
⋅
10
−
6
8.5
⋅
10
−
7
0.1
,
0.025
8.9
⋅
10
−
4
3.7
⋅
10
−
1
5.9
⋅
10
−
5
8.3
⋅
10
−
5
5.1
⋅
10
−
5
0.01
,
0.025
2.7
⋅
10
−
5
2.0
⋅
10
−
4
9.1
⋅
10
−
7
2.5
⋅
10
−
6
1.5
⋅
10
−
7
0.1
,
0.0125
8.9
⋅
10
−
4
4.3
⋅
10
−
1
5.8
⋅
10
−
5
8.3
⋅
10
−
5
2.7
⋅
10
−
4
0.01
,
0.0125
2.7
⋅
10
−
5
3.9
⋅
10
−
2
7.5
⋅
10
−
7
2.4
⋅
10
−
6
1.2
⋅
10
−
7Now, consider a nonlinear RD problem of the form(44)∂u∂t−uxx=u1−u,x,t∈0,1×0,T,with the boundary and initial values, as given in (15).Figure3 presents the numerical solutions evolved by the implicit method, the CN method, and the variable-θ method, when Δt=0.01 and Δx=0.025 (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 fu=u1−u. 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.Figure 3
Propagation of the numerical solutions for (44): (a) the implicit method, (b) the CN method, and (c) the variable-θ method, for 0≤t≤T=1.0, when Δt=0.01 and Δx=0.025.
(a)
(b)
(c)Figure 4
Numerical solutions foru1 of the Gierer–Meinhardt model at the spatial resolution Δx=Δy=1/32: (a) the variable-θ method with Δt=0.05, (b) the CN method with Δt=0.05, and (c) the CN method with Δt=0.005; (d–f) are the aerial views of (a–c), respectively.
(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:(45)fu=F1−u1−u1u22,u1u22−F+ku2T,for any constants F and k. Let Ω=0,1. We assign two sets of model constants and initial and boundary conditions as follows [28]:(46a)D=10−4,10−6T,F=0.035,k=0.049,(46b)u1x,0=1−12sin100πx,u2x,0=14sin100πx,x∈0,1,(46c)u10,t=u11,t=1,u20,t=u21,t=0,t∈0,T,where (46b) is a mid-pulse initial condition, and(47a)D=10−4,5⋅10−5T,F=0.025,k=0.0544,(47b)u1x,0=1−12cos100πx2,u2x,0=14cos100πx2,x∈0,1,(47c)∂u1∂x0,t=∂u1∂x1,t=0,∂u2∂x0,t=∂u2∂x1,t=0,t∈0,T,where (47b) is a left-pulse initial condition.In Figure5, we present the propagation of the numerical solution of u2 by the variable-θ method associated with (46a)–(46c) at the grid sizes Δt=0.01 and Δx=0.004 over 0≤t≤T=2000. 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.Figure 5
Propagation of the numerical solution ofu2 for the 1D Gray–Scott model by the variable-θ method: (a) the wave-splitting (self-replication of the pulse) and (b) its aerial view over 0≤t≤T=2000, when Δt=0.01 and Δx=0.004.
(a)
(b)Figure6 depicts the propagation of the numerical solution of u2 by the variable-θ method associated with (47a)–(47c) over 0≤t≤T=5000, 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).Figure 6
Propagation of the numerical solution ofu2 for the 1D Gray–Scott model by the variable-θ method: (a) the pulse traveling and reflecting and (b) its aerial view over 0≤t≤T=5000, when Δt=0.01 and Δx=0.004.
(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 Ikℓ represents the numerical solution obtained with the mesh resolution Δt,Δx=10−k,0.01/2ℓ−1. For example, the image I23 is associated with the mesh resolution Δt,Δx=1/100,1/400. 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).Figure 7
The wave-splitting (self-replication of the pulse) by the variable-θ method over 0≤t≤2000 with various Δt,Δx. The image Ikℓ represents the numerical solution obtained with the mesh resolution Δt,Δx=10−k,0.01/2ℓ−1. (a) (I11). (b) (I12). (c) (I13). (d) (I21). (e) (I22). (f) (I23). (g) (I31). (h) (I32). (i) (I33).
(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Δx’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Δx’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 fu given in (45). We choose problem coefficients as follows [15]:(48)Ω=−1,1×−1,1,D=0.001,0.001T,F=1,k=0.For the purpose of error analysis, we select a smooth solutionu^=u^1,u^2 defined as(49)u^1x,y,t=cos2tcos2πxcosπy,u^2x,y,t=cos2tcosπxcos2πy,and replace the reaction kinetics fu with fu^u:(50)fu^u≔fu+∂u^∂t−DΔu^−fu^.Then,u^=u^1,u^2 in (49) would be the exact solution of ∂u/∂t−DΔu−fu=fu^u with the initial condition u0=u^x,y,0.Table2 summarizes the L∞-error E∞T with T=1.0 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 Δx=Δy. Since the solution including the initial condition is smooth over the entire time interval 0,T=1, 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 T=1.0 for the Gray–Scott model, when Δt=0.025 and Δx=Δy=0.05.Table 2
L
∞-error E∞T and the elapsed time (CPU) for the numerical solution of the Gray–Scott model in the 2D space: T=1.0.
ε
=
10
−
6
Implicit
CN
Variable-θ
Δ
t
,
Δ
x
=
Δ
y
E
∞
T
CPU
E
∞
T
CPU
E
∞
T
CPU
u
1
u
2
u
1
u
2
u
1
u
2
0.05
,
0.1
4.3
⋅
10
−
2
4.0
⋅
10
−
2
0.036 s
9.6
⋅
10
−
4
9.4
⋅
10
−
4
0.047 s
9.6
⋅
10
−
4
9.4
⋅
10
−
4
0.062 s
0.025
,
0.05
2.0
⋅
10
−
2
1.7
⋅
10
−
2
0.225 s
2.5
⋅
10
−
4
2.3
⋅
10
−
4
0.280 s
2.5
⋅
10
−
4
2.3
⋅
10
−
4
0.394 s
0.0125
,
0.025
8.0
⋅
10
−
3
8.0
⋅
10
−
3
1.629 s
6.4
⋅
10
−
5
5.5
⋅
10
−
5
2.102 s
6.4
⋅
10
−
5
5.5
⋅
10
−
5
2.798 s
Conv. order
1.2
1.2
1.9
2.0
1.9
2.0Figure 8
The numerical solution of the variable-θ method and its error at T=1.0 for the Gray–Scott model, when Δt=0.025 and Δx=Δy=0.05: (a) u1, (b) u2, (c) e1=u1−u1^, and (d) e2=u2−u2^.
(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,Δt,Δx=Δy=0.05,0.1. That is, the constant c0 in (41) is evaluated using the experimentally optimal ω^Δt0,h0 with Δt0,h0=0.05,0.1 and then (40) is utilized to estimate ω^Δt,h for other grid sizes Δt,h=Δx=Δy. 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 u1 and u2) 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 Ω=−1,1×−1,1⊂ℝ2 with the following reaction kinetics and parameters:(51)D=ε2,κμT,fu=u12u2−u1,1μu12ε−u2T,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 Figure9. In this section, we restrict our attention to the dynamics of u1 of the model.(52)ε=0.04,μ=0.1,κ=0.0128,u1x,y,0=121+0.001∑k=120coskπy2sech2x2+y22ε,u2x,y,0=cosh1−x2+y23cosh1.Figure 9
The initial values in (52) for the Gierer–Meinhardt model with ϵ=0.04 at the mesh resolution Δx=Δy=1/32: (a) u1x,y,0 and (b) u2x,y,0.
(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 u1 of the Gierer–Meinhardt model with the spatial resolution Δx=Δy=1/32. When the time step size is set Δt=0.05, 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 Δt=0.05, 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 Δt=0.005, 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 Δt, 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 accuracyRemark 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.
Figure10 presents numerical solutions and their errors of u1 for the Gray–Scott model (48)–(50) at T=1.0 by the variable-θ method, varying Δx=Δy with fixed Δt=0.25. Compared with Figure 8, the solutions show stability and a good accuracy, although the time step size is as large as Δt=0.25. 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 Δt≥0.2, 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 un+1/2 is approximated by the extrapolation scheme (18). However, when Δt≤0.1, our method produces stable solutions for all choices of spatial grid sizes. Figure 11 presents numerical solutions of u1 for the Gierer–Meinhardt model at T=340 by the variable-θ method with fixed Δt=0.1 and various Δx=Δy. Note that for Gierer–Meinhardt model, the pattern forming is slow down as the spatial grid size becomes smaller, as shown in I2 and I3 of Figure 11; this tendency has been observed for all other choices of Δt≤0.1. This is another example that accuracy of the numerical solution is much more sensitive to the spatial mesh resolution than the temporal one.Figure 10
Numerical solutions foru1 of the Gray–Scott model at T=1.0 by the variable-θ method varying Δx=Δy with fixed Δt=0.25. The image Iℓ represents the numerical solution obtained with the mesh resolution Δt,Δx=Δy=0.25,0.02/2ℓ−1 and Iℓ′ represents the error of Iℓ. (a) (I1). (b) (I2). (c) (I3). (d)I1′. (e)I2′. (f)I3′.
(a)
(b)
(c)
(d)
(e)
(f)Figure 11
Numerical solutions foru1 of the Gierer–Meinhardt model at T=340 by the variable-θ method, varying Δx=Δy with fixed Δt=0.1. The image Iℓ represents the numerical solution obtained with the mesh resolution Δt,Δx=Δy=0.1,0.02/2ℓ−1 and Iℓ′ represents the aerial view of Iℓ. (a) (I1). (b) (I2). (c) (I3). (d)I1′. (e)I2′. (f)I3′.
(a)
(b)
(c)
(d)
(e)
(f)


## 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 0,2-Padé (IPC-[0,2]) method [26] and Lawson and Morris (LM) local extrapolation method [23].Table1 presents the L∞-error E∞T at T=1.0 when the five methods are applied for the numerical solution of (15), with various Δt and Δx. 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).Table 1
L
∞-error E∞T for (15) at T=1.0.
Δ
t
,
Δ
x
Implicit
CN
IPC-0,2
LM
Variable-θ
0.1
,
0.05
9.0
⋅
10
−
4
2.4
⋅
10
−
1
6.0
⋅
10
−
5
8.4
⋅
10
−
5
8.4
⋅
10
−
5
0.01
,
0.05
2.8
⋅
10
−
5
4.8
⋅
10
−
7
1.6
⋅
10
−
6
3.2
⋅
10
−
6
8.5
⋅
10
−
7
0.1
,
0.025
8.9
⋅
10
−
4
3.7
⋅
10
−
1
5.9
⋅
10
−
5
8.3
⋅
10
−
5
5.1
⋅
10
−
5
0.01
,
0.025
2.7
⋅
10
−
5
2.0
⋅
10
−
4
9.1
⋅
10
−
7
2.5
⋅
10
−
6
1.5
⋅
10
−
7
0.1
,
0.0125
8.9
⋅
10
−
4
4.3
⋅
10
−
1
5.8
⋅
10
−
5
8.3
⋅
10
−
5
2.7
⋅
10
−
4
0.01
,
0.0125
2.7
⋅
10
−
5
3.9
⋅
10
−
2
7.5
⋅
10
−
7
2.4
⋅
10
−
6
1.2
⋅
10
−
7Now, consider a nonlinear RD problem of the form(44)∂u∂t−uxx=u1−u,x,t∈0,1×0,T,with the boundary and initial values, as given in (15).Figure3 presents the numerical solutions evolved by the implicit method, the CN method, and the variable-θ method, when Δt=0.01 and Δx=0.025 (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 fu=u1−u. 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.Figure 3
Propagation of the numerical solutions for (44): (a) the implicit method, (b) the CN method, and (c) the variable-θ method, for 0≤t≤T=1.0, when Δt=0.01 and Δx=0.025.
(a)
(b)
(c)Figure 4
Numerical solutions foru1 of the Gierer–Meinhardt model at the spatial resolution Δx=Δy=1/32: (a) the variable-θ method with Δt=0.05, (b) the CN method with Δt=0.05, and (c) the CN method with Δt=0.005; (d–f) are the aerial views of (a–c), respectively.
(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:(45)fu=F1−u1−u1u22,u1u22−F+ku2T,for any constants F and k. Let Ω=0,1. We assign two sets of model constants and initial and boundary conditions as follows [28]:(46a)D=10−4,10−6T,F=0.035,k=0.049,(46b)u1x,0=1−12sin100πx,u2x,0=14sin100πx,x∈0,1,(46c)u10,t=u11,t=1,u20,t=u21,t=0,t∈0,T,where (46b) is a mid-pulse initial condition, and(47a)D=10−4,5⋅10−5T,F=0.025,k=0.0544,(47b)u1x,0=1−12cos100πx2,u2x,0=14cos100πx2,x∈0,1,(47c)∂u1∂x0,t=∂u1∂x1,t=0,∂u2∂x0,t=∂u2∂x1,t=0,t∈0,T,where (47b) is a left-pulse initial condition.In Figure5, we present the propagation of the numerical solution of u2 by the variable-θ method associated with (46a)–(46c) at the grid sizes Δt=0.01 and Δx=0.004 over 0≤t≤T=2000. 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.Figure 5
Propagation of the numerical solution ofu2 for the 1D Gray–Scott model by the variable-θ method: (a) the wave-splitting (self-replication of the pulse) and (b) its aerial view over 0≤t≤T=2000, when Δt=0.01 and Δx=0.004.
(a)
(b)Figure6 depicts the propagation of the numerical solution of u2 by the variable-θ method associated with (47a)–(47c) over 0≤t≤T=5000, 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).Figure 6
Propagation of the numerical solution ofu2 for the 1D Gray–Scott model by the variable-θ method: (a) the pulse traveling and reflecting and (b) its aerial view over 0≤t≤T=5000, when Δt=0.01 and Δx=0.004.
(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 Ikℓ represents the numerical solution obtained with the mesh resolution Δt,Δx=10−k,0.01/2ℓ−1. For example, the image I23 is associated with the mesh resolution Δt,Δx=1/100,1/400. 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).Figure 7
The wave-splitting (self-replication of the pulse) by the variable-θ method over 0≤t≤2000 with various Δt,Δx. The image Ikℓ represents the numerical solution obtained with the mesh resolution Δt,Δx=10−k,0.01/2ℓ−1. (a) (I11). (b) (I12). (c) (I13). (d) (I21). (e) (I22). (f) (I23). (g) (I31). (h) (I32). (i) (I33).
(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Δx’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Δx’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 fu given in (45). We choose problem coefficients as follows [15]:(48)Ω=−1,1×−1,1,D=0.001,0.001T,F=1,k=0.For the purpose of error analysis, we select a smooth solutionu^=u^1,u^2 defined as(49)u^1x,y,t=cos2tcos2πxcosπy,u^2x,y,t=cos2tcosπxcos2πy,and replace the reaction kinetics fu with fu^u:(50)fu^u≔fu+∂u^∂t−DΔu^−fu^.Then,u^=u^1,u^2 in (49) would be the exact solution of ∂u/∂t−DΔu−fu=fu^u with the initial condition u0=u^x,y,0.Table2 summarizes the L∞-error E∞T with T=1.0 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 Δx=Δy. Since the solution including the initial condition is smooth over the entire time interval 0,T=1, 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 T=1.0 for the Gray–Scott model, when Δt=0.025 and Δx=Δy=0.05.Table 2
L
∞-error E∞T and the elapsed time (CPU) for the numerical solution of the Gray–Scott model in the 2D space: T=1.0.
ε
=
10
−
6
Implicit
CN
Variable-θ
Δ
t
,
Δ
x
=
Δ
y
E
∞
T
CPU
E
∞
T
CPU
E
∞
T
CPU
u
1
u
2
u
1
u
2
u
1
u
2
0.05
,
0.1
4.3
⋅
10
−
2
4.0
⋅
10
−
2
0.036 s
9.6
⋅
10
−
4
9.4
⋅
10
−
4
0.047 s
9.6
⋅
10
−
4
9.4
⋅
10
−
4
0.062 s
0.025
,
0.05
2.0
⋅
10
−
2
1.7
⋅
10
−
2
0.225 s
2.5
⋅
10
−
4
2.3
⋅
10
−
4
0.280 s
2.5
⋅
10
−
4
2.3
⋅
10
−
4
0.394 s
0.0125
,
0.025
8.0
⋅
10
−
3
8.0
⋅
10
−
3
1.629 s
6.4
⋅
10
−
5
5.5
⋅
10
−
5
2.102 s
6.4
⋅
10
−
5
5.5
⋅
10
−
5
2.798 s
Conv. order
1.2
1.2
1.9
2.0
1.9
2.0Figure 8
The numerical solution of the variable-θ method and its error at T=1.0 for the Gray–Scott model, when Δt=0.025 and Δx=Δy=0.05: (a) u1, (b) u2, (c) e1=u1−u1^, and (d) e2=u2−u2^.
(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,Δt,Δx=Δy=0.05,0.1. That is, the constant c0 in (41) is evaluated using the experimentally optimal ω^Δt0,h0 with Δt0,h0=0.05,0.1 and then (40) is utilized to estimate ω^Δt,h for other grid sizes Δt,h=Δx=Δy. 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 u1 and u2) 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 Ω=−1,1×−1,1⊂ℝ2 with the following reaction kinetics and parameters:(51)D=ε2,κμT,fu=u12u2−u1,1μu12ε−u2T,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 Figure9. In this section, we restrict our attention to the dynamics of u1 of the model.(52)ε=0.04,μ=0.1,κ=0.0128,u1x,y,0=121+0.001∑k=120coskπy2sech2x2+y22ε,u2x,y,0=cosh1−x2+y23cosh1.Figure 9
The initial values in (52) for the Gierer–Meinhardt model with ϵ=0.04 at the mesh resolution Δx=Δy=1/32: (a) u1x,y,0 and (b) u2x,y,0.
(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 u1 of the Gierer–Meinhardt model with the spatial resolution Δx=Δy=1/32. When the time step size is set Δt=0.05, 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 Δt=0.05, 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 Δt=0.005, 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 Δt, 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 accuracyRemark 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.
Figure10 presents numerical solutions and their errors of u1 for the Gray–Scott model (48)–(50) at T=1.0 by the variable-θ method, varying Δx=Δy with fixed Δt=0.25. Compared with Figure 8, the solutions show stability and a good accuracy, although the time step size is as large as Δt=0.25. 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 Δt≥0.2, 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 un+1/2 is approximated by the extrapolation scheme (18). However, when Δt≤0.1, our method produces stable solutions for all choices of spatial grid sizes. Figure 11 presents numerical solutions of u1 for the Gierer–Meinhardt model at T=340 by the variable-θ method with fixed Δt=0.1 and various Δx=Δy. Note that for Gierer–Meinhardt model, the pattern forming is slow down as the spatial grid size becomes smaller, as shown in I2 and I3 of Figure 11; this tendency has been observed for all other choices of Δt≤0.1. This is another example that accuracy of the numerical solution is much more sensitive to the spatial mesh resolution than the temporal one.Figure 10
Numerical solutions foru1 of the Gray–Scott model at T=1.0 by the variable-θ method varying Δx=Δy with fixed Δt=0.25. The image Iℓ represents the numerical solution obtained with the mesh resolution Δt,Δx=Δy=0.25,0.02/2ℓ−1 and Iℓ′ represents the error of Iℓ. (a) (I1). (b) (I2). (c) (I3). (d)I1′. (e)I2′. (f)I3′.
(a)
(b)
(c)
(d)
(e)
(f)Figure 11
Numerical solutions foru1 of the Gierer–Meinhardt model at T=340 by the variable-θ method, varying Δx=Δy with fixed Δt=0.1. The image Iℓ represents the numerical solution obtained with the mesh resolution Δt,Δx=Δy=0.1,0.02/2ℓ−1 and Iℓ′ represents the aerial view of Iℓ. (a) (I1). (b) (I2). (c) (I3). (d)I1′. (e)I2′. (f)I3′.
(a)
(b)
(c)
(d)
(e)
(f)


## 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:(45)fu=F1−u1−u1u22,u1u22−F+ku2T,for any constants F and k. Let Ω=0,1. We assign two sets of model constants and initial and boundary conditions as follows [28]:(46a)D=10−4,10−6T,F=0.035,k=0.049,(46b)u1x,0=1−12sin100πx,u2x,0=14sin100πx,x∈0,1,(46c)u10,t=u11,t=1,u20,t=u21,t=0,t∈0,T,where (46b) is a mid-pulse initial condition, and(47a)D=10−4,5⋅10−5T,F=0.025,k=0.0544,(47b)u1x,0=1−12cos100πx2,u2x,0=14cos100πx2,x∈0,1,(47c)∂u1∂x0,t=∂u1∂x1,t=0,∂u2∂x0,t=∂u2∂x1,t=0,t∈0,T,where (47b) is a left-pulse initial condition.In Figure5, we present the propagation of the numerical solution of u2 by the variable-θ method associated with (46a)–(46c) at the grid sizes Δt=0.01 and Δx=0.004 over 0≤t≤T=2000. 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.Figure 5
Propagation of the numerical solution ofu2 for the 1D Gray–Scott model by the variable-θ method: (a) the wave-splitting (self-replication of the pulse) and (b) its aerial view over 0≤t≤T=2000, when Δt=0.01 and Δx=0.004.
(a)
(b)Figure6 depicts the propagation of the numerical solution of u2 by the variable-θ method associated with (47a)–(47c) over 0≤t≤T=5000, 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).Figure 6
Propagation of the numerical solution ofu2 for the 1D Gray–Scott model by the variable-θ method: (a) the pulse traveling and reflecting and (b) its aerial view over 0≤t≤T=5000, when Δt=0.01 and Δx=0.004.
(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 Ikℓ represents the numerical solution obtained with the mesh resolution Δt,Δx=10−k,0.01/2ℓ−1. For example, the image I23 is associated with the mesh resolution Δt,Δx=1/100,1/400. 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).Figure 7
The wave-splitting (self-replication of the pulse) by the variable-θ method over 0≤t≤2000 with various Δt,Δx. The image Ikℓ represents the numerical solution obtained with the mesh resolution Δt,Δx=10−k,0.01/2ℓ−1. (a) (I11). (b) (I12). (c) (I13). (d) (I21). (e) (I22). (f) (I23). (g) (I31). (h) (I32). (i) (I33).
(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Δx’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Δx’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 fu given in (45). We choose problem coefficients as follows [15]:(48)Ω=−1,1×−1,1,D=0.001,0.001T,F=1,k=0.For the purpose of error analysis, we select a smooth solutionu^=u^1,u^2 defined as(49)u^1x,y,t=cos2tcos2πxcosπy,u^2x,y,t=cos2tcosπxcos2πy,and replace the reaction kinetics fu with fu^u:(50)fu^u≔fu+∂u^∂t−DΔu^−fu^.Then,u^=u^1,u^2 in (49) would be the exact solution of ∂u/∂t−DΔu−fu=fu^u with the initial condition u0=u^x,y,0.Table2 summarizes the L∞-error E∞T with T=1.0 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 Δx=Δy. Since the solution including the initial condition is smooth over the entire time interval 0,T=1, 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 T=1.0 for the Gray–Scott model, when Δt=0.025 and Δx=Δy=0.05.Table 2
L
∞-error E∞T and the elapsed time (CPU) for the numerical solution of the Gray–Scott model in the 2D space: T=1.0.
ε
=
10
−
6
Implicit
CN
Variable-θ
Δ
t
,
Δ
x
=
Δ
y
E
∞
T
CPU
E
∞
T
CPU
E
∞
T
CPU
u
1
u
2
u
1
u
2
u
1
u
2
0.05
,
0.1
4.3
⋅
10
−
2
4.0
⋅
10
−
2
0.036 s
9.6
⋅
10
−
4
9.4
⋅
10
−
4
0.047 s
9.6
⋅
10
−
4
9.4
⋅
10
−
4
0.062 s
0.025
,
0.05
2.0
⋅
10
−
2
1.7
⋅
10
−
2
0.225 s
2.5
⋅
10
−
4
2.3
⋅
10
−
4
0.280 s
2.5
⋅
10
−
4
2.3
⋅
10
−
4
0.394 s
0.0125
,
0.025
8.0
⋅
10
−
3
8.0
⋅
10
−
3
1.629 s
6.4
⋅
10
−
5
5.5
⋅
10
−
5
2.102 s
6.4
⋅
10
−
5
5.5
⋅
10
−
5
2.798 s
Conv. order
1.2
1.2
1.9
2.0
1.9
2.0Figure 8
The numerical solution of the variable-θ method and its error at T=1.0 for the Gray–Scott model, when Δt=0.025 and Δx=Δy=0.05: (a) u1, (b) u2, (c) e1=u1−u1^, and (d) e2=u2−u2^.
(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,Δt,Δx=Δy=0.05,0.1. That is, the constant c0 in (41) is evaluated using the experimentally optimal ω^Δt0,h0 with Δt0,h0=0.05,0.1 and then (40) is utilized to estimate ω^Δt,h for other grid sizes Δt,h=Δx=Δy. 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 u1 and u2) 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 Ω=−1,1×−1,1⊂ℝ2 with the following reaction kinetics and parameters:(51)D=ε2,κμT,fu=u12u2−u1,1μu12ε−u2T,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 Figure9. In this section, we restrict our attention to the dynamics of u1 of the model.(52)ε=0.04,μ=0.1,κ=0.0128,u1x,y,0=121+0.001∑k=120coskπy2sech2x2+y22ε,u2x,y,0=cosh1−x2+y23cosh1.Figure 9
The initial values in (52) for the Gierer–Meinhardt model with ϵ=0.04 at the mesh resolution Δx=Δy=1/32: (a) u1x,y,0 and (b) u2x,y,0.
(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 u1 of the Gierer–Meinhardt model with the spatial resolution Δx=Δy=1/32. When the time step size is set Δt=0.05, 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 Δt=0.05, 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 Δt=0.005, 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 Δt, 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 accuracyRemark 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.
Figure10 presents numerical solutions and their errors of u1 for the Gray–Scott model (48)–(50) at T=1.0 by the variable-θ method, varying Δx=Δy with fixed Δt=0.25. Compared with Figure 8, the solutions show stability and a good accuracy, although the time step size is as large as Δt=0.25. 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 Δt≥0.2, 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 un+1/2 is approximated by the extrapolation scheme (18). However, when Δt≤0.1, our method produces stable solutions for all choices of spatial grid sizes. Figure 11 presents numerical solutions of u1 for the Gierer–Meinhardt model at T=340 by the variable-θ method with fixed Δt=0.1 and various Δx=Δy. Note that for Gierer–Meinhardt model, the pattern forming is slow down as the spatial grid size becomes smaller, as shown in I2 and I3 of Figure 11; this tendency has been observed for all other choices of Δt≤0.1. This is another example that accuracy of the numerical solution is much more sensitive to the spatial mesh resolution than the temporal one.Figure 10
Numerical solutions foru1 of the Gray–Scott model at T=1.0 by the variable-θ method varying Δx=Δy with fixed Δt=0.25. The image Iℓ represents the numerical solution obtained with the mesh resolution Δt,Δx=Δy=0.25,0.02/2ℓ−1 and Iℓ′ represents the error of Iℓ. (a) (I1). (b) (I2). (c) (I3). (d)I1′. (e)I2′. (f)I3′.
(a)
(b)
(c)
(d)
(e)
(f)Figure 11
Numerical solutions foru1 of the Gierer–Meinhardt model at T=340 by the variable-θ method, varying Δx=Δy with fixed Δt=0.1. The image Iℓ represents the numerical solution obtained with the mesh resolution Δt,Δx=Δy=0.1,0.02/2ℓ−1 and Iℓ′ represents the aerial view of Iℓ. (a) (I1). (b) (I2). (c) (I3). (d)I1′. (e)I2′. (f)I3′.
(a)
(b)
(c)
(d)
(e)
(f)


## 6. Conclusions

TheCrank–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.



---

*Source: 5163704-2020-04-30.xml*