A Time-Split MacCormack Scheme for Two-Dimensional Nonlinear Reaction-Diffusion Equations

TL;DR

This paper introduces a three-level explicit time-split MacCormack scheme for efficiently solving two-dimensional nonlinear reaction-diffusion equations, demonstrating stability, convergence, and improved computational cost.

Contribution

It presents a novel time-split MacCormack scheme with proven stability and convergence for nonlinear reaction-diffusion equations, reducing computational costs.

Findings

The scheme is stable under CFL conditions.

Convergence rate is validated through numerical experiments.

The method reduces computational costs compared to traditional schemes.

Abstract

A three-level explicit time-split MacCormack scheme is proposed for solving the two-dimensional nonlinear reaction-diffusion equations. The computational cost is reduced thank to the splitting and the explicit MacCormack scheme. Under the well known condition of Courant-Friedrich-Lewy (CFL) for stability of explicit numerical schemes applied to linear parabolic partial differential equations, we prove the stability and convergence of the method in $L^{\infty}(0,T;L^{2})$-norm. A wide set of numerical evidences which provide the convergence rate of the new algorithm are presented and critically discussed.

Equations

$$u_t-a\Delta u=f(u),(x,y)\in \Omega ,t\in (0,T];$$

$$u(x,y,0)=u_0(x,y),(x,y)\in \overline{\Omega };$$

$$u(x,y,t)=\phi (x,y,t),(x,y)\in \partial \Omega ,t\in (0,T];$$

$$u_{ij}^{\ast }=L_x(\Delta t_x)u_{ij}^n,$$

$$u_{ij}^{\ast }=L_y(\Delta t_y)u_{ij}^n.$$

$$u_{ij}^{n+1}={\left[L_y\left(\frac{\Delta t}{2m}\right)\right]}^m{L}_x(\Delta t){\left[L_y\left(\frac{\Delta t}{2m}\right)\right]}^m{u}_{ij}^n.$$

$${\delta }_t{u}_{ij}^{\ast }=\frac{u_{ij}^{\ast }-u_{ij}^n}{k/2};{\delta }_t{u}_{ij}^{\ast \ast }=\frac{u_{ij}^{\ast \ast }-u_{ij}^{\ast }}{k};{\delta }_t{u}_{ij}^{n+1}=\frac{u_{ij}^{n+1}-u_{ij}^{\ast \ast }}{k/2};{\delta }_x{u}_{i+\frac{1}{2},j}^n=\frac{u_{i+1,j}^n-u_{ij}^n}{h};$$

$${\delta }_y{u}_{i,j+\frac{1}{2}}^n=\frac{u_{i,j+1}^n-u_{ij}^n}{h};{\delta }_x^2{u}_{ij}^n=\frac{{\delta }_x{u}_{i+\frac{1}{2},j}^n-{\delta }_x{u}_{i-\frac{1}{2},j}^n}{h};{\delta }_y^2{u}_{ij}^n=\frac{{\delta }_y{u}_{i,j+\frac{1}{2}}^n-{\delta }_y{u}_{i,j-\frac{1}{2}}^n}{h}.$$

$$\Vert u^n{\Vert }_{L^2}=h{\left(\underset{i,j=1}{\sum ^{M-1}}|u_{ij}^n{|}^2\right)}^{\frac{1}{2}};\Vert {\delta }_x{u}^n{\Vert }_{L^2}=h{\left(\underset{j=1}{\sum ^{M-1}}\underset{i=0}{\sum ^{M-1}}|{\delta }_x{u}_{i+\frac{1}{2},j}^n{|}^2\right)}^{\frac{1}{2}};$$

$$\Vert {\delta }_y{u}^n{\Vert }_{L^2}=h{\left(\underset{j=0}{\sum ^{M-1}}\underset{i=1}{\sum ^{M-1}}|{\delta }_y{u}_{i,j+\frac{1}{2}}^n{|}^2\right)}^{\frac{1}{2}};\Vert {\delta }_{\lambda }^2{u}^n{\Vert }_{L^2}=h{\left(\underset{i,j=1}{\sum ^{M-1}}|{\delta }_{\lambda }^2{u}_{ij}^n{|}^2\right)}^{\frac{1}{2}},$$

$$(u^n,v^n)=h^2\underset{i,j=1}{\sum ^{M-1}}{u}_{ij}^n{v}_{ij}^n;<{\delta }_x{u}^n,{\delta }_x{v}^n{>}_x=h^2\underset{j=1}{\sum ^{M-1}}\underset{i=0}{\sum ^{M-1}}{\delta }_x{u}_{i+\frac{1}{2},j}^n{\delta }_x{v}_{i+\frac{1}{2},j}^n;$$

$$<{\delta }_y{u}^n,{\delta }_y{v}^n{>}_y=h^2\underset{j=0}{\sum ^{M-1}}\underset{i=1}{\sum ^{M-1}}{\delta }_y{u}_{i,j+\frac{1}{2}}^n{\delta }_y{v}_{i,j+\frac{1}{2}}^n.$$

$$|u^n{|}_{H^1}={\left(\Vert {\delta }_x{u}^n{\Vert }_{L^2}^2+\Vert {\delta }_y{u}^n{\Vert }_{L^2}^2\right)}^{\frac{1}{2}}\text{and}\Vert u^n{\Vert }_{H^1}={\left(\Vert u^n{\Vert }_{L^2}^2+\Vert {\delta }_x{u}^n{\Vert }_{L^2}^2+\Vert {\delta }_y{u}^n{\Vert }_{L^2}^2\right)}^{\frac{1}{2}}.$$

$$u_{ij}^{n+1}=L_y(k/2)L_x(k)L_y(k/2)u_{ij}^n.$$

$$u_{ij}^{\ast }=L_y(k/2)u_{ij}^n;u_{ij}^{\ast \ast }=L_x(k)u_{ij}^{\ast }=L_x(k)L_y(k/2)u_{ij}^n\text{and}{u}_{ij}^{n+1}=L_y(k/2)u_{ij}^{\ast \ast }.$$

$$u_{ij}^{\overline{\ast }}=u_{ij}^n+\frac{k}{2}{u}_t{)}_{ij}^n+\frac{k^2}{8}{u}_{2t}{)}_{ij}^n+O(k^3);u_{ij}^{\overline{\overline{\ast }}}=u_{ij}^n+\frac{k}{2}{u}_t{)}_{ij}^{\overline{\ast }}+\frac{k^2}{8}{u}_{2t}{)}_{ij}^{\overline{\ast }}+O(k^3).$$

$$u_t-a{u}_{yy}=f(u),\text{whichisequivalentto}{u}_t=a{u}_{yy}+f(u).$$

$$u_{2t}={\left(a{u}_{yy}+f(u)\right)}_t=a^2{u}_{4y}+a{\left(f(u)\right)}_{yy}+\left(a{u}_{yy}+f(u)\right)f^{{}^{\prime }}(u).$$

$$u_{ij}^{\overline{\ast }}=u_{ij}^n+\frac{k}{2}[a{u}_{yy}+f(u){]}_{ij}^n+\frac{k^2}{8}{\left[a^2{u}_{4y}+a{\left(f(u)\right)}_{yy}+\left(a{u}_{yy}+f(u)\right)f^{{}^{\prime }}(u)\right]}_{ij}^n+O(k^3);$$

$$u_{ij}^{\overline{\overline{\ast }}}=u_{ij}^n+\frac{k}{2}[a{u}_{yy}+f(u){]}_{ij}^{\overline{\ast }}+\frac{k^2}{8}{\left[a^2{u}_{4y}+a{\left(f(u)\right)}_{yy}+\left(a{u}_{yy}+f(u)\right)f^{{}^{\prime }}(u)\right]}_{ij}^{\overline{\ast }}+O(k^3).$$

