Article
Kyungpook Mathematical Journal 2021; 61(3): 591-612
Published online September 30, 2021
Copyright © Kyungpook Mathematical Journal.
Higher Order Uniformly Convergent Numerical Scheme for Singularly Perturbed Reaction-Diffusion Problems
Worku Tilahun Anilay, Gemechis File Duressa, Mesfin Mekuria Woldaregay*
Department of Mathematics, Jimma University, Jimma, Ethiopia
e-mail : workutil12@gmail.com
Department of Mathematics, Jimma University, Jimma, Ethiopia
e-mail : gammeef@gmail.com
Department of Applied Mathematics, Adama Science and Technology University, Adama, Ethiopia
e-mail : msfnmkr02@gmail.com
Received: April 10, 2020; Revised: January 14, 2021; Accepted: February 8, 2021
Abstract
In this paper, a uniformly convergent numerical scheme is designed for solving singularly perturbed reaction-diffusion problems. The problem is converted to an equivalent weak form and then a Galerkin finite element method is used on a piecewise uniform Shishkin mesh with linear basis functions. The convergence of the developed scheme is proved and it is shown to be almost fourth order uniformly convergent in the maximum norm. To exhibit the applicability of the scheme, model examples are considered and solved for different values of a singular perturbation parameter ε and mesh elements. The proposed scheme approximates the exact solution very well.
Keywords: Finite Element, Fitted Mesh, Parameter Uniform, Singularly Perturbed.
1. Introduction
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
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
A numerical solution
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
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,
positive constant independent of ε and
2. Description of the Problem
Consider a class of singularly perturbed reaction-diffusion problems of the form:
subjected to the boundary conditions
where
2.1. Properties of the exact solution
Let us denote the differential equation in (2.1)-(2.2) by the differential operator
Lemma 2.1. Assumethat
Consequently,
The next lemma gives a bound on the derivatives of the solution of the considered problem.
Lemma 2.2. Let
Sharper bounds on the solution and its derivatives are required in the proof of the error estimates. To do that, the solution
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. Numerical Scheme Formulation
3.1. Finite element method on Shishkin mesh
Let
where
and
so that
Since the considered problem exhibits two boundary layers on the left and right side, the domain
where τ is the Shishkin mesh transition parameter given by
The mesh size
Using the mesh points in (3.2), we construct a set of piecewise linear basis function
for
Let us consider a typical element
where
where the coefficients
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
where
where
called the stiffness matrix and
called the load vector.
Since the base functions are linear, each element
Here, the subscripts 1 and 2 are labels of the endpoint nodes on a typical element
The linear system of equations for the entire mesh becomes
where
and
Applying the boundary conditions
and the two auxiliary equation corresponding to nodes 1 and
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
So, the reduced stiffness matrixin (3.18) satisfies the criteria of
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
where
so that on each
Using Taylor series expansion the interpolation error on
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
In the second case, we have
Now using (3.23) it follows that
Case 2: For the outer layer region. For
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
Let denote
So, we have
since
In similar manner on double number of mesh elements
Combining (3.31) and (3.32) for eliminating the term
which is the Richardson extrapolated approximate solution. The error for the approximate solution in (3.33) becomes
4. Numerical Results and Discussion
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)
with the boundary conditions
with the boundary conditions
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
-
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 andN=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.
5. Conclusion
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.
Acknowledgements.
The author would like to thank the referees for their constructive comments and suggestions.
References
- E. B. Becker, G. F. Carey and J. T. Oden, Finite Elements, An Introduction: Volume I, Printice-Hall Inc, Englewood Cliffs, New jersey, 1981.
- M. Chandru, T. Prabha & V. Shanthi, A Hybrid Difference Scheme for a Second-Order Singularly Perturbed Reaction-Diffusion Problem with Non-smooth Data, Int.J. Appl. Comput. Math, 1(2015), 87-100. https://doi.org/10.1007/s40819-014-0004-8
- C. Clavero, R. K. Bawa and S. Natesan, A robust second-order numerical method for global solution and global normalized flux of singularly perturbed self-adjoint boundary-value problems, Int. J. Comput. Math., 86(10-11)(2009), 1731-1745.
- G. File and F. Edosa, Higher order compact finite difference method for singularly perturbed one dimensional reaction-diffusion problems, J. Nigerian Math. Soc., 36(3)(2017), 491-502.
- G. File and Y. N. Reddy, Domain Decomposition Method for Solving Singular Perturbation Problems, International Journal of Applied Science and Engineering, 11(4)(2013), 433-448.
- F. W. Gelu, G. F. Duressa and T. A.Bullo, Sixth-order compact finite difference method for singularly perturbed 1D reaction diffusion problems, J. of Taibah University for Science, 11(2)(2017), 302-308.
- F. W. Galu, G. F. Duressa and T. A. Bullo, Tenth Order Compact Finite Difference Method for Solving Singularly Perturbed 1D Reaction-Diffusion Equations, International Journal of Applied Science and Engineering, 8(3)(2016), 15-24.
- M. K. Kadalbajoo and P. Arora, B-Spline collocation method for a singularly perturbed reaction-diffusion problem using artificial viscosity, Int. J. Comput. Methods, 6(1)(2009), 2341.
- N. Kopteva and E. O’Riordan, Shishkin meshes in the numerical solution of singularly perturbed differential equations, Int. J. Numer. Anal. Model., 1(1)(2004) 1-18.
- M. Kumar and S. C. Rao, High order parameter-robust numerical method for singularly perturbed reaction-diffusion problems, Appl. Math. Comput., 216(4)(2010) 1036-1046.
- J. J. Miller, H. MacMullen, E. O’Riordan and G. I. Shishkin, A second-order parameter-uniform overlapping Schwarz method for reaction-diffusion problems with boundary layers, J. Comput. Appl. Math., 130(1-2)(2001), 231-244.
- J. J. Miller, E. O’Riordan and G. I. Shishkin, Fitted numerical methods for singular perturbation problems: error estimates in the maximum norm for linear problems in one and two dimensions, World Scientific, New Jersey, 2012.
- S. Natesan, R. K. Bawa and C. A. Clavero, Uniformly convergent compact numerical scheme for the normalized flux of singularly perturbed reaction-diffusion problems, Int. J. Inf. Syst. Sci., 3(2)(2007), 207-221.
- S. Natesan and R. K. Bawa, Second-order numerical scheme for singularly perturbed reaction-diffusion Robin problems, SIAM J Numer Anal, 2(3-4)(2007), 177-192.
- A. H. Nayfeh, Introduction to perturbation techniques, John Wiley & Sons, New York, 2011.
- K. Phaneendra, K. M. Latha and Y. N. Reddy, A Fitted arithmetic average three-point finite difference method for singularly Perturbed two-point boundary value Problems with dual layers, International Journal of Applied Science and Engineering, 13(3)(2015), 205-216.
- S. C. Rao and M. A. Kumar, A uniformly convergent exponential spline difference scheme for singularly perturbed reaction-diffusion problems, Neural Parallel Sci. Comput., 18(2)(2010) 121-136.
- J. N. Reddy, An introduction to the finite element method, McGraw-Hill, New York, 1993.
- E. R. O’Malley, Singular perturbation methods for ordinary differential equations, Springer-Verlag, New York, 1991.
- S. Russell, Sparse grid methods for singularly perturbed problems, PhD thesis, National University of Ireland, 2016.
- M. K. Singh, S. Natesan, Richardson extrapolation technique for singularly perturbed system of parabolic partial differential equations with exponential boundary layers, Appl. Math. Comput., 333(2018) 254-75.
- F. Wondimu, G. File and T. Aga, Fourth order compact finite difference method for solving singularly perturbed 1D reaction diffusion equations with Dirichlet boundary conditions, Momona Ethiopian J. of Science, 8(2)(2016) 168-181.
- L. Zhilin, Q. Zhonghua and T. Tao, Numerical Solutions of Differential Equations: Introduction to Finite Difference and Finite Element Methods, Cambridge University Press, New york, 2018.