The convection-diffusion-reaction equation is consists of three processes [19]. The first process is convection and is due to the movement of materials from one region to another. The second process is diffusion and is due to the movement of materials from a region of high concentration to a region of low concentration. The third process is reaction and is due to the decay, absorption and reaction of substances with other components. The convection-diffusion-reaction PDE provides a very useful and important mathematical model in a wide range of applications in sciences and engineering. Applications include the water quality problem in river networks [11], simulation of oil extraction from under-ground reservoirs [9], convective heat transport problems with large Peclet numbers [8], electromagnetic field problems in moving media [13], financial modeling of option pricing [2], turbulence models [16], drift diffusion equations of semiconductor device modeling [23], atmospheric pollution [24], fluid flow with high Reynolds numbers [20] and ground water transport [1]. In many of these applications, the unknown variables in the governing PDEs represent physical quantities that cannot take negative values such as pollutants, population, and concentration of chemical compounds [3].

Differential equations whose highest order derivative(s) is multiplied by a small perturbation parameter ε, 0 < ε ≪ 1 are called singularly perturbed differential equations [21]. Solutions of singularly perturbed problems, unlike regular problems, have a boundary and/or interior layers. The boundary layer is a narrow sub-domain specified by the parameter on which the solution varies by a finite value. The derivatives of the solution in this sub-domain grow without bound as ε tends to zero.

In the case of singularly perturbed problems, the use of numerical methods developed for solving regular problems leads to errors in the solution that depend on the value of the parameter ε. Errors of the numerical solution depend on the distribution of mesh points and become small only when the effective mesh-size in the layer is much less than the value of the parameter ε [18]. Such numerical methods turn out to be inapplicable for singularly perturbed problems. Because of this, there is an interest in the development of numerical methods where solution errors are independent of the parameter or that converge ε-uniformly. When the solutions of a PDE are ε-uniformly convergent, the methods and solutions are called robust [10]. Some ε-uniform numerical schemes developed for the considered problem can be found in [4, 5, 6, 12, 14, 28].

It is well known that classical numerical methods for solving singular perturbation problems are unstable and fail to give accurate results when the perturbation parameter ε is small. The accuracy and convergence of the methods need attention, because the treatment of singular perturbation problems is not trivial, and the solution depends on perturbation parameter and mesh size h [7]. This suggests that numerical treatment of singularly perturbed 1-D parabolic convection-diffusion problems should be improved. The presence of the singular perturbation parameter ε, leads to oscillations or divergence in the computed solutions while using classical numerical methods. To avoid these oscillations or divergence, an unacceptably large number of mesh points are required when ε is very small, which is not practical. So, to overcome this drawback associated with classical numerical methods, we develop a method based on the method of line (MOL) using a non-standard finite difference method in a spatial direction together with an implicit Runge-Kutta method of order two and three in the temporal direction, which treat the problem without creating an oscillation. Thus, this paper presents an accurate and ε-uniformly convergent numerical method for solving singular perturbation 1-D parabolic convection-diffusion problem.

The paper is organized as follows. In Section 1 a brief introduction about the problem is given, in Section 2 the definition of the problem and the behavior of its analytical solution is given. In Section 3, the discretization of the spatial domain and techniques of non-standard finite differences are discussed, and the ε-uniform convergence of the semi-discrete problem is proved. Next, the Runge-Kutta method is used for the system of IVPs resulting from the spatial discretization and the convergence of the discrete scheme is discussed. In Section 4, numerical examples and results are given to validate the theoretical analysis and finally in Section 5, the conclusion of the work done is given.

Notation: Throughout this paper M and K denote the number of mesh points in the space and time directions respectively. The symbol C is a positive constant independent of the perturbation parameter and mesh parameters M and K. The norms ‖.‖and‖.‖ΩM×QK are used to denote maximum norms defined as ‖g‖=maxx,t|g(x,t)|,(x,t)∈D and ‖gi,j‖ΩM×QK=maxi,j|g(xi,tj)|,0≤i≤M,0≤j≤K.

A singularly perturbed 1-D parabolic convection-diffusion problem on the domain D=Ωx×Q=(0,1)×(0,T] is given by:

where ε is the perturbation parameter such that 0 < ε ≪ 1, the coefficient functions a(x), b(x) and the source function f(x,t) are assumed sufficiently smooth. This condition guarantees the existence of a unique solution for the problem in Eq.(2.1).