$$u_{yy,ij}^n={\delta }_y^2{u}_{ij}^n+O(h^2);{\left(f(u)\right)}_{yy,ij}^n={\delta }_y^2\left(f(u_{ij}^n)\right)+O(h^2);u_{4y,ij}^n={\delta }_y^2({\delta }_y^2{u}_{ij}^n)+O(h^2);$$

$$u_{yy,ij}^{\overline{\ast }}={\delta }_y^2{u}_{ij}^{\overline{\ast }}+O(h^2);{\left(f(u)\right)}_{yy,ij}^{\overline{\ast }}={\delta }_y^2\left(f(u_{ij}^{\overline{\ast }})\right)+O(h^2);u_{4y,ij}^{\overline{\ast }}={\delta }_y^2({\delta }_y^2{u}_{ij}^{\overline{\ast }})+O(h^2),$$

$$u_{ij}^{\overline{\ast }}=u_{ij}^n+\frac{k}{2}[a{\delta }_y^2{u}_{ij}^n+f(u_{ij}^n)]+k^2{\rho }_{ij}^n+O(k^3+k{h}^2);$$

$$u_{ij}^{\overline{\overline{\ast }}}=u_{ij}^n+\frac{k}{2}[a{\delta }_y^2{u}_{ij}^{\overline{\ast }}+f(u_{ij}^{\overline{\ast }})]+k^2{\rho }_{ij}^{\overline{\ast }}+O(k^3+k{h}^2),$$

$${\rho }_{ij}^{\alpha }=\frac{1}{8}\left[a^2{\delta }_y^2({\delta }_y^2{u}_{ij}^{\alpha })+a{\delta }_y^{2}f(u_{ij}^{\alpha })+a{\delta }_y^2{u}_{ij}^{\alpha }{f}^{{}^{\prime }}(u_{ij}^{\alpha })+f(u_{ij}^{\alpha })f^{{}^{\prime }}(u_{ij}^{\alpha })\right],$$

$$f(u_{ij}^{\overline{\ast }})=f(u_{ij}^n)+\frac{k}{2}{u}_{t,ij}^n{f}^{{}^{\prime }}(u_{ij}^n)+O(k^2).$$

$$u_{t,ij}^n=a{u}_{yy,ij}^n+f(u_{ij}^n)=a{\delta }_y^2{u}_{ij}^n+f(u_{ij}^n)+O(h^2).$$

$$f(u_{ij}^{\overline{\ast }})=f(u_{ij}^n)+\frac{k}{2}\left[a{\delta }_y^2{u}_{ij}^n+f(u_{ij}^n)\right]f^{{}^{\prime }}(u_{ij}^n)+O(k^2+k{h}^2).$$

$$u_{ij}^{\overline{\overline{\ast }}}=u_{ij}^n+\frac{k}{2}\left(a{\delta }_y^2{u}_{ij}^n+f(u_{ij}^n)\right)+k^2(2{\rho }_{ij}^n+{\rho }_{ij}^{\overline{\ast }})+O(k^3+k{h}^2).$$

$$\frac{u_{ij}^{\overline{\ast }}+u_{ij}^{\overline{\overline{\ast }}}}{2}=u_{ij}^n+\frac{k}{2}\left(a{\delta }_y^2{u}_{ij}^n+f(u_{ij}^n)\right)+\frac{1}{2}{k}^2(3{\rho }_{ij}^n+{\rho }_{ij}^{\overline{\ast }})+O(k^3+k{h}^2),$$

Full text

A Time-Split MacCormack Scheme for Two-Dimensional Nonlinear Reaction-Diffusion Equations

Eric Ngondiep Tel.: +966506048689. E-mail addresses: ericngondiep@gmail.com or engondiep@imamu.edu.sa (Eric Ngondiep). 

(Department of Mathematics and Statistics, College of Science, Al-Imam Muhammad Ibn Saud Islamic University (IMSIU), 90950 Riyadh, Kingdom of Saudi Arabia

 

Hydrological Research Centre, Institute for Geological and Mining Research, 4110 Yaounde-Cameroon)

Abstract. A three-level explicit time-split MacCormack scheme is proposed for solving the two-dimensional nonlinear reaction-diffusion equations. The computational cost is reduced thank to the splitting and the explicit MacCormack scheme. Under the well known condition of Courant-Friedrich-Lewy (CFL) for stability of explicit numerical schemes applied to linear parabolic partial differential equations, we prove the stability and convergence of the method in $L^{\infty}(0,T;L^{2})$-norm. A wide set of numerical evidences which provide the convergence rate of the new algorithm are presented and critically discussed.  

 

 Keywords: $2$D nonlinear reaction-diffusion equations, locally one-dimensional operators (splitting), explicit MacCormack scheme, a three-level explicit time-split MacCormack method, stability and convergence rate.

AMS Subject Classification (MSC). 65M10, 65M05.

1 Introduction and motivation

A large number of biological problems of significant interest are modeled by parabolic equations [9]. The general framework is a set of biological entities (either ions, molecules, proteins or cells) that interact with each other and diffuse within a given domain. So it becomes possible to build some models via reaction-diffusion equations. For example, the dendritic spines possess a twitching motion which are described by the reaction-diffusion models [12]. In this paper, we consider the following two-dimensional reaction-diffusion equations,

[TABLE]

with the initial condition

[TABLE]

and the boundary condition

[TABLE]

where $a$ is the diffusive coefficient, $f\in\mathcal{C}^{1}(\mathbb{R})$ is a Lipschitz function, $\Omega=(0,1)^{2},$ $\Delta$ denotes the Laplacian operator, $\partial\Omega$ is the boundary of $\Omega$ and $u_{t}$ designates $\frac{\partial u}{\partial t}$. The initial condition $u_{0}$ and the boundary condition $\varphi$ are assumed to be regular enough and satisfy the requirement $\varphi(x,y,0)=u_{0}(x,y),$ for every $(x,y)\in\partial\Omega,$ so that the initial value problem $(\ref{1})$-$(\ref{3}),$ admits a smooth solution.

In the last decades [23, 19, 32], MacCormack approach which is a predictor-corrector, finite difference scheme has been used to solve certain classes of nonlinear partial differential equations (PDEs). There exist both explicit and implicit versions of the method, but the explicit predates the implicit by more than a decade, and it is considered as one of the milestones of computational fluid dynamics. Both versions facilitate the solution of parabolic and hyperbolic equations by marching forward in time [19, 20, 21]. The popularity of MacCormack explicit method is due in part to its simplicity and ease of implementation. The predictor and corrector phases each uses forward differencing for first-order time derivatives, with alternate one-side differencing for first-order space derivatives. This is especially convenient for systems of equations with nonlinear advertive jacobian matrices associated with one-side explicit schemes, such as Lax-Wendroff approach (for instance, see [13, 25, 32]). However, the explicit MacCormack is not a suitable method for solving high Reynolds numbers flows, where the viscous regions become very thin (see [2], P. $630$). To overcome this difficulty, MacCormack [18] developed a hybrid version of his scheme, known as the MacCormack rapid solver method. The new algorithm is an explicit-implicit method. For example, in a search of an efficient solution, the authors [26, 27, 31, 29] applied this hybrid method to some complex PDEs (such as: mixed Stokes-Darcy model and $2D$ incompressible Navier-Stokes equations) and they obtained satisfactory results regarding both stability and convergence rate of the method. It is worth noticing to mention that the rapid solver algorithm has a good stability condition and it is too faster than a large set of numerical methods for solving steady and unsteady flows at high to low Reynolds numbers [18]. So, the hybrid method of MacCormack will be used to solve the $2$D reaction-diffusion equations $(\ref{1})$-$(\ref{3}),$ in our future works.

