A robust adaptive grid method for singularly perturbed Burger-Huxley equations

1.   Introduction

In this paper, we consider the following one-dimensional unsteady singularly perturbed Burger-Huxley equation

where $0<\varepsilon\ll 1$ is a perturbation parameter and $\alpha,\;\beta\geq 0$, $\gamma\in (0,1)$ are given constants. Such equation describes the interaction between convection, diffusion and reaction.

It is well known that the Burger-Huxley equation describes many phenomena such as busting oscillation[9], interspike[24], population genetics[1], bifurcation and chaos[33]. In the past few decades, various analytical methods were proposed to solve the Burger-Huxley equation. For instance, by using Hirota method, Satsuma [30] obtained an exact solitary solution for this equation. Wang et.al., [31] constructed an exact solitary wave solution of the generalized Burger-Huxley equation. In[32], Wazwaz derived some travelling wave solutions for the generalized forms of Burgers, Burgers-KdV and Burger-Huxley equation by using the standard tanh-function technique. Recently, many researches paid attention to solve Burger-Huxley equation by the numerical methods such as a domian decomposition method[15,14,16], variational iteration method[2], finite difference scheme[27,28], spectral methods[8,17,18], collocation method[26,20,25]. However, for suitable values of $\alpha$, $\beta$, $\gamma$ and $0<\varepsilon\ll 1$, the solution of Burger-Huxley problem, which is called singularly perturbed Burger-Huxley problem, usually has one or two boundary layers. Due to the presence of boundary layer(s), the above presented methods are in question and known to be inadequate to approximate the exact solution. Therefore, to obtain $\varepsilon$-uniformly convergent methods, Kaushik and Sharma[19] developed a parameter-uniform convergent finite difference scheme for the problem (1). Gupta and Kadalbajoo[13] constructed a numerical scheme that comprises of implicit-Euler method which is first order uniformly accurate to discretize in temporal direction on a uniform mesh and a monotone hybrid finite difference operator to discretize the spatial variable with a piecewise uniform Shishkin mesh. For the spatial direction, their method has been shown to be first order parameter uniform convergent in the outside region of boundary layer, and almost second order parameter uniform convergent in the boundary layer region.

Obviously, the methods presented in [19,13] are belong to the lay-adapted mesh approach. As far as we know, the convergence results for the layer-adapted mesh approach is more easy to be obtained. But this special mesh approach requires a priori information about the location and width of the boundary layer. So, it is very necessary to study the adaptive moving grid approach by equidistributing a positive monitor function. Since this method has an advantage that it can be applied using little or no a priori information.

In order to serve this purpose, we will devise a uniformly convergent numerical scheme for the problem (1). We first use the classical backward-Euler formula on a uniform mesh to approximate the time derivative. The resulting nonlinear singularly perturbed ordinary differential equation is linearized by using quasi-lineariztion technique. Then, an upwind finite difference scheme is used for space discretization. The a posteriori error bound is derived to design an adaptive spatial grid generation algorithm. Finally, the convergence analysis is given to show that the presented adaptive grid method is first-order in the time and spatial directions, respectively.

2.   Bounds on the solution and its time derivatives

In this section, we will give error estimates for the solution of continuous problem and its time derivatives. Let $D = (0,1)\times (0,T]$, the initial boundary $\Gamma_i = \{(x,t):t = 0,x\in [0,1]\}$, left boundary $\Gamma_l = \{(x,t):x = 0,t\in [0,T]\}$ and right boundary $\Gamma_r = \{(x,t):x = 1,t\in [0,T]\}$, then, $\partial D = \Gamma_l\bigcup \Gamma_i\bigcup \Gamma_r$. For any given function $g(x,t)\in C^0(\overline{D})$, the maximum norm is defined as follows

First, the operator $\mathcal{L}_{x,\varepsilon}$ defined in (1) satisfies the following maximum principle.

Lemma 2.1. (Maximum principle) Let $v(x,t)\in C^{2,1}(\overline{D})$. If $v(x,t)\geq 0$, $\forall (x,t)\in \partial D$ and $\mathcal{L}_{x,\varepsilon} v(x,t)\geq 0$, $\forall (x,t)\in D$, then $v(x,t)\geq 0$, $\forall (x,t)\in \overline{D}$.

Proof. The proof can be seen in Lemma 1 of [13].

Furthermore, based on the above maximum principle, the following lemma [13] gives the uniform stability estimate.

Lemma 2.2. Let $u(x,t)$ be the solution of the problem (1), then we have

Finally, to derive the convergence and stability of the time discrete scheme, we give the bounds on the time derivatives as follows:

Lemma 2.3. Let $u(x,t)$ be the solution of the problem (1). Then there exists a constant $C$, independent of $\varepsilon$, such that

Proof. It follows from Lemma 2.2, we can obtain $\left|u(x,t)\right|\leq C$. Next, one can obtain the first-order derivative bound with respect to time variable as follows.

First, at $t = 0$, we have $u(x,0) = u_0(x)$, which gives $\frac{\partial u(x,0)}{\partial x} = u'_0(x)$ and $\frac{\partial^2 u(x,0)}{\partial x^2} = u''_0(x)$, $x\in [0,1]$. Then, we get

On the boundaries $x = 0$ and $x = 1$, we obtain

Therefore, for sufficiently smooth functions $u_0(x)$, $S_0(t)$ and $S_1(t)$, there exists a constant $C_1$, such that

For $\forall (x,t)\in (0,1)\times (0,T]$, it follows from (1) that

As the operator $\mathcal{L}_{x,\varepsilon}$ is uniformly stable, using Lemma 2.2, we get

Similarly, we get the bound for the second-order derivative $\frac{\partial^2 u(x,t)}{\partial t^2}$.

3.   Time semi-discretization and quasi-linearization

3.1.   Time semi-discretization

Let

be a uniform mesh in time direction, where $M$ denotes the number of mesh intervals in the time direction. Here, discretizing the problem (1) with respect to temporal variable by using the backward Euler method, we can obtain the following non-linear ordinary differential equations

where $u^n$ is the semi-discrete approximation to the exact solution $u(x,t)$ of the continuous problem (1) at time level $t_n = n \Delta t$, $0\leq n\leq M-1$.

Next, to analyze the uniform convergence of the solution $u^n$ to the exact solution $u(x,t_n)$, we will do the stability analysis and derive the consistency result of the scheme (2). Obviously, the operator $(I+\Delta t \bar{\mathcal{L}}_{x,\varepsilon})$ satisfies the maximum principle, which ensures the stability of the scheme (2).

Let $e_{n+1} = u(x,t_{n+1})-\widehat{u}^{n+1}(x)$ be the local error of semi-discretization scheme (2), where $\widehat{u}^{n+1}(x)$ is the solution obtained after one step of semi-discrete scheme (2) by taking the exact value $u(x,t_n)$ instead of $u^n(x)$ as the starting data. Then, we obtain the following boundary value problem

Lemma 3.1. Let $e_{n+1}$ be the local error of semi-discretization scheme (2), then, we have

Proof. Since the solution of the problem (1) is smooth enough, by using Taylor series expansion, we have

where $\xi\in (t^n,t^{n+1})$.

Then by applying Taylor expansion to linearize the nonlinear differential operator $\bar{\mathcal{L}}_{x,\varepsilon}$, the desired estimate (4) follows from Lemma 1 of [7].

Based on the above results, we can obtain the following global error estimate.

Theorem 3.2.

Thus, the semi-discretization (2) is uniformly convergent of first-order in time.

3.2.   Quasilinearization