In this paper, we assume the case a(x)≥ α > 0 and b(x)≥ β > 0 which enables the existence of the boundary layer on the right side of the domain.

2.1. Properties of Analytical Solution

Next, we see some of the properties of analytical solutions to the problem.

Let u0(x)∈C2[0,1] and μ0,μ1∈C1[0,T]. Impose the compatibility conditions u0(0)=μ0(0),  u0(1)=μ0(1) and

so that the data matches at the two corners points (0,0) and (1,0).

Let a(x), b(x) and f(x,t) be continuous functions on the domain D=Ω × Q. Then Eq.((2.1)) has unique solution u(x,t)∈ C² (D) [25].

Since we considered a right boundary layer problem using compatibility conditions, we deduce that there exists a constant C independent of ε such that ∀(x,t)∈D¯=[0,1]×[0,T], we have the following conditions that guarantee the existence of a constant C independent of ε such that ∀(x,t)∈D¯ given as:

For details, the interested reader can refer to page 105 of [25].

To show the bounds of the solution u(x,t) of Eq.(2.1), without loss of generality we assume the initial condition to be zero. Since u₀(x) is sufficiently smooth, using the property of the norm we can prove the following lemma.

Lemma 2.1.The bound on the solution u(x,t) of the continuous problem Eq.(2.1) is given by

Proof. From inequality |u(x,t)−u(x,0)|−|u(x,t)−u0(x)|≤Ct , we have

Since t ∈ [0,T] and u₀(x) is bounded it implies |u(x,t)| ≤ C.

Lemma 2.2.(Continuous maximum principle)Letψ∈C2,1(D¯)be such thatψ≥0,∀(x,t)∈∂D. Then for the differential operatorL=∂∂t−ε∂2∂x2+a(x)∂∂x+b(x),Lψ(x,t)>0,∀(x,t)∈Dimplies thatψ(x,t)≥0,∀(x,t)∈D¯.

Proof. Let (x*,t*) be such that ψ(x*,t*)=min(x,t)∈D¯ψ(x,t) and suppose that ψ(x*,t*)<0. It is clear that (x*,t*)∉∂D. So we have

Since

which implies ψx(x*,t*)=0,, ψt(x*,t*)=0 and ψxx(x*,t*)≥0 we get Lψ(x*,t*)<0 which contradicts the assumption made above. So we have Lψ(x*,t*)>0,∀(x,t)∈D. Hence ψ(x,t)≥0,∀(x,t)∈D¯..

Lemma 2.3.(Stability estimate)Let u(x, t) be the solution of problem in Eq.(2.1). Then we have the bound

where ‖f‖=maxx,t∈D|f(x,t)|..

Proof. Define barrier functions ϑ±(x,t) as

At the initial value:

At the boundary points:

For the differential operator:

which implies that Lϑ±(x,t)≥0,∀(x,t)∈D. Hence by maximum principle we have,

Hence, the proof is completed.

Lemma 2.4.The bound on the derivative of the solution u(x,t) of Eq.(2.1) with respect to x is given by

Proof. Interested reader can see the proof on [5].

3.1. Discretization in Spatial Direction

On the spatial domain [0, 1], we introduce uniform mesh with mesh length Δx=h such that ΩxM={xi}i=0M,x0=0,xM=1,h=1/M where M is the number of mesh points in the spatial discretization. For a problem in the form of Eq.(2.1), we consider the sub-equation obtained by neglecting the variable t. Following the non-standard finite difference formulation boundary value problems in [17] and [22] we get

First we rewrite Eq.(3.1) equivalently as a system of coupled first order differential equations:

Solving Eq.(3.3) we obtain y=exp(a(x)εx) in the discrete form: yi=exp(a(xi)εxi). To get the discrete difference scheme for y, we approximate Eq.(3.2) as

where U_i is denoted for the approximation of u(x) at grid point x_i while using the spatial discretization points.

Using the upwind finite difference for the first derivative, we obtain the scheme as

Now combining the Eqs.(3.3), (3.4) and (3.5) we solve for ρi2. We obtain the denominator function as:

Using ρi2 and Eq.(3.5) into the discretization of the main equation in Eq.(2.1), we obtain