Armed with the information gleaned from both MacCormack and MacCormack rapid solver methods, we can now analyze a time-split MacCormack technique applied to problem $(\ref{1})$-$(\ref{3}).$ Firstly, it’s worth noting to recall that the problem considered in this paper has been solved in literature by a wide set of explicit, implicit and coupled explicit-implicit numerical schemes. While some explicit methods usually suffer the severely restricted temporal step size [17, 35], the fully implicit methods although unconditionally stable, provide a large system of nonlinear equations at every time level [3, 15]. These systems lead to a considerable computational cost in practical applications. A possible improvement is to use the second time discretization such as the linearized Cank-Nicolson, implicit-explicit and collocation approaches [10, 14, 5, 30, 34, 33, 38, 9, 16]. Unfortunately, at least two starting values are needed to begin these algorithms. These values can be obtained by the initial and boundary conditions and an additional iterative method. To overcome this difficulty, the authors [37] applied a two-level linearized compact ADI scheme to problem $(\ref{1})$-$(\ref{3}).$ The main results of their work (namely Theorems $1$-$2$) have been proved under the assumptions that the time step $\Delta t=k,$ mesh size $\Delta x=\Delta y=h$ and the ratio $\frac{\Delta t}{h^{2}}\leq C,$ must be sufficiently small ([37], page $9$, line above $(3.26)$ and page $10,$ Theorem $2$). These requirements in general are less restrictive and can make the method even more impractical. Furthermore, this paper represents an extension of the work in [28].

The time-split MacCormack approach we study for the initial-boundary value problem $(\ref{1})$-$(\ref{3})$ is new, a three-level explicit predictor-corrector method, second order accurate in time and fourth order convergent in space, under the time step restriction: $\frac{ak}{h^{2}}\leq\frac{1}{2},$ and it is motivated by this time step restriction (indeed, lots of explicit schemes for solving equation $(\ref{1})$-$(\ref{3}),$ are stable under the well-known condition of Courant-Friedrich-Lewy: $\frac{4ak}{h^{2}}\leq 1$) and its efficiency and effectiveness. From this observation, it is obvious that: (a) a time-split MacCormack approach is more practical, (b) although the new algorithm and a two-level linearized compact ADI method have the same order of convergence, the linearized compact ADI scheme requires substantially more computer times to solve problem $(\ref{1})$-$(\ref{3}),$ than does a time-split MacCormack. An explicit time-split MacCormack algorithm [21, 28, 22] ”splits” the original MacCormack scheme into a sequence of one-dimensional operations, thereby achieving a good stability condition. In other words, the splitting makes it possible to advance the solution in each direction with the maximum allowable time step. This is particularly advantageous if the allowable time steps $\Delta t_{x}$ and $\Delta t_{y}$, are much different because of differences in the mesh spacings $\Delta x$ and $\Delta y.$ In order to explain this method, we will make use of the $1$D difference operators $L_{x}(\Delta t_{x})$ and $L_{y}(\Delta t_{y})$. Setting $u_{ij}^{n}=u(x_{i},y_{j},t^{n}),$ the $L_{x}(\Delta t_{x})$ operator applied to $u_{ij}^{n}$,

[TABLE]

is by definition equivalent to the two-step predictor-corrector MacCormack formulation. The $L_{y}(\Delta t_{y})$ operator is defined in a similar manner, that is,

[TABLE]

These expressions make use of a dummy time index, which is denoted by the asterisk. Now, letting $\Delta t_{x}=\Delta t$ and $\Delta t_{y}=\frac{\Delta t}{2m},$ where $m$ is a positive integer, a second order accurate scheme can be constructed by applying the $L_{x}$ and $Ly$ operators to $u_{ij}^{n},$ in the following way:

[TABLE]

This sequence is quite useful for the case $\Delta y<<\Delta x.$ In general, a scheme formed by a sequence of these operators is: $(1)$ stable, if the time step of each operator does not exceed the allowable time step for that operator; $(2)$ consistent, if the sums of the time steps for each of the operators are equal: and $(3)$ second-order accurate, if the sequence is symmetric.

In this paper, we are interested in a numerical solution of the initial-boundary value problem $(\ref{1})$-$(\ref{3}),$ using a time-split MacCormack approach. Specifically, the work is focused on the following four items:

1.

full description of a three-level explicit time-split MacCormack scheme for solving the nonlinear reaction-diffusion equations $(\ref{1})$-$(\ref{3});$

2.

stability analysis of the numerical scheme;

3.

error estimates of the method;

4.

a wide set of numerical examples which provide the convergence rate, confirms the theoretical results and shows the efficiency and effectiveness of the method.

Items $1$, $2$ and $3$, are our original contributions since as far as we know, there is no available work in literature which solves the reaction-diffusion model $(\ref{1})$-$(\ref{3}),$ using a time-split MacCormack method.

The paper is organized as follows: Section 2 considers a detailed description of a three-level explicit time-split MacCormack method applied to problem $(\ref{1})$-$(\ref{3}).$ In section $\ref{sec3},$ we study the stability of the numerical scheme under the condition given above, while section 4 analyzes the error estimates and the convergence of the method. A large set of numerical examples which provides the convergence rate of the new algorithm and confirms the theoretical result (on the stability) are presented and discussed in section $\ref{sec5}.$ We draw the general conclusion and present our future works in section $\ref{sec6}.$

2 Full description of a time-split MacCormack method

This section deals with the description of a three-level explicit time-split MacCormack method applied to two-dimensional nonlinear reaction-diffusion equations $(\ref{1})$-$(\ref{3}).$

Let $N$ and $M$ be two positive integers. Set $k:=\Delta t=\frac{T}{N};$ $h:=\Delta x=\Delta y=\frac{1}{M},$ be the time step and mesh size, respectively. Put $t^{n}=kn,$ $t^{*}=(n+r)k,$ $t^{**}=(n+s)k,$ where $0<r<s<1,$ so $t^{*}\in(t^{n},t^{n+1}),$ $t^{**}\in(t^{*},t^{n+1});$ $n=0,1,2,...,N-1;$ $x_{i}=ih;$ $y_{j}=jh;$ $0\leq i,j\leq M$. Also, let $\Omega_{k}=\{t^{n},0\leq n\leq N\};$ $\overline{\Omega}_{h}=\{(x_{i},y_{j}),0\leq i,j\leq M\};$ $\Omega_{h}=\overline{\Omega}_{h}\cap\Omega$ and $\partial\Omega_{h}=\overline{\Omega}_{h}\cap\partial\Omega.$

Consider $\mathcal{U}_{h}=\{u_{ij}^{n},n=0,1,...,N;\text{\,}i,j=0,1,2,...,M\}$ be the space of grid functions defined on $\Omega_{h}\times\Omega_{k}.$ Letting

[TABLE]

[TABLE]

Using this, we define the following norms and scalar products.

[TABLE]

[TABLE]

where $\lambda=x$ or $y.$ Furthermore, the scalar products are defined as

[TABLE]

and

[TABLE]

The space $H^{1}(\Omega)$ is endowed with the norm $|\cdot|_{H^{1}}$ (respectively, $\|\cdot\|_{H^{1}}$) defined as

[TABLE]

It is worth noticing to mention that a time-split MacCormack [22, 21] ”splits” the original MacCormack scheme into a sequence of $1$D operators, thereby achieving a less restrictive stability condition. In order words, the splitting makes it possible to advance the solution in each direction with the maximum allowable time step ([2], page $231$).

In other to give a detailed description of this method, we consider the $1$D difference operators $L_{x}(\Delta t_{x})$ and $L_{y}(\Delta t_{y})$ defined by equations $(\ref{8})$ and $(\ref{9})$, respectively. Following the approach presented in ([2], page $231$), a second-order accurate scheme can be constructed by applying the $L_{x}$ and $L_{y}$ operators to $u_{ij}^{n}$ in the following manner:

[TABLE]

Using these tools, we are able to provide a three-level explicit time-split MacCormack method for solving the initial-boundary value problem $(\ref{1})$-$(\ref{3})$. Putting $\Delta t_{x}=k,$ $\Delta t_{y}=\frac{k}{2}$ and $\Delta x=\Delta y:=h,$ it comes from equations $(\ref{8})$, $(\ref{9})$ and $(\ref{10})$ that

[TABLE]

In the following, we should find explicit expressions of equations $u_{ij}^{*}=L_{y}(k/2)u_{ij}^{n}$ and $u_{ij}^{**}=L_{x}(k)u_{ij}^{*}.$ This will help to give an explicit formula of the equation $u_{ij}^{n+1}=L_{y}(k/2)u_{ij}^{**},$ which represents a ”one-step time-split MacCormack algorithm”. For the sake of simplicity, we use both notations: $u)_{ij}^{n}=u_{ij}^{n}$ and $[u+v]_{ij}^{n}=u_{ij}^{n}+v^{n}_{ij}$.

The application of the Taylor series expansion about $(x_{i},y_{j},t^{n})$ at the predictor and corrector steps with time step $k/2$ yields

[TABLE]

From the definition of the operator $L_{y}(k/2),$ let consider the equation

[TABLE]

Using equation $(\ref{11a}),$ it is not difficult to see that

[TABLE]

This fact together with equation $(\ref{12a})$ provide

[TABLE]

[TABLE]

Now, expanding the Taylor series about $(x_{i},y_{j},t^{n})$ with mesh size $h$ using central difference representation to get

[TABLE]

[TABLE]

where $\delta_{y}^{2}w^{l}_{ij}$ is given by relation $(\ref{4}).$ Substituting equations $(\ref{14})$ into equations $(\ref{12})$ and $(\ref{13})$ to obtain

[TABLE]

and

[TABLE]

where

[TABLE]

where $\alpha=n,\overline{*}.$ The term $f(u_{ij}^{\overline{*}})$ should be expressed as a function of $f(u_{ij}^{n}),$ $f^{{}^{\prime}}(u_{ij}^{n})$ and $u_{t,ij}^{n}.$ Applying the Taylor expansion about $(x_{i},y_{j},t^{n})$ with time step $k/2$ using forward difference representation to get

[TABLE]

But, it comes from equation $(\ref{11a})$ and relations $(\ref{14})$ that

[TABLE]

This fact, together with equation $(\ref{17})$ result in

[TABLE]

Plugging equations $(\ref{15}),$ $(\ref{16})$ and $(\ref{18}),$ straightforward computations give

[TABLE]

Taking the average of $u_{ij}^{\overline{*}}$ and $u_{ij}^{\overline{\overline{*}}}$ to get

[TABLE]

where $\rho_{ij}^{\alpha}$ is given by relation $(\ref{20}).$

On the other hand, to define the operator $L_{x}(k),$ we should consider the equation

[TABLE]

It comes from equation $(\ref{21}),$ that

[TABLE]

Applying the Taylor series expansion about $(x_{i},y_{i},t^{*})$ (where $t^{*}\in(t^{n},t^{n+1}),$ is the time used at the beginning of the next step in a time-split MacCormack scheme) with mesh size $h$ using central difference representation, we obtain

[TABLE]

where $\delta_{x}^{2}u^{l}_{ij}$ is defined by equation $(\ref{4}).$ Also, expanding the Taylor series at the predictor and corrector steps about $(x_{i},y_{j},t^{*})$ with time step $k$ using forward difference, it is not difficult to observe that

[TABLE]

A combination of equations $(\ref{23})$, $(\ref{22})$, $(\ref{21})$ and $(\ref{21a})$ provides

[TABLE]

In order to obtain a simple expression of $\delta_{x}^{2}u^{\overline{**}}_{ij}$, we should use the first equation in $(\ref{25})$. Tracking the infinitesimal term in this equation, direct computations give

[TABLE]

The truncation of this error term does not compromise the result. This fact, together with relation $(\ref{25})$ yield

[TABLE]

Taking the average of $u_{ij}^{\overline{**}}$ and $u_{ij}^{\overline{\overline{**}}}$, it is not hard to see that

[TABLE]

In way similar, starting with the one-dimensional equation: $u_{t}-au_{yy}=f(u)$, expanding the Taylor series about $(x_{i},y_{j},t^{**})$ (where $t^{**}$ represents the time used at the last step in a time-split MacCormack approach) at the predictor and corrector steps with time step $k/2$ and mesh size $h,$ using forward difference representations to get

[TABLE]

where we set $\mu=**,\overline{n+1}$, and

[TABLE]

To construct a three-level explicit time-split MacCormack method for solving the nonlinear reaction-diffusion equation $(\ref{1})$-$(\ref{3}),$ we must follow the ideas presented in the literature to construct the explicit MacCormack scheme[18, 19, 21, 22]. Specifically, we should neglect the terms of second order together with the infinitesimal term $O(k^{3}+kh^{2})$ in equations $(\ref{20a}),$ $(\ref{31})$ and $(\ref{33})$. In addition, the terms $u_{ij}^{*},$ $u_{ij}^{**}$ and $u^{n+1}_{ij}$ must be defined as the average of predicted and corrected values, that is,

[TABLE]

Thus, equations

[TABLE]

are by definition equivalent to

[TABLE]

Since the operator $L_{y}(k/2)L_{x}(k)L_{y}(k/2)$ is symmetric, this fact together with relations $(\ref{20a})$, $(\ref{31})$ and $(\ref{33})$ show that the obtained method is a three-level technique, an explicit predictor-corrector scheme, second order accurate in time and fourth order convergent in space. This theoretical result is confirmed by a wide set of numerical examples (we refer the readers to section 5). From the definition of the linear operators $"\delta_{x}^{2}"$ and $"\delta_{y}^{2}"$ given by $(\ref{4})$, equation $(\ref{37})$ can be rewritten as follows. For $n=0,1,...,N-1;$

[TABLE]

[TABLE]

[TABLE]

with the initial and boundary conditions. For $i,j=0,1,...,M,$

[TABLE]

[TABLE]

[TABLE]

which represent a detailed description of a three-level explicit time-split MacCormack method applied to problem $(\ref{1})$-$(\ref{3}).$

In the rest of this paper, we prove the stability, the error estimates and the convergence rate of a three-level time-split MacCormack approach under the time step restriction

[TABLE]

where $a$ is the diffusive coefficient given in equation $(\ref{1}).$ We recall that $k$ is the time step and $h$ is the grid size. Estimate $(\ref{44})$ is well known in literature as CFL condition for stability of the explicit schemes when solving linear parabolic equations. We assume that the analytical solution $\overline{u}\in L^{\infty}(0,T;L^{2}(\Omega))\cap H^{1}(0,T;H^{3}(\Omega))\cap H^{2}(0,T;H^{1}(\Omega))\cap H^{2}(0,T;L^{2}(\Omega))\cap L^{2}(0,T;H^{4}(\Omega)),$ that is, there exists a positive constant $\widetilde{C},$ independent of both time step $k$ and mesh size $h,$ so that

[TABLE]

3 Stability analysis of a three-level time-split MacCormack scheme

This section considers a deep analysis of the stability of a three-level time-split MacCormack scheme $(\ref{40})$-$(\ref{24})$ for solving equations $(\ref{1})$-$(\ref{3}).$

Theorem 3.1**.**

Let $u$ be the solution provided by the numerical scheme $(\ref{40})$-$(\ref{24})$. Under the time step restriction $(\ref{44}),$ the following estimate holds

[TABLE]

