Many physical phenomena are modeled by parameter dependent differential equations. The behaviour of the solution depends on the magnitude of the parameter. A differential equation in which the highest order derivative term is multiplied by a small positive parameter ε∈(0,1] is called a singularly perturbed differential equation and the parameter ε is called singular perturbation parameter [15], [19]. Singularly perturbed problems (SPPs) arises in the modeling of fluid dynamics, elasticity, quantum mechanics, reaction-diffusion process, chemical-reactor theory, plasma dynamics, meteorology, diffraction theory, aerodynamics, oceanography, modeling of semi-conductor, hydrodynamics and many other similar areas [4].

The solutions of singularlyperturbed differential equation exhibit a multi-scale character; i.e. there are thin layers or regions of the domain where the solution changes rapidly or jumps suddenly, forming boundary layers, while away from the layers the solution behaves regularly or changes slowly in the outer region. As a result such problems are called boundary layer problems [5]. Due to the multi-scale character of the solution, classical numerical methods such as FDM, FEM and the collocation method are inefficient. This happens since the error estimates of numerical methods depend explicitly on the derivatives of the solution and the derivatives are unbounded as ε→0 [20]. Early numerical solutions of singularly perturbed differential equations were obtained using a standard finite difference method on a uniform mesh. In this approach, as the perturbation parameter decreases in magnitude the mesh needs to be refined sufficiently to capture the boundary layer(s). Hence, such methods are inefficient and inaccurate.

A numerical solution U is said to be ε-uniformly convergent to the exact solution u, if there exist a positive integer N₀ and positive numbers C and p, where N₀, C and p are all independent of N and ε, such that, for all N>N₀, sup0<ε≤1‖U−u‖≤CN−p. Here p is the order of the method and N is number of mesh elements [11], [12].

Different authors have developed several numerical methods for solving singularly perturbed reaction-diffusion problems. Some of these can be found in [4], [5], [6, 7], and [22]. The results in these papers are good in terms of the value of ε and mesh size h. The main drawback of the schemes in the above listed papers is that the methods fail to give good results as ε→0, which means that the methods are not ε-uniformly convergent. The recent papers [2], [3], [10], [13, 14], [16] and [17] contain ε-uniform numerical methods developed for solving singularly perturbed reaction-diffusion problems.

In this paper, we developed a fitted mesh finite element method with Richardson extrapolation for solving second order singularly perturbed two point boundary value problems of the reaction-diffusion type. The proposed method involves dividing the domain of the solution into a finite number of elements and using variational concepts to construct an approximate solution over the collection of elements. The main reason behind seeking an approximate solution on a collection of sub-domains is the fact that it is easier to represent a complicated function as a collection of simple polynomials [1], [18]. Most numerical methods developed so far for solving singularly perturbed reaction-diffusion problems are not parameter-uniformly convergent. So developing a higher order numerical method whose convergence does not depend on the perturbation parameter has a great importance to the scientific research area[9]. This paper deals with formulating a higher order uniformly convergent method to find numerical solution of singularly perturbed 1D reaction-diffusion problems.

Notations. In this paper, N denotes the number of mesh elements and C is a

positive constant independent of ε and N. The notation ‖⋅‖ denotes maximum norm.

Consider a class of singularly perturbed reaction-diffusion problems of the form:

subjected to the boundary conditions

where 0<ε≪1 is the singular perturbation parameter and α and γ are given constants. Assume that the functions f(x) and b(x) are sufficiently smooth and b(x)≥β>0 for some constant β. As the problem is of the reaction-diffusion type, the condition b(x)≥β>0 ensure the existence of dual boundary layers near the two end points of the solution domain.

2.1. Properties of the exact solution

Let us denote the differential equation in (2.1)-(2.2) by the differential operator L as L=−εd2dx2+b(x). The considered singularly perturbed problem in (2.1)-(2.2) satisfies the maximum principle.

Lemma 2.1. Assumethat ψ(0)≥0 and ψ(1)≥0. Then, Lψ(x)≥0,∀x∈Ω, implies that ψ(x)≥0, ∀x∈Ω¯.

