Abstract

This paper considers the numerical treatment of singularly perturbed time-dependent convection-diffusion-reaction equation. The diffusion term of the equation is multiplied by a small perturbation parameter (ε), which takes an arbitrary value in the interval (0, 1]. For small values of ε, the solution of the equation exhibits an exponential boundary layer which makes it difficult to solve it analytically or using classical numerical methods. We proposed numerical schemes using the Crank–Nicolson method in time derivative discretization and the nonstandard finite difference method (exact finite difference method) in space derivative discretization on a uniform and piecewise uniform Shishkin mesh. The existence of unique discrete solutions and the stability of the schemes are discussed and proved. Uniform convergence of the schemes is proved. The formulated schemes converge uniformly with linear order of convergence. The method on Shishkin mesh possesses boundary layer resolving property. We validated the methods by considering two numerical examples for different values of and mesh length.

1. Introduction

Differential equations are commonly used for modeling processes in various applied mathematics areas. It relates one or more functions and their derivatives; the functions represent physical quantities, the derivatives represent their rates of change, and the differential equations define the relationship between the two. Many authors are working on the analytical and numerical solutions of different types of differential equations that arise from different application areas[1, 2].

The convection-diffusion-reaction equations provide a very useful and important mathematical model in a wide range of disciplines. The application areas include fluid flow with high Reynolds numbers, convective heat transport problems with large Peclet numbers, turbulence models, simulation of oil extraction from underground reservoirs, financial modeling of option pricing, the water quality problem in river networks, drift diffusion equations of semiconductor device modeling, atmospheric pollution; and so on [3]. In many of these applications, the unknown variables in the governing equation represent physical quantities that cannot take negative values, such as concentrations of chemical compounds, pollutants, and population [4]. In the listed scientific areas, the models are most often described by singularly perturbed differential equations.

Convection dominated reaction-diffusion equations or singularly perturbed convection-diffusion-equations are differential equations in which the diffusion term (or the highest order derivative term) of the equation is multiplied by a small perturbation parameter . In general, the solutions of such equations exhibit boundary layer behaviour. Thus, the solutions of singularly perturbed equations possess a multiscale property. The presence of both slow and fast phenomena makes the problem stiff from the standpoint of an approximate solution.

To obtain an approximate solution for singularly perturbed problems, the methods are mainly classified into two classes: (i) the asymptotic (or analytical) approach and (ii) the numerical approach. The numerical approach provides quantitative information about the solution of a particular problem, while the asymptotic approach provides the qualitative behavior of the solution of a family of equations and semiquantitative information about the solution of any particular member of the family [5].

There are two strategies for designing numerical schemes that give small errors inside the boundary layer region. The first approach is the class of fitted mesh methods, which consists of choosing a fine mesh in the boundary layer region and a coarse mesh in the outer layer region. The second approach is the fitted operator methods, in which the mesh remains uniform and the difference schemes reflect the qualitative behaviour of the solution inside the boundary layer region. A nice discussion of using one or both of the above strategies can be found in [3, 69].

The fitted operator methods can be categorized into two classes based on the techniques of numerical scheme formulation. That is, the exponentially fitted operator FDM and the nonstandard FDM. In nonstandard FDM the discretized mesh is uniform and a complicated denominator function is developed for the diffusion part of the equation using the exact finite difference technique to handle the effect of the perturbation parameter. This method can be extended to solve higher dimensional and nonlinear singularly perturbed problems. Different authors have published papers on the designing and implementing of nonstandard FDM on uniform mesh for different classes of model equations [10]. In the exponentially fitted operator FDM, the discretized mesh is uniform and the exponential fitting parameter is induced using the asymptotic solution behaviour to control the rapid growth or decay of the solution in the boundary layer region. The main drawback of the exponentially fitted operator FDM is that it does not work for nonlinear equations and higher dimensional equations.

In the case of the fitted mesh methods, one uses nonuniform grids, which are fine in the boundary layer region and coarse outside the layer. In designing appropriate numerical methods using the fitted mesh method, one must have some a priori information about sensitive regions where the solution undergoes some rapid changes. Different mesh generation techniques have been developed so far. The well-known meshes are the exponentially graded Bakhvalov mesh and the piecewise uniform Shishkin mesh [6].

Standard numerical methods on uniform mesh fail to give accurate results and create oscillations in the computed solution due to the presence of a small perturbation parameter on the diffusion term [11, 12]. To avoid the oscillations and obtain an accurate solution using the standard numerical methods on uniform mesh, it requires an unacceptably large number of mesh points when is very small. So, this becomes unpractical due to round-off error and the limited size of the computer’s memory. To obtain an accurate and uniformly convergent numerical solution on a small number of mesh points, authors use the fitted mesh methods or the fitted operator methods.

Nowadays, the numerical treatment of singularly perturbed convection-diffusion-equations has attracted great attention [4]. The authors in [1216] used the fitted mesh methods for solving different classes of singularly perturbed problems, and the authors in [1719] used the fitted operator methods for developing uniformly convergent numerical schemes.