where $C$ is a positive constant independent of the time step $k$ and mesh size $h$ and $\widetilde{C}$ is given by relation $(\ref{43}).$

The following result (namely Lemmas 3.1) plays a crucial role in the proof of Theorem 3.1.

Lemma 3.1**.**

Setting $u_{ij}^{n}=u(x_{i},y_{j},t^{n})$ be the numerical solution provided by the scheme $(\ref{40})$-$(\ref{24})$, $\overline{u}_{ij}^{n}=\overline{u}(x_{i},y_{j},t^{n})$ be the exact one and let $e_{ij}^{n}=u_{ij}^{n}-\overline{u}_{ij}^{n}$ be the error. We recall that $\overline{u}_{ij}^{*}=\frac{\overline{u}_{ij}^{\overline{*}}+\overline{u}_{ij}^{\overline{\overline{*}}}}{2},$ $\overline{u}_{ij}^{**}=\frac{\overline{u}_{ij}^{\overline{**}}+\overline{u}_{ij}^{\overline{\overline{**}}}}{2}$, satisfy relations $(\ref{20a})$ and $(\ref{31})$, respectively. $u_{ij}^{*}$ and $u_{ij}^{**}$ are given by equations $(\ref{40})$ and $(\ref{41})$, respectively. The following equalities hold:

[TABLE]

and

[TABLE]

where the operators $\delta_{x}$ and $\delta_{y}$ are given by relation $(\ref{4}).$

Proof.

(of Lemma 3.1). Firstly, it is not hard to observe that

[TABLE]

for $i=0,1,...,M$ and $j=1,2,...,M-1.$ We should prove only equation $(\ref{72}).$ The proof of relation $(\ref{72a})$ is similar.

It follows from the definition of the operator $\delta_{x}^{2}$ and the scalar product $<\cdot,\cdot>_{x}$ given by $(\ref{4})$ and $(\ref{6}),$ respectively, that

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

It comes from the boundary condition $(\ref{24})$ that $e_{Mj}^{n}=e_{0j}^{n}=0.$ So $(\delta_{x}e_{M-\frac{1}{2},j}^{n})e_{Mj}^{n}=0$ and $(\delta_{x}e_{\frac{1}{2},j}^{n})e_{0j}^{n}=0.$ This fact, together with equation $(\ref{73})$ provide

[TABLE]

[TABLE]

∎

Proof.

(of Theorem 3.1). A combination of equations $(\ref{20a})$, $(\ref{35})$ and $(\ref{40})$ gives

[TABLE]

where $\rho_{ij}^{\alpha}$ is defined by $(\ref{20})$. Utilizing the definition of the operator $"\delta_{y}^{2}",$ equation $(\ref{45})$ is equivalent to

[TABLE]

Of course, the aim of this study is to give a general picture of the stability analysis of the numerical scheme $(\ref{40})$-$(\ref{24})$. Since the formulas can become quite heavy, for the sake of readability, we should neglect the higher order terms in time step $k$ and grid spacing $h.$ However, the truncation of the infinitesimal terms does not compromise the result on the stability analysis. Using this, equation $(\ref{46})$ becomes

[TABLE]

Taking the square, it holds

[TABLE]

[TABLE]

Now, using equality $2(a-b)b=a^{2}-b^{2}-(a-b)^{2}$ and inequality $(a\pm b)^{2}\leq 2(a^{2}+b^{2}),$ for any $a,b\in\mathbb{R},$ and by simple computations, it is not hard to see that

[TABLE]

[TABLE]

[TABLE]

$f\in\mathcal{C}^{1}(\mathbb{R})$ is a Lipschitz function, so there is a positive constant $C$ independent of the time step $k$ and the mesh size $h$ so that

[TABLE]

From inequality $(\ref{53})$, it is easy to see that

[TABLE]

A combination of estimates $(\ref{47})$-$(\ref{54})$ results in

[TABLE]

[TABLE]

Summing this up from $i,j=1,2,...,M-1,$ and rearranging terms, this provides

[TABLE]

[TABLE]

which implies

[TABLE]

[TABLE]

Multiplying both sides of inequality $(\ref{56})$ by $h^{2},$ and using equation $(\ref{72a})$ to get

[TABLE]

From the time step restriction $(\ref{44})$, $-1+\frac{2ak}{h^{2}}\leq 0,$ utilizing this, it follows

[TABLE]

In way similar, combining equations $(\ref{33}),$ $(\ref{35})$ and $(\ref{42})$, it is not hard to show that

[TABLE]

We must find a similar estimate associated with $\|e^{**}\|_{L^{2}(\Omega)}^{2}$ and $\|e^{*}\|_{L^{2}(\Omega)}^{2}.$ Using equations $(\ref{31}),$ $(\ref{35})$ and $(\ref{41})$, it holds

[TABLE]

for $i=1,2,...,M-1,$ and $j=0,1,...,M.$ Taking the square, we obtain

[TABLE]

which implies

[TABLE]

Now, summing relation $(\ref{59a})$ for $i,j=1,2,...,M-1,$ and multiplying the obtained equation by $h^{2},$ this results in

[TABLE]

which implies

[TABLE]

which is equivalent to

[TABLE]

Utilizing equality $(\ref{72}),$ this gives

[TABLE]

It comes from the time step restriction $(\ref{44})$ that $\frac{2ak}{h^{2}}\leq 1.$ So, estimate $(\ref{59b})$ provides

[TABLE]

Now, plugging inequalities $(\ref{57})$, $(\ref{59})$ and $(\ref{59c}),$ straightforward calculations yield

[TABLE]

Summing this up from $n=0,1,2,..,p-1,$ for any nonnegative integer $p$ satisfying $1\leq p\leq N,$ to get

[TABLE]

It comes from the initial condition given in $(\ref{24})$, that $e^{0}_{ij}=0,$ for $0\leq i,j\leq M.$ Applying the discrete Gronwall Lemma, estimate $(\ref{60})$ gives

[TABLE]

But $k=\frac{T}{N},$ so $Ckp=CT\frac{p}{N}\leq CT$ (since $p\leq N$). This fact, together with estimate $(\ref{61})$ result in

[TABLE]

Taking the square root, it is easy to see that

[TABLE]

We have that $\|u^{p}\|_{L^{2}(\Omega)}-\|\overline{u}^{p}\|_{L^{2}(\Omega)}\leq\|u^{p}-\overline{u}^{p}\|_{L^{2}(\Omega)}=\|e^{p}\|_{L^{2}(\Omega)}.$ A combination of this inequality together with relation $(\ref{62})$ yields

[TABLE]

[TABLE]

Since $\overline{u}$ is the exact solution, the proof of Theorem 3.1 is completed thanks to estimate $(\ref{43}).$ ∎

4 Convergence of the method

This section deals with the error estimates of a three-level time-split MacCormack method $(\ref{40})$-$(\ref{24})$ applied to equations $(\ref{1})$-$(\ref{3}),$ under the time step restriction $(\ref{44}).$ We assume that the exact solution $\overline{u}$ satisfies estimate $(\ref{43})$. We recall that

[TABLE]

is the space of grid functions defined on $\Omega_{h}\times\Omega_{k},$ where $\Omega_{k}=\{t^{n},\text{\,}0\leq n\leq N\}$ and

$\Omega_{h}=\{(x_{i},y_{j}),\text{\,}0\leq i,j\leq M\}\cap\Omega.$

Let introduce the following discrete norms

[TABLE]

and

[TABLE]

Theorem 4.1**.**

Suppose $u$ be the solution provided by a three-level time-split MacCormack approach $(\ref{40})$-$(\ref{24})$. Under the time step restriction $(\ref{44}),$ the error term $e=u-\overline{u},$ satisfies

[TABLE]