Proof. Assume that there exists x*∈Ω¯=[0,1] such that ψ(x)∈C2(Ω)∪C0(∂Ω) and ψ(x*)=minΩ¯ψ(x)<0. Thus, from the conditions given on the boundaries, it is clear that x*∉{0,1} which implies that x*∈(0,1). It follows from the assumption that ddxψ(x*)=0 and d2dx2ψ(x*)≥0.

Consequently, Lψ(x*)=−εd2dx2ψ(x*)+b(x)ψ(x*)<0. This is a contradiction and hence we can conclude that ψ(x)≥0, ∀x∈Ω¯.

The next lemma gives a bound on the derivatives of the solution of the considered problem.

Lemma 2.2. Let u(x)∈C2(Ω)∪C0(∂Ω) be the solution of the problem in (2.1)-(2.2). Then, its derivatives satisfies the bound

Proof. The proof is carried out by constructing barrier functions and applying the maximum principle. The details of the proof is given in Appendix section.

Sharper bounds on the solution and its derivatives are required in the proof of the error estimates. To do that, the solution u of (2.1)-(2.2) is decomposed as u=v+wL+wR, where v denote the regular component of the solution and wL and w_R respectively denotes the left and right singular components of the solution. This decomposition is known as Shishkin decomposition [12].

Lemma 2.3. The derivatives of the regular component solution satisfies the bound

and the derivatives of the singular components solution satisfies the bound

The detail proof is given in Appendix section.

3.1. Finite element method on Shishkin mesh

Let H01 denote the set of all functions whose order 1 or less is square integrable over Ω¯=[0,1] and vanishes at the end point of the domain. Using the weak formulation, we rewrite (2.1)-(2.2) as

where

and

so that B(u,v) is a bilinear function in u and v and L(v) is a linear functional of v, ∀v∈H01. We call (3.1) is weak form (or variational form) of (2.1)-(2.2). Interested reader can refer the weak formulation of boundary value problems in [1].

Since the considered problem exhibits two boundary layers on the left and right side, the domain Ω¯=[0,1] is subdivided into boundary layer regions (0,τ) and (1-τ,1) and outer layer region (τ,1-τ). In this paper, we use the Shishkin piece-wise uniform mesh in order to have more mesh points in the boundary layer regions. Let N be the number of elements on the domain Ω¯, then we define mesh discretization as Ω¯τN={xi}i=1N+1 with hi+1=xi+1−xi for i=1,2,3,...,N+1. The sub-domains (0,τ) and (1-τ,1) are discretized into N4 mesh elements and the sub-domain (τ,1-τ) is discretized into N2 mesh elements as

where τ is the Shishkin mesh transition parameter given by

The mesh size h_i in each subdomain is given by

Using the mesh points in (3.2), we construct a set of piecewise linear basis function ϕi(x)∈H01 of the form;

for i=1,2,...,N, which are often called the hat functions.

Let us consider a typical element Ωτe=[x1e,x2e] from the domain ΩτN using the above discretization. Here e denote the element number and x1e and x2e denote the left and right end of the element Ωτe. Thus, the weak form of the problem over each element is given as

where ue(x) is denoted for the restriction of u(x) on the element Ωτe=[x1e,x2e]. Representing the numerical solution by the linear combination of the basis function on each element as

where the coefficients cje are unknowns to be determined, Ne is the number of nodes in Ωτe and ϕje(x) are the basis functions on the element Ωτe. Differentiating (3.7) once gives

Plugging the approximate solution in (3.7) and its derivative in (3.8) into (3.8) gives

Applying Galerkin method i.e. taking the test function v(x)=ϕie(x) for i=1,2 gives

where ϕie(x) are a piecewise linear base functions. We rewrite (3.10) as

where

called the stiffness matrix and

called the load vector.

Since the base functions are linear, each element Ωτe has two nodes (i.e. degree of freedom) and there are two equations per element of the form

Here, the subscripts 1 and 2 are labels of the endpoint nodes on a typical element Ωeτ. These subscripts are to be relabeled upon assembling the elements so as to coincide with appropriate nodes 1, 2, 3,..., N+1 in the final mesh. The equations on the entire elements assembled as follows. Since the mesh contain N elements and N+1 nodes, we have N+1 equations in N+1 degree of freedom describing the assembled system of elements.

The linear system of equations for the entire mesh becomes

where

and