In this section, we will introduce the Newtom-Raphson-Kantorovich approximation approach to linearize the nonlinear differential equation (2). For the nonlinear ordinary differential equation $u'' = f(u,u',x)$, Bellman and Kalaba[3] developed the following linearization scheme

where $\left\{u_{(k)}\right\}_{k = 0}^{\infty}$ are a sequence of approximate solutions of $u'' = f(u,u',x)$.

Here, in our case, based on the scheme (6), the time semi-discretisation (2) can be written as

with the boundary conditions

where $u_{(k)}$ is the nominal solution of the problem (2) with a reasonable initial guess $u_{(0)}$ satisfying the initial condition $u_0(x)$ and $k = 0,1,2,\cdots,$ is the iteration index.

For the sake of convenience, let $u_{(k+1)} = \widetilde{u}$, the above equation can be written as

where

Here, for $n\geq 0$, assuming that $a^{n+1}(x)$, $b^{n+1}(x)$ and $f^{n+1}(x)$ are sufficiently smooth functions and there exist two positive constants $p$ and $q$ such that

Under these conditions, the solution of above problem (10)-(12) has a unique solution that usually has a boundary layer at $x = 1$ for $\varepsilon\rightarrow 0$, see, e.g., [29]. For $k = 0,1,2,\cdots$ and $n\geq 0$, we will solve the sequence of second order singularly perturbed linear boundary value problem (10)-(12) instead of solving the original nonlinear problem (2). Furthermore, for the solution $u^{n+1}(x)$ of original nonlinear problem (2), we should prove that

whereas numerically, we require that

where $\varepsilon$ is a small tolerance value to terminate the computation.

Now, to obtain (13) or (14), the following theorem shows that not only the convergence of this sequence $\left\{u_{(k)}^{n+1}\right\}_{k = 0}^{\infty}$ is quadratic, but also its proportionality constant is independent of $k$.

Theorem 3.3. Let $\left\{u_{(k)}^{n+1}\right\}_{k = 0}^{\infty}$ be the sequence produced by quasilinearization technique at $(n+1)$ time level. Then, there exists a positive constant $C$, independent of $k$, such that

Proof. The proof can be seen in Theorem 5 of [13].

4.   Spatial discretization

In this section, we shall develop a finite difference scheme of (10)-(12) on the following non-uniform spatial mesh

Here, the spatial step size is denoted by

Firstly, for a given mesh function $v(x_i,t_n) = v_i^n$, we define some difference operators as follows:

Then, on $\overline{\Omega}_x^N$, the finite difference spatial discretization of (10)-(12) takes the form

where $\bar{\mathcal{L}}_{x,\varepsilon}^N$ is the discretization of the differential operator $\bar{\mathcal{L}}_{x,\varepsilon}$, which takes the following form:

Here, $U_i^n$ is the numerical solution of exact solution $u(x_i,t_n)$ and

Remark 1. Here, by choosing a reasonable initial approximation $u_{(0)}^{n+1}$, we iteratively solve the linear boundary value problem (10)-(12) to obtain the solution $u_{(k)}^{n+1}$. At every iteration, we opt to numerically solve the problem (10)-(12) using the upwind finite difference scheme (15), where we use the center finite difference scheme to approximate the derivatives $\frac{\partial u_{(k)}^{n+1}}{\partial x}$ given in $b^{n+1}(x)$ and $f^{n+1}(x)$.

5.   Adaptive spatial grid algorithm

It is well known that the monitor functions are used by many authors [21,22,4,5] to obtain adaptive grid algorithms that produce layer-resolving meshes in solving singularly perturbed ordinary differential equations. In recent years, there has been tremendous interest in developing the adaptive grid approach of singularly perturbed parabolic equations [10,11,12]. It should be pointed out that the authors in [10,11] obtained a different spatial mesh on each time level $t_n$ by equidistributing a monitor function $M(u(x,t_n),x)$. Meanwhile, the authors[12] used the idea of equidistributing the spatial grids for a fixed time level $T_0$($0<T_0\leq T$) and obtained a spatial grid for all the levels. However, it is hard to find the fixed time level $T_0$ which can reflect the information of boundary layer. This drawback motivates to construct a spatial grid which is suitable for all time levels, in which the same order of convergence and the error estimation can be established.

Here, for singularly perturbed Burger-Huxley equation, in order to construct such a nonuniform spatial grid, we use the idea of equidistribution of a monitor function which is given by

Thus, a grid is said to be equidistributing $M(u(x,t),x)$, if

Let $\widehat{U}^{n}(x)\in C[0,1]$ be a piecewise linear interpolant function through the knots $(x_i,U_i^n)$. Then, in practical computation, we choose the following monitor function

which is the discrete analogue of the monitor function (16).

Therefore, to obtain the equidistribution grid and the corresponding numerical solution, we construct the following iteration algorithm:

Step 1. Take $k = 0$, let the uniform spatial mesh $\left\{x_i^{(k)} = i/N,i = 0,1,\cdots,N\right\}$ be the initial mesh.

Step 2. For $k = 1,2,\cdots$, assuming that the mesh $\left\{x_i^{(k)}\right\}$ is given, compute the discrete solution $U_i^{n,(k)}$ that satisfy (15) on $\left\{x_i^{(k)}\right\}$. For each $i$, let $h_i^{(k)} = x_i^{(k)}-x_{i-1}^{(k)}$ and

be the maximum distance between the points $\left(x_{i-1}^{(k)}, U_{i-1}^{n,(k)}\right)$ and $\left(x_i^{(k)}, U_i^{n,(k)}\right)$. Then, the total maximum arc-length can be written as

Step 3. Test mesh: Let $C_0$ be a user chosen constant with $C_0> 1$. If the stopping criterion

holds true, then go to Step 5. Otherwise, continue to Step 4.

Step 4. Generate a new mesh by equidistributing arc-length of current computed solution: choose point $0 = x_0^{(k+1)}<x_1^{(k+1)}<\cdots<x_N^{(k+1)} = 1$ such that

where $i = 1,2,\cdots,N$ and $\widehat{U}^{n,(k)}(x)\in C[0,1]$ is a piecewise linear interpolant function through the knots $\left(x_i,U_i^{n,(k)}\right)$. Return to Step 3.

Step 5. Take $\left\{x_{i}^{(k)}\right\}$ as the final mesh and calculate the corresponding solution $U_i^n$, then stop.

6.   Error analysis

In this section, in order to derive an a posteriori error estimate for the fully discrete scheme of (15), we will list some preliminary results as follows.

Firstly, we rewrite (11) in the form

where $\phi^{n+1}(x) = f^{n+1}(x) = f_{(k)}(x,t^{n+1})$, $n\geq 0$.

It is easy to see that the operator $(I+\Delta t\widetilde{\mathcal{L}}_{x,\varepsilon})$ satisfies the maximum principle, which implies $\left\|\widetilde{u}^{n+1}\right\|_{\infty}\leq C$. In addition, under the sufficient smoothness of the function $f^{n+1}(x)$, the function $g^n(x)+f^{n+1}(x)$ is $\varepsilon$-uniformly bounded in the spatial domain.

Denote

the Equation (20) can be written into the following system of equations

where $A(x) = \text{diag}(a^1(x),\cdots,a^M(x))$, $B(x) = \text{diag}(b^1(x),\cdots,b^M(x))$ are $M\times M$ diagonal matrixes, respectively.

Obviously, the simple upwind finite difference scheme of (22) can be written into the following matrix form

where $\widetilde{\bf U}_i = (\widetilde{U}_i^1,\cdots,\widetilde{U}_i^M)^T$ is the approximation solution of $\widetilde{\bf {u}}(x_i)$.

Next, based on the above preparation work, we can obtain the following result:

Theorem 6.1. Let $\widetilde{\bf {u}}(x)$ be the solution of (22), $\widetilde{\bf {U}}_i$ be the solution of (24) on an arbitrary nonuniform mesh $\overline{\Omega}_x^N$ and ${\widehat{{\bf U}}}(x)$ be piecewise linear interpolant function vector through knots $(x_i,\widetilde{\bf {U}}_i)$. Then we have

Proof. The proof is similar to Theorem 3.1 in [23].

Furthermore, based on the properties of discrete Green's function[22], we can obtain the following convergence result:

Theorem 6.2. Let $\widetilde{\bf {u}}(x)$ be the solution of (22), $\widetilde{\bf {U}}_i$ be the solution of (24) on a nonuniform grid $\left\{x_i\right\}_{i = 0}^N$ generated by the above adaptive spatial grid algorithm. Then we have the following bound

Proof. The proof is similar to Theorem 4.2 in [6].

From Theorem 6.2, an easy induction gives the following result:

Theorem 6.3. For each time level $n = 1,\cdots,M$, let $\widetilde{u}^{n}(x)$ be the solution of problem (10)-(12) and $\widetilde{U}_i^{n}$ be the solution of discrete system (24) calculated on a grid $\left\{x_i\right\}_{i = 0}^N$ generated by the above adaptive spatial grid algorithm. Then we have the following bound

Corollary 1. If we take $N^{-q}\leq C\Delta t$ with $0<q<1$, then we have

Now, we give the main result of our paper.

Theorem 6.4. Let $u(x,t)$ be the exact solution of the Burgers-Huxley equation (1), $\widetilde{u}^{n}(x)$ be the solution of the semi-discrete problem (10)-(12) after the time discretization and quasilinearization process and $U_i^{n}$ be the solution of fully discrete scheme (15) calculated on a grid $\left\{x_i\right\}_{i = 0}^N$ generated by the above adaptive spatial grid algorithm and time level $t_n = n \Delta t$. Assume that $N^{-q}\leq C\Delta t$ with $0<q<1$. Then we have the following bound

Proof. Let $E_i^n = u(x_i,t_n)-U_i^n$, $0\leq i\leq N$ be the global error at time level $t_n$, $1\leq n\leq M$. Then, splitting the global error $\{E_i^n\}$ as follows:

From (4) and Corollary 1, we have

For $\forall \varepsilon>0$, it follows from Theorem 3.3 that

Then, by applying the stability of the fully discrete scheme (15), we obtain

Finally, let $\varepsilon = C \Delta t(\Delta t+N^{-1+q})>0$, from (31)-(34), we obtain the following recurrence relation:

and hence, the result (30) follows from it.

7.   Numerical examples and discussion

In this section, we will give some numerical results obtained by the fully discrete scheme (15) for two test problems on the rectangular mesh $\overline{\Omega}_x^N\times S_t^M$, where $\overline{\Omega}_x^N$ is the equidistribution mesh generated from the numerical algorithm. It is well known that the constant $C_0$ in the above adaptive spatial mesh algorithm gives an intimation of how close we are from the equidistribution of the monitor function. In all the numerical experiments, we take $C_0 = 2$ and begin with $N = 32$ and the time step $\Delta t = 0.05$ and we multiply $N$ by two and divide $\Delta t$ by two.

For small enough values of the parameters $\varepsilon$, the exact solution of the following examples are not available. Therefore, we use the double mesh principal to calculate the maximum point-wise error, which are defined as

where $U^{N,M}(x_i,t_n)$ denote the numerical solution obtained on the mesh with $N$ mesh intervals in the spatial direction and $M$ mesh intervals in the time direction. From these values, we compute the corresponding order of convergence by

In addition, for the above Newton linearization process, we use the following convergence criterion for the numerical solution

Here, we chose $0$ as the initial guess in all cases and the iterations were stopped when the absolute error tolerance is achieved.

For comparison purposes, we use the upwind differences scheme (15) on the piecewise-uniform Shishkin mesh, which is constructed as follows[19,13]: Let $N$ be divisible by 2 and $\tau$ be the transition parameter which determines the point of transition from a fine mesh to the coarse mesh and is defined as

where $\theta$ is the constant whose value depends upon the method applied. Here, we choose $\theta = 1$. Then, we divide $[0,1]$ into two subdomains $[0,1-\tau]$ and $[1-\tau,1]$. Finally, we place $N/2$ number of subintervals in each of the subdomains and obtain the Shishkin mesh as follows:

Example 5.1. In this example, we consider the following singularly perturbed Burgers-Huxley equation

In Table 1, we give the maximum point-wise error $E_\varepsilon^{N,M}$ for $\varepsilon = 2^{-j}$, $j = 0,2,\cdots,18$. Meanwhile, the orders of convergence associated with the fully discrete scheme (15) are listed in Table 2. It is shown from these results that, with the increase of $N$ and $M$, the convergence order of the presented numerical method is more and more close to 1. To compare our results, the computed errors and the corresponding orders of convergence for the well-known layer-adapted meshes, Shishkin mesh, are given in Tables 3-4. Obviously, as can be seen from Tables 1-4 that adaptive grid method presented in this paper produced better results than that produced by using the Shishkin mesh.

In addition, the numerical solution is plotted in Figure 1 for $N = 64$, $M = 40$ and $\varepsilon = 2^{-14}$. Meanwhile, Figure 2 shows the computed solution at different time levels. It is shown from Figures 1-2 that the numerical solution of Example 5.1 has a boundary layer near $x = 1$ with the time variable $t\rightarrow 1$. Furthermore, in order to make the reader's understanding of the meshes generation by the above spatial grid algorithm, the process of grid movement after each iteration is plotted in Figure 3 for $N = 64$, $M = 40$ and $\varepsilon = 2^{-12}$, which should be read from bottom to top. The left of this figure is labeled with the value of $C_0$ for which the stopping criterion (19) becomes an equation. It is shown from Figure 3 that taking any smaller value of $C_0$ will increase the number of iterations. Obviously, the mesh starts to move toward the boundary $x = 1$ after each iteration. Thus, the presented adaptive grid method has the advantage that without any prior information of the boundary layer.

Example 5.2. Consider the following singularly perturbed Burger's equation:

For $\varepsilon = 2^{-j}$, $j = 0,2,\cdots,18$, the calculated maximum point-wise errors and the corresponding order convergence for Example 5.2 are listed in Tables 5-6, respectively. The numerical results given in Table 6 reveal the first-order uniform convergence with the increase of $N$ and $M$. Furthermore, we have also compared the numerical results using the presented method to that using the Shishkin mesh, see Tables 7-8. It is shown from Tables 5-8 that the presented adaptive grid method is better than the Shishkin mesh method. Figure 4 displays the numerical solution of Example 5.2 at different time levels. Figure 5 shows the mesh movement after each iteration. These figures also show the existence of the boundary layer at $x = 1$ with time variable $t\rightarrow 1$.

8.   Conclution

The present paper discusses a parameter-uniform numerical method for the singularly perturbed Burgers-Huxley equation. At first, for discretizing the time derivative, we use the classical backward-Euler method and for the spatial discretization the simple upwind scheme is used on a nonuniform grid. An iterative method based on Newton-Raphson-Kantorovich technique is presented to process the nonlinear term of the Burgers-Huxley equation. Parameter-uniform error estimations are derived for the numerical solution. It is also shown from the numerical results that the proposed numerical method has a better numerical accuracy compared to the Shishkin grid method. In all, for the numerical methods of the singularly perturbed Burgers-Huxley equation, the adaptive grid method is very effective.