The proof of Theorem 4.1 requires some intermediate results (namely Lemmas 3.1, 4.1 and 4.2).

Lemma 4.1**.**

Consider $v\in H^{4}(\Omega),$ be a function satisfying $v|_{[x_{i},x_{i+1}]}\in\mathcal{C}^{6}[x_{i},x_{i+1}],$ for $i=0,1,2,...,M-1.$ Then, it holds

[TABLE]

where $\theta^{(4)}_{i}\in(x_{i-1},x_{i})$, $\theta^{(3)}_{i}\in(x_{i},x_{i+1})$ and $v_{mx}$ denotes the derivative of order $m$ of $v.$ Furthermore, for $i=2,3,...,M-2,$

[TABLE]

where $\theta^{(2)}_{i}\in(x_{i-2},x_{i-1})$, $\theta^{(4)}_{i}\in(x_{i-1},x_{i}),$ $\theta^{(3)}_{i}\in(x_{i},x_{i+1})$ and $\theta^{(1)}_{i}\in(x_{i+1},x_{i+2}).$

Proof.

(of Lemma 4.1) Expanding the Taylor series about $x_{i}$ with grid spacing $h$ using both forward and backward differences to obtain

[TABLE]

where $\theta^{(1)}_{i}\in(x_{i+1},x_{i+2});$

[TABLE]

where $\theta^{(2)}_{i}\in(x_{i-2},x_{i-1});$

[TABLE]

where $\theta^{(3)}_{i}\in(x_{i},x_{i+1});$

[TABLE]

where $\theta^{(4)}_{i}\in(x_{i-1},x_{i}).$

In way similar, applying the Taylor expansion for both derivative and higher order derivatives of $v$ to obtain

[TABLE]

where $\theta^{(5)}_{i}\in(x_{i},x_{i+1});$

[TABLE]

where $\theta^{(6)}_{i}\in(x_{i-1},x_{i});$

[TABLE]

where $\theta^{(7)}_{i}\in(x_{i},x_{i+1});$

[TABLE]

where $\theta^{(8)}_{i}\in(x_{i-1},x_{i});$

[TABLE]

where $\theta^{(9)}_{i}\in(x_{i},x_{i+1}),$ $\theta^{(10)}_{i}\in(x_{i-1},x_{i});$

[TABLE]

where $\theta^{(11)}_{i}\in(x_{i},x_{i+1}),$ $\theta^{(12)}_{i}\in(x_{i-1},x_{i});$

[TABLE]

where $\theta^{(13)}_{i}\in(x_{i},x_{i+1}),$ $\theta^{(14)}_{i}\in(x_{i-1},x_{i}).$

Now, adding equations $(\ref{78})$-$(\ref{79})$ side by side, this gives

[TABLE]

[TABLE]

Subtracting $(\ref{105})$ from $(\ref{104})$ and adding side by side $(\ref{106})$ and $(\ref{107}),$ using also equations $(\ref{108}),$ $(\ref{110})$ and $(\ref{112}),$ simple calculations provide

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Combining equations $(\ref{104})$-$(\ref{117})$, straightforward computations result in

[TABLE]

[TABLE]

[TABLE]

Since $\theta^{(3)}_{i},\theta^{(5)}_{i},\theta^{(7)}_{i},\theta^{(9)}_{i},\theta^{(11)}_{i},\theta^{(13)}_{i}\in(x_{i},x_{i+1})$ and $\theta^{(4)}_{i},\theta^{(6)}_{i},\theta^{(8)}_{i},\theta^{(10)}_{i},\theta^{(12)}_{i},\theta^{(14)}_{i}\in(x_{i-1},x_{i}),$ without loss of generality, we can assume that $\theta^{(3)}_{i}=\theta^{(5)}_{i}=\theta^{(7)}_{i}=\theta^{(9)}_{i}=\theta^{(11)}_{i}=\theta^{(13)}_{i}$ and $\theta^{(4)}_{i}=\theta^{(6)}_{i}=\theta^{(8)}_{i}=\theta^{(10)}_{i}=\theta^{(12)}_{i}=\theta^{(14)}_{i}.$ Using this, relation $(\ref{112})$ becomes

[TABLE]

This completes the proof of the first item of Lemma $\ref{l2}.$ Now, let prove the second item of Lemma $\ref{l2}.$

Plugging equations $(\ref{76})$ and $(\ref{78}),$ $(\ref{77})$ and $(\ref{79}),$ $(\ref{78})$ and $(\ref{79}),$ respectively, it is not hard to see that

[TABLE]

[TABLE]

and

[TABLE]

[TABLE]

A combination of equations $(\ref{100})$-$(\ref{101})$ yields

[TABLE]

[TABLE]

Substituting $(\ref{114})$-$(\ref{117})$ into $(\ref{103})$,simple computations result in

[TABLE]

[TABLE]

[TABLE]

which is equivalent to

[TABLE]

[TABLE]

[TABLE]

Assuming that $\theta^{(3)}_{i}=\theta^{(5)}_{i}=\theta^{(7)}_{i}=\theta^{(9)}_{i}=\theta^{(11)}_{i}=\theta^{(13)}_{i}$ and $\theta^{(4)}_{i}=\theta^{(6)}_{i}=\theta^{(8)}_{i}=\theta^{(10)}_{i}=\theta^{(12)}_{i}=\theta^{(14)}_{i},$ equation $(\ref{119})$ becomes

[TABLE]

which is equivalent to

[TABLE]

This ends the proof of Lemma $\ref{l2}.$ ∎

Lemma 4.2**.**

The term $\rho_{ij}^{n}$ given by equation $(\ref{20})$ can be bounded as

[TABLE]

where $\widehat{C}_{l},$ $l=1,2,3,$ are positive constant independent of the time step $k$ and the mesh size $h.$

Proof.

(of Lemma 4.2). It comes from relation $(\ref{20})$ that

[TABLE]

From the definition of the operator $"\delta_{y}^{2}",$ this is equivalent to

[TABLE]

[TABLE]

Combining this together with Lemma $\ref{l2},$ it is easy to see that

[TABLE]

[TABLE]

[TABLE]

On the other hand, $\overline{u}(x,\cdot,t)|_{[y_{j},y_{j+1}]}$, $fo\overline{u}(x,\cdot,t)|_{[y_{j},y_{j+1}]}\in\mathcal{C}^{6}([y_{j},y_{j+1}])$, for every $x\in(0,1)$, $t\in(0,T)$ and $j=0,1,2,...,M-1;$ $\||\overline{u}|\|_{L^{\infty}(0,T;L^{2}(\Omega))}\leq\widetilde{C}$ (according to estimate $(\ref{43})$) and $f^{{}^{\prime}}$ (the derivative of $f$) is continuous. Taking the absolute value of $\rho_{ij}^{n},$ there exist positive constants $\widehat{C}_{l},$ $l=1,2,3,$ independent of the time step $k$ and the mesh grid $h$ so that

[TABLE]

This completes the proof of Lemma $\ref{l3}.$ ∎

Armed with the results provided by Lemmas 3.1, 4.1 and $\ref{l3},$ we are ready to prove Theorem $\ref{t2}.$

Proof.

(of Theorem 4.1) We recall that the error term provided by the scheme $(\ref{40})$-$(\ref{24})$ is denoted by $e_{ij}^{n}=u_{ij}^{n}-\overline{u}_{ij}^{n},$ where $\overline{u}$ satisfies equations $(\ref{20a}),$ $(\ref{31})$ and $(\ref{33})$ and $u$ is given by relations $(\ref{40})$-$(\ref{24})$. So, it comes from equation $(\ref{46})$ that

[TABLE]

which is equivalent to

[TABLE]