Applying the boundary conditions uN(0)=α and uN(1)=γ then N-1 unknown nodal values c₂, c₃, c₄,..., c_N remain. The equation in (3.15) reduces to N-1 system of equations as

and the two auxiliary equation corresponding to nodes 1 and N+1 are

3.2. Stability of the scheme

From (3.16) we see that the entries of the stiffness matrix in (3.18) are given by

On each row of the matrix in (3.18), we obtain Ki,i−1<0, Ki,i+1<0 and Ki,i>0 since ϕi(x) and ϕi+1(x) are linear polynomials with oppositely signed slopes. In addition, for each row i=2,3,...,N it satisfies the diagonal dominance ([23], page 152)

So, the reduced stiffness matrixin (3.18) satisfies the criteria of M-matrix. Hence, it is nonsingular. This guarantees the stability of the scheme. So, the unknown nodal values c2, c3, c4,..., cN are solved easily using Thomas algorithm [1].

3.3. Convergence analysis of the scheme

In this section, we prove the uniform convergence of the developed scheme using the maximum norm. The next theorem gives the ε-uniform or parameter uniform error bound of the scheme in the maximum norm.

Theorem 3.1. Let U_i be the numerical solution and u(x_i) be exact solution of the problem in (2.1)-(2.2) on the Shishkin mesh Ω¯τN. Then

where C is a constant independent of ε and N.

Proof. The bound is obtained separately on each element Ωi=[xi−1,xi]. Note that for any function g on Ωi, we write the approximate solution using the basis functions ϕi(x) as

so that on each Ωi we have

Using Taylor series expansion the interpolation error on Ωi given by

giving the error in the weak form as (for detail see on [23], Theorem 7.11)

Using Lemma 1.2, we obtain

Since the coefficient matrix is invertible from stability condition, so we have

Similar to the decomposition of continuous solution we decompose the approximate solution as

Using the decomposition in (3.24) and Lemma 1.3, we obtain

Now, the argument depends on the parameter τ. That is, whether τ=14 or τ=2εβlnN. In case τ=14≤2εβlnN which implies 1ε≤Cβ(lnN)2. Here since the mesh is uniform, we have hi=N−1 implying that hi2=N−2. Thus using (3.23), it follows that

In the second case, we have τ=2εβlnN≤14 which implies 1ε≥Cβ(lnN)2. We consider separately the inner layer regions and the outer layer region. Case 1: For the inner layer (boundary layer) regions. That is for 1≤i≤N4 and 3N4+1≤i≤N then

Now using (3.23) it follows that

Case 2: For the outer layer region. For N4+1≤i≤3N4 then τ≤1−xi so that e−(1−xi)βε≤N−2. Similarly e−xiβε≤N−2 since τ≤xi which implies that

Then combining the outer layer region bound, left and right boundary layer bounds in (3.25), (3.26) and (3.27) gives

Since the right side of (3.28) is independent of ε, taking suprimum on both sides of (3.28) we obtain the required bound.

3.4. Richardson extrapolation

The aim of the Richardson extrapolation is to increase the order of convergence by combining the discrete solutions calculated on different meshes. For detail analysis of Richardson extrapolation on fitted mesh one can refer [21]. To apply the technique, we solve the discrete problem (3.18) on the fine mesh element Ω¯τ2N with 2N mesh-intervals on piecewise Shishkin mesh having the same transition points τ. The corresponding mesh widths in Ω¯τ2N becomes

Let denote UiN and Ui2N an approximate solution on N and 2N number of mesh elements respectively. From (3.23) we have

So, we have

since 1ε≤Cβ(lnN)2.

In similar manner on double number of mesh elements 2N

Combining (3.31) and (3.32) for eliminating the term N−2(lnN)2 gives

which is the Richardson extrapolated approximate solution. The error for the approximate solution in (3.33) becomes

To verify the established theoretical results in this paper, we perform an experiment using the proposed numerical scheme on the problem of the form given in (2.1)-(2.2)

Example 4.1. Consider the singularly perturbed problem

with the boundary conditions u(0)=0, u(1)=0. Its exact solution isu(x)=1−xe(x−1ε)−(1−x)e(−xε).

Example 4.2. Consider the singularly perturbed problem

with the boundary conditions u(0)=0, u(1)=0. Its exact solution is u(x)=1−e(−xε)−e(−(1−x)ε).