The nonstandard FDM on Shishkin mesh has better accuracy and order of convergence than the equivalent standard FDM on Shishkin mesh [20]. In addition, the uniform convergence analysis of the nonstandard FDMs were restricted to uniform mesh [7]. In [21], He and Wang developed a new form of the nonstandard FDM for solving a convection dominated diffusion equation. In the study of [20], Woldaregay developed numerical schemes using the nonstandard FDM on uniform mesh and Shishkin mesh for solving singularly perturbed boundary value problem. Motivated by the results in [20, 21], in this paper we proposed a uniformly convergent numerical scheme for solving the problem in (1). The proposed scheme uses the Crank–Nicolson method for time derivative discretization and the midpoint upwind nonstandard FDM on the uniform mesh and the Shishkin mesh for the spatial derivative discretization. Furthermore, we discuss the convergence analysis of the proposed scheme. The proposed nonstandard FDM on the Shishkin mesh gives a boundary layer resolving result.

1.1. Organization of the Paper

The remaining parts of the paper are organised as follows: in Section 2, the problem under consideration is defined; in Section 3, the numerical schemes and their convergence analysis are given. In Section 4, the numerical result and discussion are included; in Section 5, the conclusion is given.

1.2. Notations

and denotes the number of mesh intervals in space and time discretization, respectively; C is denoted for positive constant independent of and . The norm used in this article is the maximum or supremum norm.

2. Continuous Problem

On the domain with the boundary where , we considered a singularly perturbed time dependent convection-diffusion-reaction equation.where and is the perturbation parameter.

In (1), the term is called the diffusion term, is called the convection term and is called the reaction term [9]. The functions , and are assumed to be sufficiently smooth and bounded. We assume that the coefficient functions and satisfieswhere and are constants, for the existence of a boundary layer on the right side of the spatial domain.

The existence and uniqueness of the solution of (1) can be established by assuming that the data is Holder continuous and by imposing appropriate compatibility conditions at the corner points. Using the assumptions of sufficiently smoothness of and . The required compatibility condition at the corner points are and

So that the data matches at the corner points and . For and be continuous on domain , then (1) has unique solution .

Lemma 1. The solution of (1) satisfies the following estimate:where is constant independent of .

Proof. The result follows directly from the compatibility condition. See the detailed proof in [22].
The reduced problem obtained by setting in (1) is given by the following equation:For small values of , the solution of (1) is very close to the solution of (5).

Lemma 2 (The maximum principle). Let be the solution of (1) which satisfies and implies that .

Proof. The proof is by contradiction. Suppose that there exist such that . It is clear that i.e. . Since , which implies , and . Giving that , which is a contradiction to the assumption made . Therefore, .

Lemma 3 (Stability estimate). The solution of (1) satisfies the boundedwhere is lower bound of .

Proof. By defining the barrier functions as , on the boundary lines we obtain , and
On the domain , we have the following equation:With the hypothesis of the maximum principle, we obtain . Hence, the proof is completed.
Using Shishkin decomposition in [12, 23], the solution of (1) can be decomposed into regular and singular component as follows:where denotes the regular component and denotes the singular component. The regular component satisfies the following nonhomogeneous equation:and the singular component satisfies the following homogeneous equation:

Lemma 4 (see [14]). The derivatives of the regular and singular component solutions respectively satisfies the boundand the derivatives with respect to and of the solution of (1) satisfies the boundedwhere is lower bound of .

3. The Numerical Schemes

3.1. Time Semidiscretization

We divide the time domain into subintervals using the grid points , where . We denote for the approximation of in temporal discretization at grid point. Using the Crank–Nicolson finite difference method for semidiscretizing (1), we obtain a system of boundary value problems

for , where, . The semidiscrete scheme in (13) satisfies the maximum principle which is stated as follows.

Lemma 5 (Semidiscrete maximum principle). Let be smooth function on . If and , then .

Proof. Similar techniques as the proof in Lemma 2 gives the required result.
Next, we find a bound for the truncation error in temporal semidiscretization. Let us denote the local truncation error for each time step by .

Lemma 6 (Local error estimate). The local error in the temporal semidiscretization is bounded as follows:

Proof. Using the Taylor series expansion we obtain the local truncation error of the semidiscrete scheme as follows:Since and and . Using the semidiscrete maximum principle we obtain the required bound.
The global error estimate up to time step is bounded as follows:

3.2. Spatial Discretization

In this subsection, we approximate the spatial derivatives using nonstandard FDM on uniform mesh and on the Shishkin mesh. Moreover, we discuss the convergence analysis of the schemes.

3.2.1. Exact Finite Difference

For constructing the exact finite difference scheme for convection diffusion equation, we follow the techniques developed in [4, 7]. Consider a constant coefficient homogeneous convection diffusion equationswhere and . The characteristic equation of (17) is given by the following equation:and its solution is given as . The constant coefficient equation in (17) has two independent solutions which is given as follows:

Consider a uniform mesh , where is the number of mesh intervals in spatial discretization. We want to develop a difference equation which has the same general solution as the differential equation in (18) has at the mesh point , which is given by the following equation:

From the theory of difference equations for second order linear difference equations, we obtain the following equation:rearranging gives the following equation:

Simplifying and substituting the values of , we obtain the following equation:which is an exact difference scheme for (18). Rearranging, we obtain the following equation:

The denominator of the second derivative approximation is and the denominator of the first derivative approximation is .

3.2.2. Spatial Discretization via Nonstandard FDM on Uniform Mesh

We modify for the variable coefficient differential equations of the form in (13) as the following form:where is denoted for . Using the denominator function and applying the midpoint upwind finite difference discretization the required difference schemes becomes the following equation:

for , whereand and denoted for , and , respectively.

3.2.3. Stability and Convergence Analysis on Uniform Mesh

Here, we want to show that the discrete scheme in (27) satisfies the maximum principle, stability estimates, and uniform convergence.

Lemma 7. Let be any mesh function satisfying . Then, implies that .

Proof. Suppose there exist such that . Suppose that which implies . So, we have that and . So, we obtain for . Thus, the supposition , for is wrong. Hence, we obtain .
Using the discrete maximum principle, we want to prove the discrete scheme satisfies the uniform stability result.

Lemma 8. The solution of the discrete scheme in (27) satisfy the bound.

Proof. Let and define barrier functions by . At the boundary points, we have the following equation:On the discretized domain , we obtain the following equation:By the discrete maximum principle in Lemma 7, we obtain . Hence the required bound is satisfied.

Lemma 9. The discrete solution of the scheme in (27) satisfy the bound

Proof. The proof is similar to the above Lemma 8 proof by construction a barrier functions .
Let us denote the forward and backward finite differences operators in spatial variable as follows:respectively and the second order finite difference operator as follows:

Lemma 10. For a fixed mesh and for , it holdswhere .

Lemma 11 (see [24]). For all , the following estimate holds

Theorem 1. The truncation error of nonstandard FDM on uniform mesh satisfies the bound

Proof. The spatial discretization truncation error is given by the following equation:Using the result in Lemma 11, the truncation error bound becomes the following equation:Using the bound in Lemma 4 into (39), we obtain the following equation:

Theorem 2. The nonstandard FDM on uniform mesh satisfies the error bound

Proof. Using the results in Lemma 10 and applying the bound in Lemma 9, we obtain the error bound.

Theorem 3. Let and be the solution of the problem in (1) and (27), respectively, then the nonstandard FDM on uniform mesh satisfies the error bound

Proof. Combining the error bound for the temporal and spatial discretization in (16) and Theorem 2 gives the required bound.

3.3. Spatial Discretization via Nonstandard FDM on Shishkin Mesh

The development of nonstandard FDM and the uniform convergence of the schemes were restricted on a uniform mesh. In this paper, we extend the scheme and uniform convergence of the scheme on piecewise-uniform Shishkin mesh. Let be the discretization of the spatial domain where be the number of grid points be an even positive integer. For each , we define . Since the considered problem exhibits the right boundary layer, we set a mesh transition parameter as follows:

Which divides the domain [0, 1] in to outer layer region and boundary layer region , where is a constant that represents the order of the scheme and is the lower bound of . Here, we set a uniform mesh in with mesh spacing and similarly a uniform mesh is placed in with the mesh spacing . Hence, the mesh point is given as follows:

The nonstandard FDM are commonly applied on uniform mesh, so using it on layer adapted mesh requires the following change on the denominator function:

For the mesh function at the grid points , we denote the second derivatives term approximation as follows:

For , the midpoint upwind nonstandard FDM on Shishkin mesh is given as follows:where

3.3.1. Stability and Convergence Analysis on the Shishkin Mesh

Let be the mesh width in the outer layer region and be the mesh width in layer region . Then, the mesh widths satisfies the following conditions:

Lemma 12 (The discrete maximum principle). Let be any mesh function satisfying . Then, implies that .

Proof. The proof is similar to the proof of Lemma 7.

Lemma 13. The solution of the discrete scheme in (47) satisfy the bound.

Proof. Let and define barrier functions by . At the boundary points, we have the following equation:On the discretized spatial domain , we have the following equation:Using the maximum principle in Lemma 12, we obtain . Hence, the required bound is satisfied.

Lemma 14. The discrete solution of the scheme in (47) satisfy the bound

Proof. The proof is similar to the above-given lemma proof by the construction of barrier functions .

3.3.2. Decomposition of the Discrete Solution

Here, we decompose the numerical solution into regular and singular components using the Shishkin decomposition in [12] as follows:where the regular part satisfied the following nonhomogeneous equation:and the singular component satisfies the following homogeneous equation:

The error in the computed solution can also be decomposed as follows:

First, let us consider the case which is a uniform mesh case say the mesh size is . The truncation error becomes the following equation:

Since

From the truncation error in uniform mesh, we obtain the following equation:

Using Lemma 14, we obtain

For the case , we estimate the error in the regular and singular solutions separately. Here, implies that