where $C_{r}$ is a parameter that does not depend neither on the time step $k$ nor the grid spacing $h$ and $\rho_{ij}^{\alpha}$ is defined by $(\ref{20})$. Taking the square, it is not hard to see that

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Applying the inequalities: $2ab\leq a^{2}+b^{2},$ $(a\pm b)^{2}\leq 2(a^{2}+b^{2})$ and $(a\pm b\pm c)^{2}\leq 3(a^{2}+b^{2}+c^{2}),$ for every $a,b,c\in\mathbb{R},$ together with the time step restriction $(\ref{44})$ (that is, $2ak\leq h^{2}$), relation $(\ref{67})$ becomes

[TABLE]

[TABLE]

[TABLE]

[TABLE]

which implies

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

From estimates $(\ref{53})$-$(\ref{54})$, we have that

[TABLE]

This fact, together with estimate $(\ref{68})$ results in

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Utilizing the time step restriction $(\ref{44})$, $\frac{2ak}{h^{2}}\leq 1,$ this implies

[TABLE]

[TABLE]

[TABLE]

Summing this up from $i,j=1,2,...M-1,$ provides

[TABLE]

[TABLE]

[TABLE]

Combining the boundary condition $(\ref{24})$, $e^{n}_{Mj}=e^{n}_{0j}=0,$ for all $j=0,1,...,M,$ Lemmas 3.1 and $\ref{l3},$ and multiplying both sides of inequality $(\ref{56a})$ by $h^{2},$ straightforward computations yield

[TABLE]

[TABLE]

[TABLE]

Since $h=\frac{1}{M},$ $k\leq 1+k^{2}$ and $h^{2}\leq 1+h^{4},$ this becomes

[TABLE]

[TABLE]

which implies

[TABLE]

[TABLE]

where we absorbed all the constants into a constant $\widehat{C}_{4}.$

Similarly, one shows that

[TABLE]

[TABLE]

where all the constants have been absorbed into a constant $\widehat{C}_{5},$ and

[TABLE]

[TABLE]

where all the constants have been absorbed into a constant $\widehat{C}_{6}.$

Now, setting

[TABLE]

and

[TABLE]

plugging estimates $(\ref{121})$-$(\ref{123}),$ straightforward calculations yield

[TABLE]

[TABLE]

Absorbing all the constants into a constant $\widehat{C}_{7},$ this yields

[TABLE]

[TABLE]

Summing this up from $n=0,1,2,..,p-1,$ for any nonnegative integer $p$ such that $1\leq p\leq N,$ we obtain

[TABLE]

[TABLE]

It comes from the initial condition given in $(\ref{24})$, that $e^{0}_{ij}=0,$ for $0\leq i,j\leq M.$ Applying the Gronwall Lemma, estimate $(\ref{126})$ provides

[TABLE]

But $k=\frac{T}{N},$ so $\widehat{C}_{7}kp=\widehat{C}_{7}T\frac{p}{N}\leq\widehat{C}_{7}T$ (since $p\leq N$). This fact, together with inequality $(\ref{127})$ result in

[TABLE]

where $\varphi_{3}(k,h)^{2}=k^{-1}\varphi_{2}(k,h),$ $\varphi_{2}(k,h)$ is given by equation $(\ref{125}).$ Taking the square root, it is easy to see that

[TABLE]

It comes from equality $\varphi_{3}(k,h)^{2}=k^{-1}\varphi_{2}(k,h),$ and equation $(\ref{125})$ that

[TABLE]

[TABLE]

where $\widetilde{C}_{8}$ is a positive constant independent of $k$ and $h,$ and $\varphi_{4}(k,h)$ tends to zero when $k,h\rightarrow 0.$ Taking the maximum over $p$ of estimate $(\ref{128})$, for $0\leq p\leq N$, the proof of Theorem 4.1 is completed thanks to equation $(\ref{66})$. ∎

5 Numerical experiments and Convergence rate

In this section we construct an exact solution to the initial-boundary value problem $(\ref{1})$-$(\ref{3})$ for a specific source term $f$. Furthermore, using Matlab we perform some numerical experiments in bidimensional case. In that case we obtain satisfactory results, so our algorithm performances are not worse for multidimensional problems. We consider two cases which are physical examples associated with the diffusive coefficient $a=1,$ together with the example introduced in [8]. We confirm the predicted convergence rate from the theory (see Section 2, Page $6$, last paragraph). This convergence rate is obtained by listing in Tables $1$-$6$ the errors between the computed solution and the exact one with different values of mesh size $h$ and time step $k,$ satisfying $k=\frac{1}{2}h^{2}.$ Finally, we look at the error estimates of our proposed method for the parameter $T=1.$

Assuming that the exact solution to problem $(\ref{1})$-$(\ref{3})$ is of the form $\overline{u}(x,y,t)=\left[1+\exp(ct+dx+by)\right]^{-n},$ where $n$ is an integer. By simple calculations, it holds

[TABLE]

[TABLE]

and

[TABLE]

In way similar

[TABLE]

Combining equations $(\ref{1n})$-$(\ref{3n}),$ it is not hard to see that

[TABLE]

[TABLE]

Setting $c=-(d^{2}+b^{2}),$ this becomes

[TABLE]

[TABLE]

[TABLE]

$\bullet$: Case 1: $n=1$. In this case, equation $(\ref{4n})$ becomes

[TABLE]

Now, taking $2(d^{2}+b^{2})=1,$ this gives $b^{2}=\frac{1}{2}-d^{2}.$ Since $b^{2}$ must be strictly greater than zero, this implies $d^{2}<\frac{1}{2}.$ For $d=\pm\frac{\sqrt{3}}{3},$ this implies $b=\pm\frac{\sqrt{6}}{6}$ and $c=-\frac{1}{2}.$ Letting $f(\overline{u})=(1-\overline{u})\overline{u}^{2},$ our exact solution is given by $\overline{u}(x,y,t)=\left[1+\exp\left(-\frac{1}{2}t+\frac{\sqrt{3}}{3}x+\frac{\sqrt{6}}{6}y\right)\right]^{-1},$ for $t\in[0,1]$ and $(x,y)\in[0,1]^{2}.$ The initial and boundary conditions are determined by this solution.

$\bullet$: Case 2: $n=-1$. It comes from equation $(\ref{4n})$ that

[TABLE]

Since $\overline{u}=1+\exp(ct+dx+by),$ so $-\exp(ct+dx+by)=1-u.$ Taking $-2(d^{2}+b^{2})=-1,$ it holds $d=\pm\frac{\sqrt{3}}{3},$ $b=\pm\frac{\sqrt{6}}{6}$ and $c=-\frac{1}{2}.$ Setting $f(\overline{u})=1-\overline{u},$ we consider the exact solution defined as $\overline{u}(x,y,t)=1+\exp\left(-\frac{1}{2}t+\frac{\sqrt{3}}{3}x+\frac{\sqrt{6}}{6}y\right),$ for $t\in[0,1]$ and $(x,y)\in[0,1]^{2}.$ The initial and boundary conditions are determined by this solution.