The maximum absolute error is calculated using

The ε-uniform error is calculated using

The rate of convergence is given as

The maximum absolute error of Example 4.1 and 4.2 before Richardson extrapolation is given in Tables 1 and 5 respectively for different values of singular perturbation parameter. In these tables the result are computed for the inner boundary layers(I) and outer boundary layer(O) region separately. Similarly, in Tables 2 and 6 the maximum absolute error of Example 4.1 and 4.2 after Richardson extrapolation is given respectively. On these tables one can observe that, as ε goes small the maximum absolute error for inner layer and outer layer region becomes stable and bounded. This indicates that maximum absolute error of the scheme is independent of the perturbation parameter ε, implying that the scheme is ε-uniformly convergent. One can observe from these tables, the method with Richardson extrapolation is more accurate than that of before the Richardson extrapolation. In Table 3, we compare the result of the proposed scheme with the results in [8]. As one observes in this table the proposed scheme gives more accurate result. In Tables 4 and 7, the rate of convergence of the proposed scheme before and after Richardson extrapolation is given. The proposed scheme before extrapolation is almost second order uniformly convergent and the scheme after the extrapolation gives almost fourth order uniform convergence. The computed result is in good agreement with the theoretical finding.

Table 1 . Maximum absolute error of Example 4.1 in inner layer(I) and outer layer(O) region before Richardson extrapolation..

ε↓		N=24	25	26	27	28	29
10-3	I	2.0159e-02	7.0588e-03	2.5475e-03	8.6985e-04	2.5648e-04	6.4050e-05
	O	3.0700e-03	9.2902e-04	2.1628e-04	2.7825e-05	1.8248e-06	4.5649e-07
10-4	I	1.8932e-02	6.6543e-03	2.3997e-03	8.1752e-04	2.6683e-04	8.4387e-05
	O	3.9895e-03	2.4233e-03	1.0726e-03	3.5717e-04	1.1995e-04	3.9066e-05
10-5	I	1.8545e-02	6.5265e-03	2.3530e-03	8.0097e-04	2.6155e-04	8.2713e-05
	O	4.1588e-03	2.5426e-03	1.2422e-03	5.2323e-04	1.8505e-04	6.4679e-05
10-6	I	1.8422e-02	6.4861e-03	2.3383e-03	7.9574e-04	2.5988e-04	8.2183e-05
	O	4.2122e-03	2.5745e-03	1.2593e-03	5.3739e-04	2.1131e-04	7.7506e-05
10-7	I	1.8383e-02	6.4733e-03	2.3336e-03	7.9408e-04	2.5935e-04	8.2016e-05
	O	4.2292e-03	2.5846e-03	1.2647e-03	5.3998e-04	2.1258e-04	7.9837e-05
EN		1.8365e-02	6.4674e-03	2.3314e-03	7.9332e-04	2.5911e-04	8.1938e-05

Table 2 . Maximum absolute error of Example 4.1 in inner layer(I) and outer layer(O) region after Richardson extrapolation..

ε↓		N=24	25	26	27	28	29
10-3	I	7.0588e-03	8.6985e-04	6.4050e-05	4.0017e-06	2.5010e-07	1.5632e-08
	O	9.2902e-04	2.7825e-05	4.5649e-07	2.8536e-08	1.7837e-09	1.6811e-10
10-4	I	6.6543e-03	8.1752e-04	8.4387e-05	7.8780e-06	6.8770e-07	5.7247e-08
	O	2.4233e-03	3.5717e-04	3.9066e-05	4.0666e-06	4.1394e-07	4.0840e-08
10-5	I	6.5265e-03	8.0097e-04	8.2713e-05	7.7212e-06	6.7399e-07	5.6058e-08
	O	2.5426e-03	5.2323e-04	6.4679e-05	6.7397e-06	4.1874e-07	4.0911e-08
10-6	I	6.4861e-03	7.9574e-04	8.2183e-05	7.6716e-06	6.6967e-07	5.5682e-08
	O	2.5745e-03	5.3739e-04	7.7506e-05	6.7231e-06	4.1681e-07	4.4345e-08