In this discretization Eq.(2.1) reduces to semi-discrete form as:

where U0(t)=μ0(t),UM(t)=μ1(t), t∈[0,T] and L^h is a difference operator. The system of IVP in Eq.(3.6) can be written in compact form as:

where A is a tridiagonal coefficient matrix of size M−1×M−1 and Ui(t) and Fi(t) are M-1 size column vectors. The entries of A and F are given as:

and

respectively.

Now we need to show the semi-discrete operator L^h also satisfies the maximum principle and the uniform stability estimate.

Lemma 3.1.(Semi-discrete maximum principle)The operator defined by the discrete scheme in Eq.(3.6) satisfies a semi-discrete maximum principle. That is, SupposeU0(t)≥0,UM(t)≥0. Then LhUi(t)≥0,∀i=1(1)M−1 implies that Ui(t)≥0,∀i=0(1)M.

Proof. Suppose there exists s∈0,1,2,...M such that Us(t)=min0≤i≤MUi(t). Suppose that Us(t)<0 which implies s≠0,M. We have Us+1−Us>0 and Us−Us−1<0. We also have

Using the assumption, we get LhUi(t)<0 for i=1(1)M−1. Thus the supposition Ui(t)<0,i=1(1)M−1 is wrong. Hence Ui(t)≥0,∀i=0(1)M.

Lemma 3.2. The solution U_i(t) of the semi-discrete problem in Eq.(3.6) satisfies the following bound.

Proof. Let q=β−1max|LhUi(t)|+max{u0(xi),μ0(t),μ1(t)} and define the barrier function Ψi±(t) by: Ψi±(t)=q±Ui(t).

At the boundary points we have

On the discretized domain 0 < i < M, we have

From Lemma 3.1, using the semi-discrete maximum principle, we obtain

3.2. Convergence Analysis for Semi-discrete Scheme

In the above two lemmas we proved that the semi-discrete operator L^h satisfyies the maximum principle and the uniform stability estimate. In the next theorems we prove the ε-uniform convergence of the spatial discretization.

Theorem 3.1.Suppose the coefficients functions a(x), ;b(x) and the source function f(x,t) in Eq.(2.1) are sufficiently smooth, so that u(x,t)∈ C⁴([0,1] △ [0,T]). Then the difference of the solution u(x,t) and the semi-discrete solution U_i(t) of the Eq.(2.1) satisfies the following bound.

Proof. By considering the truncation error in spatial discretization we get:

Hence, we obtain the bound as:

Above we used the estimate ε|h2ρi2−1|≤Ch which based on the non-standard denominator function behavior used in [27]. Let γ=aih/ε,γ∈(0,∞). Then

Above we used the estimate ε|h2ρi2−1|≤Ch which based on the non-standard denominator function behavior used in [27]. Let γ=aih/ε,γ∈(0,∞), Then

where R(γ)=exp(γ)−1−γγ(exp(γ)−1), and from this we obtain the bounds limγ→0R(γ)=12 and limγ→∞R(γ)=0. Therefore R(γ)≤C,  γ∈(0,∞)

Using the boundedness of the derivatives of the solution in Lemma 2.4, with Eq.(3.11), we obtain:

Lemma 3.3.For a fixed mesh and for ε ⇒ 0, it holds

where x_j = jh,h = 1/M,∀j = 1(1)M - 1.

Proof. Consider the partition [0,1] : 0=x₀ < x₁ < ... < x_M= 1 for the interior grid points. We have

Since x₁ = h,1 - x _M-1 = h, repeated applications of L'Hospital's rule gives

This complete the proof.

Theorem 3.2.Under the hypothesis of boundedness of the semi-discrete solution from Lemma 3.3 and Theorem 3.1 above, the semi-discrete solution satisfies the following bound.

where ‖u(xi,t)−Ui(t)‖ΩM=max0≤i≤M|u(xi,t)−Ui(t)|.

Proof. This is immediate from the boundedness of the solution from Lemma 3.3 and Theorem 3.1 which the required estimates.

3.3. Discretization in Temporal Direction

On the time domain [0, T], we introduce the discretization with step size Δ t=t_j+1-t_j,j = 1(1)K so that Q^K be discretized domain where K is the number of mesh in the temporal direction. We use a lower order numerical scheme for discretizing the system of initial value problems in Eq.(3.7), and use a Runge-Kutta method developed by Bogacki and Shampine in [26]. First rewrite Eq.(3.7) in the form:

with the initial condition U(xi,0)=ϕ(xi),i=0(1)M. Here f(t,Ui(t))=−AUi(t)+Fi(t) so for each j = 1(1)K we write the scheme as:

where U_i,j denotes the approximation of U_i(t) at the grid point t_j. It is stated in [15] that, for j = 1(1)K, the local approximation U_i,j+1 to U_i(t_j+1) has third order accuracy (i.e. O(Δ t)³).

Lemma 3.4.From the above approximation method in the temporal direction, the global error estimates in this direction are given by

where E_j+1 is the global error in the temporal direction at time step (j+1)^th.

Proof. Using the local error estimate e_j up to the j^th time step, we obtain the global error estimate at the (j+1)^th time step.

Then using the boundedness of the solution, Lemma 3.4 implies

This shows that the discretization in temporal direction is consistent and global error is bounded. Now we use Eq.(3.14) to prove the ε-uniform convergence of the fully discrete scheme as

Using boundedness of the solution, Theorem 3.2, and Eq.(3.14), we obtain:

To validate the established theoretical results in this paper, we perform experiments using the proposed numerical scheme on the problem of the form given in Eq.(2.1).

Example 4.1. In this example, we consider the initial boundary value problem: ∂u∂t−ε∂2u∂x2+(2−x2)∂u∂x+xu(x,t)=10t2e−tx(1−x), , (x,t)∈(0,1)×(0,1] with initial condition u(x,0)=0, x∈[0,1], and boundary conditions u(x,0)=0, x∈[0,1], u(1,t)=0, t∈[0,1].

Example 4.2. In this example, we consider the initial boundary value problem: ∂u∂t−ε∂2u∂x2+(1+x(1−x))∂u∂x=f(x,t),  (x,t)∈(0,1)×(0,1] with initial condition u(x,0)=0,x∈[0,1], and boundary conditions u(0,t)=0, t∈[0,1], u(1,t)=0, t∈[0,1] where we choose the initial and the source functions f(x,t) are from the exact solution u(x,t)=e−t(c1+c2x−e−(1−x)ε) where c1=e−1ε and c2=1−e−1ε.

The exact solution is not known for the first example, therefore maximum nodal errors are calculated using the double mesh principle given in [27] as

where M is the number of mesh points in x and Δ t is the mesh length in the t direction. Let Ui,jM,Δt be the computed solution of the problem using mesh numbers M and Δ t, and let Ui,j2M,Δt/2 be the computed solution with twice as many (2M, 2K) mesh points which we get by adding the mid points xi+1/2 and t_j+1/2 into the mesh. For any values M and K the ε -uniform error estimate are calculated using the formula EM,Δt=maxε|EεM,Δt|.

The rate of convergence of the method is calculated using the formula

The solution of the problems given in Example 4.1 and 4.2 has a boundary layer at the right side of the x-domain (see Figures 1 and 2). The computed solutions U_i,j for different values of perturbation parameters are also shown in these figures. The numerical results displayed in tables 1 and 3 clearly indicate that the proposed method is ε-uniform convergent. From the results in these tables, we observe that the maximum point-wise error decreases as M increases for each value of ε. In addition, the maximum point-wise error is stable as ε ⇒ 0 for each M and Δt. Using computed results in these two examples, we confirm that the proposed numerical method is more accurate, stable and ε-uniform convergent with the rate of convergence one. The results in the proposed method are better than those given in [12] and [28].

In this paper, an ε-uniform numerical method has been developed for solving singularly perturbed 1-D parabolic convection-diffusion problems with a boundary layer on the right side of the domain. The developed method is based on the method of a line that constitutes the non-standard finite-difference for the spatial discretization and an implicit Runge-Kutta method of order 2 and 3 is used in the temporal direction for the system of initial value problem resulting from the spatial discretization. The stability and convergence of the proposed scheme are analyzed. Two model examples have been considered to validate the theoretical analysis by taking different values for the perturbation parameter ε. The computational results are presented in tables and figures. The proposed numerical scheme is first-order convergent. The performance of the scheme is investigated by comparing the results with prior studies. The proposed method gives more accurate and ε-uniformly convergent numerical results.