To analyze the convergence rate of our numerical scheme, we take the mesh size $h\in\{\frac{1}{2},\frac{1}{2^{2}},\frac{1}{2^{3}},\frac{1}{2^{4}},\frac{1}{2^{5}}\}$ and time step $k\in\{\frac{1}{2^{2}},\frac{1}{2^{3}},\frac{1}{2^{4}},\frac{1}{2^{5}},\frac{1}{2^{6}},\frac{1}{2^{7}},\frac{1}{2^{8}},\frac{1}{2^{9}},\frac{1}{2^{10}}\frac{1}{2^{11}}\},$ by a mid-point refinement. Under the time step restriction $(\ref{44}),$ we set $k=\frac{1}{2}h^{2}$ and we compute the error estimates: $\||E(u)|\|_{L^{2}(0,T;L^{2})},$ $\||E(u)|\|_{L^{\infty}(0,T;L^{2})}$ and $\||E(u)|\|_{L^{1}(0,T;L^{2})}$ related to the time-split method to see that the algorithm is stable, second order accuracy in time and fourth order convergent in space. In addition, we plot the approximate solution, the exact one and the errors versus $n.$ From this analysis, a three-level explicit time-split MaCormack method is both efficient and effective than a two-level linearized compact ADI approach. In fact, although the two-level linearized compact ADI scheme has the same convergent rate (see [37], Theorem $6.6,$ p. $19$) this method requires too much computer times to achieve the solution. Furthermore, when $h$ varies in the given range, we observe from Tables $1$-$6$ that the approximation errors $O(k^{\beta})+O(h^{\theta})$ are dominated by the h-terms $O(h^{\theta})$ (or $k$-terms $O(k^{\beta})$). So, the ratio $r^{m}_{u},$ where $m=1,2,\infty,$ of the approximation errors on two adjacent mesh levels $\Omega_{2h}$ and $\Omega_{h}$ is approximately $(2h)^{\theta}/h^{\theta}=2^{\theta},$ where $m$ refers to the $L^{m}(0,T;L^{2}(\Omega)$-error norm. Hence, we can simply use $r^{m}_{u}$ to estimate the corresponding convergence rate with respect to $h.$ Define the norms for the approximate solution $u,$ the exact one $\overline{u},$ and the errors $E(u),$ as follows

[TABLE]

[TABLE]

and

[TABLE]

$\bullet$ Test $1.$ Let $\Omega$ be the unit square $(0,1)\times(0,1)$ and $T$ be the final time, $T=1.$ We assume that the diffusive coefficient $a=1,$ and we choose the force $f(\overline{u})=(1-\overline{u})\overline{u}^{2},$ in such a way that the exact solution $\overline{u}$ is given by

[TABLE]

The initial and boundary conditions are given by this solution. We take the mesh size and time step: $h\in\{\frac{1}{2},\frac{1}{2^{2}},\frac{1}{2^{3}},\frac{1}{2^{4}},\frac{1}{2^{5}}\}$ and $k\in\{\frac{1}{2^{2}},\frac{1}{2^{3}},\frac{1}{2^{4}},\frac{1}{2^{5}},\frac{1}{2^{6}},\frac{1}{2^{7}},\frac{1}{2^{8}},\frac{1}{2^{9}},\frac{1}{2^{10}},\frac{1}{2^{11}}\}.$

Tables 1,2. Analyzing of convergence rate $O(h^{\theta}+\Delta t^{\beta})$ for time-split MacCormack by $r^{m}_{u},$ with varying time step $k=\Delta t$ and mesh grid $h=\Delta x$.

Case: $k=\frac{1}{2}h^{2}$.

[TABLE]

** Case: $k=h^{2}$.**

[TABLE]

$\bullet$ Test $2.$ Now, let $\Omega$ be the unit square $(0,1)^{2}$ and $T=1.$ The diffusive term $a$ is assumed equals $1.$ We choose the force $f$ such that the analytic solution $\overline{u}$ is defined as

[TABLE]

The initial and boundary conditions also are given by the exact solution $\overline{u}.$ Similar to Test $1,$ we take the mesh size and time step: $h\in\{\frac{1}{2},\frac{1}{2^{2}},\frac{1}{2^{3}},\frac{1}{2^{4}},\frac{1}{2^{5}}\}$ and $k\in\{\frac{1}{2^{2}},\frac{1}{2^{3}},\frac{1}{2^{4}},\frac{1}{2^{5}},\frac{1}{2^{6}},\frac{1}{2^{7}},\frac{1}{2^{8}},\frac{1}{2^{9}},\frac{1}{2^{10}},\frac{1}{2^{11}}\}.$

Tables 3,4. Convergence rates $O(h^{\theta}+\Delta t^{\beta})$ for time-split MacCormack by $r^{m}_{u},$ with varying spacing $h$ and time step $k$.

Case: $k=\frac{1}{2}h^{2}$.

[TABLE]

Case: $k=h^{2}$.

[TABLE]

$\bullet$ Test $3.$ Finally, let $\Omega$ be the unit square $(0,1)\times(0,1)$ and $T=1.$ We assume that $a=1,$ and the force $f$ is chosen such that the exact solution $u$ is given by

[TABLE]

The initial and boundary conditions are given by the exact solution $\overline{u}$.

Similar to both Tests $1,2$ the mesh size and time step are chosen such that: $h\in\{\frac{1}{2},\frac{1}{2^{2}},\frac{1}{2^{3}},\frac{1}{2^{4}},\frac{1}{2^{5}}\}$ and $k\in\{\frac{1}{2^{2}},\frac{1}{2^{3}},\frac{1}{2^{4}},\frac{1}{2^{5}},\frac{1}{2^{6}},\frac{1}{2^{7}},\frac{1}{2^{8}},\frac{1}{2^{9}},\frac{1}{2^{10}},\frac{1}{2^{11}}\},$ by a mid-point refinement. We compute the error estimates: $E(u)$ related to a three-level explicit time-split MacCormack approach to see that the algorithm is second order convergent in time and fourth order accurate in space. Furthermore, we plot the errors together with the energies versus $n.$ From this analysis, it is obvious that a three-level time-split scheme is efficient and effective than a two-level linearized compact ADI method which has the same order of convergence.

Tables 5,6. Convergence rates $O(h^{\theta}+\Delta t^{\beta})$ for time-split MacCormack by $r^{m}_{u},$ with varying spacing $h$ and time step $\Delta t$.

Case: $k=\frac{1}{2}h^{2}$.

[TABLE]

Case: $k=h^{2}$.

[TABLE]

The analysis on the convergence of the numerical scheme presented in Section $\ref{sec4},$ has suggested that our algorithm is first order convergent in time and fourth order accurate in space. If the result provided in Section $\ref{sec2},$ page $6$, last paragraph is to believe, this shows that the time-split MacCormack scheme is inconsistent. Surprisingly, it comes from Tests 1-3, more precisely Figures 1-3 and Tables 1-6, that the three-level explicit time-split MacCormack technique is stable, second order accurate in time and fourth order convergent in space under the time step restriction $(\ref{44}),$ which confirms the theoretical result provided in Section $\ref{sec2},$ page $6$, last paragraph. Thus, the considered method applied to initial-boundary value problem $(\ref{1})$-$(\ref{3})$ is: stable, consistent, second order convergent in time and fourth order accurate in space.

6 General conclusion and future works

We have studied in detail the stability, error estimates and convergence rate of a three-level explicit time-split MacCormack method for solving the $2$D nonlinear reaction-diffusion equation $(\ref{1})$-$(\ref{3})$. The analysis has suggested that our method is stable, consistent, second order accuracy in time and fourth order convergent in space under the time step restriction $(\ref{44})$. This convergence rate is confirmed by a large set of numerical experiments (see both Figures 1-3 and Tables 1-6). Numerical evidences also show that the new algorithm is: (1) more efficient and effective than a two-level linearized compact ADI method, (2) fast and robust tools for the integration of general systems of parabolic PDEs. However, the time-split MacCormack method is not is a satisfactory approach for solving high Reynolds number flows where the viscous region becomes very thin. For these flows, the mesh grid must be highly refined in order to accurately resolve the viscous regions. This leads to very small time steps and subsequently long computing times. To overcome this difficulty, MacCormack developed a hybrid version of his scheme, which is known as MacCormack rapid solver method [18]. This hybrid scheme is an explicit-implicit method which has been proved to be from $10$ to $100$ more faster than a time-split MacCormack algorithm (see [2], P. 632). The rapid solver method will be applied to the two-dimensional nonlinear reaction-diffusion equations in our future works.