10-7	I	6.4733e-03	7.9408e-04	8.2016e-05	7.6559e-06	6.6828e-07	5.5642e-08
	O	2.5846e-03	5.3998e-04	7.9837e-05	6.7231e-06	4.1681e-07	4.4345e-08
EN		7.0588e-03	8.6985e-04	8.4387e-05	7.8536e-06	6.6828e-07	5.5642e-08

Table 3 . Example 4.1, Comparison of maximum absolute error..

ε↓	N=24	25	26	27	28	29
	Proposed	Scheme
2-12	7.0588e-03	8.6985e-04	6.4050e-05	4.0017e-06	2.5010e-07	1.5632e-08
2-16	6.6543e-03	8.1752e-04	8.4387e-05	7.8780e-06	6.8770e-07	5.7247e-08
2-20	7.0588e-03	8.6985e-04	8.4387e-05	7.8536e-06	6.6828e-07	5.5642e-08
2-25	7.0588e-03	8.6985e-04	8.4387e-05	7.8536e-06	6.6828e-07	5.5642e-08
	Result	in [8]
2-12	4.245e-03	1.957e-03	5.262e-04	1.227e-04	3.012e-05	7.496e-06
2-16	1.302e-03	1.292e-03	1.061e-03	4.893e-04	1.315e-04	3.068e-05
2-20	3.255e-04	3.255e-04	3.255e-04	3.231e-04	2.653e-04	1.223e-04
2-25	5.754e-05	5.754e-05	5.754e-05	5.754e-05	5.754e-05	5.752e-05

Table 4 . Rate of convergence of Example 4.1 before and after Richardson extrapolation..

ε↓	N=24	25	26	27	28
	Before RE
2-12	1.5100	1.4719	1.5518	1.6157	1.6610
2-16	1.5072	1.4749	1.5513	1.6149	1.6608
2-20	1.5059	1.4758	1.5512	1.6143	1.6611
2-25	1.5061	1.4760	1.5513	1.6145	1.6607
	After RE
2-12	3.0206	3.7635	4.0005	4.0000	3.9999
2-16	3.0250	3.2762	3.4211	3.5180	3.5865
2-20	3.0206	3.3657	3.4256	3.5548	3.5862
2-25	3.0206	3.3657	3.4256	3.5548	3.5862

Table 5 . Maximum absolute error of Example 4.2 in inner layer (I) and outer layer(O) region before Richardson extrapolation..

ε↓		N=24	25	26	27	28	29
10-3	I	1.8365e-02	6.4674e-03	2.3314e-03	7.9332e-04	2.3420e-04	5.8489e-05
	O	3.3885e-03	1.1233e-03	2.6567e-04	3.5484e-05	1.8514e-06	4.6313e-07
10-4	I	1.8365e-02	6.4674e-03	2.3314e-03	7.9332e-04	2.5911e-04	8.1938e-05
	O	4.2368e-03	2.5689e-03	1.1405e-03	3.8306e-04	1.2970e-04	4.2701e-05
10-5	I	1.8365e-02	6.4674e-03	2.3314e-03	7.9332e-04	2.5911e-04	8.1938e-05
	O	4.2370e-03	2.5893e-03	1.2671e-03	5.3499e-04	1.8961e-04	6.6528e-05
10-6	I	1.8365e-02	6.4674e-03	2.3314e-03	7.9332e-04	2.5911e-04	8.1938e-05
	O	4.2370e-03	2.5893e-03	1.2672e-03	5.4118e-04	2.1300e-04	7.8195e-05
10-7	I	1.8365e-02	6.4674e-03	2.3314e-03	7.9332e-04	2.5911e-04	8.1938e-05
	O	4.2370e-03	2.5893e-03	1.2672e-03	5.4118e-04	2.1311e-04	8.0063e-05
EN		1.8365e-02	6.4674e-03	2.3314e-03	7.9332e-04	2.5911e-04	8.1938e-05

Table 6 . Maximum absolute error of Example 4.2 in inner layer (I) and outer layer(O) region after Richardson extrapolation..

ε↓		N=24	25	26	27	28	29
10-3	I	6.4674e-03	7.9332e-04	5.8489e-05	3.6546e-06	2.2841e-07	1.4276e-08
	O	1.1233e-03	3.5484e-05	4.6313e-07	2.8951e-08	1.8096e-09	1.6806e-10
10-4	I	6.4674e-03	7.9332e-04	8.1938e-05	7.6487e-06	6.6766e-07	5.5580e-08
	O	2.5689e-03	3.8306e-04	4.2701e-05	4.5447e-06	4.7328e-07	4.7806e-08
10-5	I	6.4674e-03	7.9332e-04	8.1938e-05	7.6487e-06	6.6766e-07	5.5531e-08
	O	2.5893e-03	5.3499e-04	6.6528e-05	8.0126e-06	9.5752e-07	1.1448e-07
10-6	I	6.4674e-03	7.9332e-04	8.1938e-05	7.6487e-06	6.6766e-07	5.5515e-08
	O	2.5893e-03	5.4118e-04	7.8195e-05	9.6277e-06	1.1834e-06	1.4564e-07
10-7	I	6.4674e-03	7.9332e-04	8.1938e-05	7.6487e-06	6.6765e-07	5.5590e-08
	O	2.5893e-03	5.4118e-04	8.0063e-05	1.0292e-05	1.2827e-06	1.5936e-07
EN		2.5893e-03	5.4118e-04	8.1938e-05	1.0292e-05	1.2827e-06	1.5936e-07

Table 7 . Rate of convergence of Example 4.2 before and after Richardson extrapolation..

ε↓	N=24	25	26	27	28
	Before RE
10-3	1.5057	1.4721	1.5550	1.6144	1.7602
10-4	1.5057	1.4721	1.5550	1.6144	1.6609
10-5	1.5057	1.4721	1.5550	1.6144	1.6609
10-6	1.5057	1.4721	1.5550	1.6144	1.6609
10-7	1.5057	1.4721	1.5550	1.6144	1.6609
	After RE
10-3	3.0272	3.7617	4.0004	4.0000	4.0000
10-4	3.0272	3.2753	3.4212	3.5180	3.5865
10-5	3.0272	3.2753	3.4212	3.5180	3.5877
10-6	3.0272	3.2753	3.4212	3.5180	3.5877
10-7	3.0272	3.2753	3.4212	3.5180	3.5877

In Figure 1, the profile of the solutions is given for different values of the perturbation parameter ε=2^-4, 2^-8 and 2^-12 with boundary layer formulation with layer thickness of O(ε) as ε goes small. In Figures 2 and 3, we plot the computed solutions of Example 4.1 and 4.2 on uniform mesh and Shishkin mesh. In Figures 2(a) and 3(a), the computed solution oscillates in the boundary layer regions while using uniform mesh, whereas in case of Shishkin mesh the exact and computed solution agree well in the boundary layer regions. In addition to that, one observes in Figures 2(b) and 3(b), sufficient number of mesh points exist in the boundary layer region for ε=2^-20 and N=64. This indicate layer resolving property of the developed scheme. In Figure 4, the Log-Log scale plot of maximum absolute error versus the number of meshes points are given for different values of the perturbation parameter ε ranging from 2^-12 to 2^-25.

Figure 1. Solutions profile with boundary layer formation as ε goes small on (a) Example 4.1 on (b) Example 4.2.

Figure 2. Numerical and Exact solutions of Example 4.1 for ε=2^-20 and N=64, (a) on uniform mesh, (b) on Shishkin mesh.

Figure 3. Numerical and Exact solutions of Example 4.2 for N=64 and ε= 10^-4 , (a) on uniform mesh, (b) on Shishkin mesh.

Figure 4. The graph of maximum absolute error versus the number of elements, using Log-Log scale plot on (a) Example 4.1, on (b) Example 4.2.

Uniformly convergent finite element method with Richardson extrapolation technique is presented for solving singularly perturbed reaction-diffusion problems. The stability of the scheme is discussed. Theoretically the scheme is shown ε-uniformly convergent with order of convergence almost two before Richardson extrapolation. After the extrapolation technique is applied the order of convergence of the scheme is accelerated to almost four. Model examples are considered by taking different values for ε and the results are presented in tables and graphs. The obtained result shows that the proposed method is uniformly convergent, approximate the exact solution very well and is in a good agreement with the theoretical results of the analysis. The proposed scheme also works well for variable coefficient problems.

The author would like to thank the referees for their constructive comments and suggestions.