This work is licensed under the Creative Commons Attribution 4.0 International License.
Introduction
Let us take into consideration of the following singularly perturbed semilinear reaction-diffusion boundary value problem:
Lu\left( x \right) \equiv - \varepsilon {u^{\prime\prime}}\left( x \right) + f\left( {x,u\left( x \right)} \right) = 0,\;\;\;\;0 < x < l,u\left( 0 \right) = A,\;\;\;\;u\left( l \right) = B.
where ɛ, 0 < ɛ ≤ 1 is the perturbation parameter, f is given sufficiently smooth functions that f (x, u(x)) ∈ C ([0, l]xℝ),
{{\partial f\left( {x,u} \right)} \over {\partial u}}\; \ge \alpha > 0
. The problem (1)–(2) has boundary layers at the boundary points.
In a differential equation, if a small parameter is multiplied by the highest- order derivative term in the differential equation, generally it is called the singularly perturbed problem (denoted here by ɛ).[1,2,3]
Second-order reaction-diffusion type boundary value problems with singularly perturbed occur frequently in fluid mechanics and other fields of applied mathematics. Examples of such studied problems can be seen in [4,5,6].
Since the continuous solutions of singularly perturbed problems change very quickly in certain layers, it’s numerical analysis is always important. It is well known that when the parameter of perturbation is small enough, conventional numerical methods to solve the problem do not work well. Therefore should be develop develop appropriate numerical methods for such problems, whose convergence does not depend on the perturbation parameter. It can be found in the literature that there are many numerical finite difference schemes that are stable for all values parameter of perturbation [2,3,4,5,6,7,8,9]. One of the most important ways to easily find the methods that give such results is use of finite difference schemes with exponentially fitted [11,12,13].
As a numerical study, several examples of the second order singularly perturbed convection-diffusion problems can also be seen in [5],[21].
In [13–14,16,19] introduced numeric methods and special mesh methods for various reaction-diffusion type problem. In [20] semi-linear reaction diffusion equations are discussed. The discrete and upper of solutions were investigated for the asymptotic properties. And also it is numerical solutions are investigated on pisewise mesh.
The Numerov method is undoubtedly one of the most well known methods for reaction-diffusion type equations since it has fourth-order approach and it has been widely used in practical computational methods. Recently much fitted dinite difference scheme has been studied based on Numerov’s method in [11,12,13,14,15,16,17,18]. In [11], Phaneendra et al gave a finite difference Numerov scheme with a fitted multiplier three bands for solving singularly perturbed boundary value problem. Based on Numerov method, singularly perturbed nonlinear reaction-diffusion problem were investigated in [16,17,18,19].
In this study, we present finite difference scheme based on Numerov method for (1)–(2) problem on an uniform mesh. The some properties of the exact solution is given in section 2. In section 3, a finite difference scheme on a uniform mesh is introduced which is based on Numerov’s method. In section 4, the convergence of the approximate solution was presented and it was shown that uniform convergence was achieved at the discrete maximum norm. A numerical example and its results are given in section 5.
Notation
The C symbol in the throughout the article indicates a positive constants and does not have to be the same in each occurence which is independent of ɛ and of the mesh.
Properties of the exact solution
The semilinear equation (1) can be written in the take form;
Lu\left( x \right) \equiv - \varepsilon {u^{\prime\prime}}\left( x \right) + a\left( x \right)u\left( x \right) = F\left( x \right),0 < x < l,u\left( 0 \right) = A,\;\;u\left( l \right) = B,\;
where
a\left( x \right) = {{\partial f\left( {x,\widetilde u} \right)} \over {\partial u}}
, ũ = γu, 0 < γ < 1, F (x) = − f (x, 0).
Here we will give some important properties of the solution of (3)–(4) problem, which are required in later sections for the analysis of the numerical solution. We will indicate the maximum norm of any continuous g(x) function on the interval with ‖g‖∞.
Definition 1
(Maximum Principle). Assume that v(0) ≥ 0 and v(l) ≥ 0. Then Lv(x) 0, 0 < x < l, implies that v(x) ≥ 0, for all 0 ≤ x ≤ l.
The following two lemma and its solutions are given in [22].
Lemma 1
For any v(x) function, let v(x) ∈ C [0, l] ∩ C2 (0, l). Then the following estimate is true.
\left| {v(x)} \right| \le \left| {v(0)} \right| + \left| {v(l)} \right| + {\alpha ^{ - 1}}\,\mathop {\max }\limits_{1 \le i \le N} \left| {Lv\left( x \right)} \right|,\;\;\;0 \le x \le l.\;
Proof
Let us define the Ψ(x) function as follows:
\Psi \left( x \right) = \pm v\left( x \right) + \left| {v\left( 0 \right)} \right| + \left| {v\left( l \right)} \right| + {\alpha ^{ - 1}}\,\mathop {\max }\limits_{1 \le i \le N} \left| {Lv\left( x \right)} \right|,\;\;\;\;\;\;\;\;\;\;0 \le x \le l.\;
Then the following inequalities are satisfied\Psi \left( 0 \right) \ge 0,\;\;\Psi \left( l \right) \ge 0\;\;\;and\;L\Psi \left( x \right) \ge 0.
The maximum principle gives Ψ(x) ≥ 0, for all 0 ≤ x ≤ l, and so the inequalitiy (5) holds.
Lemma 2
Let a(x), F (x) are given sufficiently smooth functions and u(x) be the solution of the problem (3)–(4). Then the following estimates hold.
\parallel |u\left( x \right)|{|_\infty } \le C,\;\;\;\;\;\;\;\;\;\;0 \le x \le l.\;\left| {{u^\prime}(x)} \right| \le C\left\{ {1 + {1 \over {\sqrt \varepsilon }}\left( {{e^{ - {{\sqrt \alpha x} \over {\sqrt \varepsilon }}}} + {e^{ - {{\sqrt \alpha \left( {l - x} \right)} \over {\sqrt \varepsilon }}}}} \right)} \right\}.
Proof
Appliying Lemma 1. to (3)–(4) we have (8).
Lv\left( x \right) = \varphi \left( x \right),v\left( 0 \right) = O\left( {{1 \over {\sqrt \varepsilon }}} \right),\;\;\;\;\;\;\;\;\;\;v\left( l \right) = O\left( {{1 \over {\sqrt \varepsilon }}} \right),wherev\left( x \right) = {u^\prime}\left( x \right),\;\;\varphi \left( x \right) = {F^\prime}\left( x \right) - {a^\prime}\left( x \right)u\left( x \right).
The solution of the problem (10)–(11) has the following form:v\left( x \right) = {v_0}\left( x \right) + {v_1}\left( x \right).
Respectively, in here, the functions v0 (x) and v1 (x) are the solutions of the following problems:\left\{ {\matrix{ {L{v_0}\left( x \right) = \varphi \left( x \right),\,0 < x < l,} \hfill \cr {{v_0}\left( 0 \right) = {v_0}\left( l \right) = 0,} \hfill \cr } } \right.\left\{ {\matrix{ {L{v_1}\left( x \right) = 0,\,0 < x < l,} \hfill \cr {{v_1}\left( 0 \right) = {v_1}\left( l \right) = 0,} \hfill \cr } } \right.from Lemma 1., for the solution of the problem (14), we have\left| {{v_0}(x)} \right| \le {\alpha ^{ - 1}}\,\mathop {\max }\limits_{1 \le s \le l} \left| {\varphi (s)} \right|.
Thus, we obtain\left| {{v_0}(x)} \right| \le C,\;\;0 \le x \le l.
Applying maximum principle to the problem (15), we get\left| {{v_1}(x)} \right| \le w\left( x \right),where w(x) is the solution of the following problem:\left\{ {\matrix{ { - \varepsilon {w^{\prime\prime}} + \alpha {w^\prime} = 0,{\kern 1pt} \;{\kern 1pt} 0 < x < l,} \cr {w\left( 0 \right) = \left| {{v_1}(0)} \right|,\;\;w\left( l \right) = \left| {{v_1}(l)} \right|.} \cr } } \right.
The solution of this problem is has the formw\left( x \right) = {1 \over {\sinh \left( {{{\sqrt \alpha l} \over {\sqrt \varepsilon }}} \right)}}\left\{ {\left| {{v_1}(0)} \right|\sinh \left( {{{\sqrt {\alpha (} l - x)} \over {\sqrt \varepsilon }}} \right) + \left| {{v_1}(l)} \right|\sinh \left( {{{\sqrt \alpha x} \over {\sqrt \varepsilon }}} \right)} \right\},and it is from thatw\left( x \right) \le {C \over {\sqrt \varepsilon }}\left( {{e^{ - {{\sqrt \alpha x} \over {\sqrt \varepsilon }}}} + {e^{ - {{\sqrt \alpha \left( {l - x} \right)} \over {\sqrt \varepsilon }}}}} \right),is hold. Then combining (16),(17) and (20) in the following inequality, it can be easely obtained:\left| {{u^\prime}(x)} \right| \le \left| {{v_0}\left( x \right)} \right| + \left| {{v_1}(x)} \right|.
Thus the proof is completed.
Discretization and Mesh
In this section, we construct a numerical scheme for solving the problem (1)–(2) on a uniform mesh. Let wh denote the uniform mesh on [0, l].
{w_h} = \left\{ {{x_i} = ih,\;\;\;i = 1,2, \ldots ,N - 1,\;\;h = 1/N} \right\},\;\;{\bar w_h} = {w_h} \cup \left\{ {x = 0,l} \right\}.
Let us show wi = w(xi) for any function w(x), and moreover any approximation of the function u(x) at point xi with yi. We will use the following notations for any mesh function {wi} defined on
{\bar w_N}
:
{w_{\bar x,i}} = {{{w_i} - {w_{i - 1}}} \over h},\;\;{w_{x,i}} = {{{w_{i + 1}} - {w_i}} \over h},\;\;\;\;\;{w_{\bar xx,i}} = {{{w_{x,i}} - {w_{\bar x,i}}} \over h} = {{{w_{i + 1}} - 2{w_i} + {w_{i - 1}}} \over {{h^2}}},\;\;\;
and
{\left\| w \right\|_{C\left( {{{\bar w}_h}} \right)}}: = \mathop {\max }\limits_{1 \le i \le N} \left| {{w_i}} \right|.
To find the difference approach that corresponds to (1), let us use the following identity and use the interpolating quadrature formulas in [8] on each intervals (xi−1, xi) and (xi, xi+1),
\chi _i^{ - 1}{h^{ - 1}}\int\limits_{{x_{i - 1}}}^{{x_{i + 1}}} Lu(x){\varphi _i}(x)dx = 0,\,\,1 \le i \le N - 1,
then we obtain the following relation:
l{w_i} \equiv - \varepsilon {\theta _i}{w_{\bar xx,i}} + f({x_i},{w_i}) = {R_i},\;\;i = 1,2, \ldots ,N - 1,
where
{R_i} = \chi _i^{ - 1}{h^{ - 1}}\left\{ {\int_{{x_{i - 1}}}^{{x_{i + 1}}} (f\left( {x,w} \right) - f({x_i},{w_i}){\varphi _i}\left( x \right)dx} \right\}.
In here θi is called fitting factor and after a simple calculation, the value of
{\theta _i} = \chi _i^{ - 1} = 1
.
if Ri omited in equation (22) then we have numerical scheme for (1)–(2)\left\{ {\matrix{ {l{w_i} \equiv - \varepsilon {w_{\bar xx,i}} + f({x_i},{w_i}) = 0,{\kern 1pt} \;{\kern 1pt} i = 1,2, \ldots ,N - 1,} \hfill \cr {{w_0} = A,\;\;{w_N} = B.} \hfill \cr } } \right.
We note that
\varphi _i^{\left( 1 \right)}
and
\varphi _i^{\left( 2 \right)}
are basis functions such that they are solutions of the following problems respectively,
-{\varepsilon \varphi _i^{(1)}}^{\prime\prime} = 0,\;\varphi _i^{(1)}\left( {{x_i}} \right) = 1,\;\varphi _i^{(1)}\left( {{x_{i - 1}}} \right) = 0, - {\varepsilon \varphi _i^{(2)}}^{\prime\prime} = 0,\;\varphi _i^{(2)}\left( {{x_i}} \right) = 1,\;\varphi _i^{(2)}\left( {{x_{i + 1}}} \right) = 0,
and which are
{\varphi _i}\left( x \right) = \left\{ {\matrix{ {\varphi _i^{\left( 1 \right)} \equiv {{\left( {x - {x_{i - 1}}} \right)} \over h};\;\;\;{x_{i - 1}} < x < {x_i},} \cr {\varphi _i^{\left( 2 \right)} \equiv {{\left( {{x_{i + 1}} - x} \right)} \over h};\;\;\;{x_i} < x < {x_{i + 1}},} \cr {0\;{\kern 1pt} \;\;\;\;\;\;\;\;\;\;\;\;\;,\;\;x \notin \left( {{x_{i - 1}},{x_{i + 1}}} \right),}} } \right.
and also
{\chi _i} = {h^{ - 1}}\int_{{x_{i - 1}}}^{{x_{i + 1}}} {\varphi _i}\left( x \right)dx = 1
.
Description of the Numerov’s method
For convenience, let yi = y(xi),
{y^{\left( n \right)}}\left( {{x_i}} \right) = y_i^{(n)}
at x = xi, for any function y(x), using the Taylor series expansion, we obtain:
{{{y_{i + 1}} - 2{y_i} + {y_{i - 1}}} \over {{h^2}}} = y_i^{\prime\prime} + {1 \over {12}}{h^2}y_i^{\prime\prime\prime\prime} + O({h^4}).
If the error term is omited, and instead of
y_i^{\prime\prime}
and
y_i^{(4)}
the corresponding terms are written, then we have the following a Numerov finite difference scheme:
\left( {{\varepsilon \over {{h^2}}} - {{{a_{i - 1}}} \over {12}}} \right){y_{i - 1}} - \left( {{{2\varepsilon } \over {{h^2}}} + {{10{a_i}} \over {12}}} \right){y_i} + \left( {{\varepsilon \over {{h^2}}} - {{{a_{i + 1}}} \over {12}}} \right){y_{i + 1}} = - {1 \over {12}}\left( {{F_{i - 1}} + 10{F_i} + {F_{i + 1}}} \right).\;i = 1,2, \ldots ,N - 1.
From [10],
{A_i} > 0,\;\;{B_i} > 0,\;\;{C_i} - {A_i} - {B_i} > 0,
are valid and the siystem has only one solution. This system can be solved by Thomas algorithm.
An Algorithm for Numerov Type Scheme
If Numerov method is applied to problem (1)–(2) we obtain following equation:
\left\{ {\matrix{ { - \varepsilon {w_{\bar xx,i}} + \left[ {1 + {1 \over {12}}{h^2}{w_{\bar xx,i}}} \right]f\left( {{x_j},{w_j}} \right) = 0,\;\;\;\;j = 1,2, \ldots ,N;} \cr {{w_0} = A,\;\;{w_{N + 1}} = B} \cr } } \right.
With the help of the Newton-Raphson-Kantorovich approach, we’ll get a new Numerov type difference scheme. Instead of f (xj, uj) function, in the equation for (28), let’s write the following equivalent,
f\left( {{x_j},w_j^{(k)}} \right) = f\left( {{x_j},w_j^{(k - 1)}} \right) + {{\partial f} \over {\partial w}}\left( {{x_j},w_j^{\left( {k - 1} \right)}} \right).(w_j^{(k)} - w_j^{(k - 1)})
then we obtain the following relation.
- \varepsilon {w_{\bar xx}}w_j^{\left( k \right)} + \left[ {1 + {1 \over {12}}{h^2}{w_{\bar xx}}} \right]\left\{ {f\left( {{x_j},w_j^{\left( {k - 1} \right)}} \right) + {{\partial f} \over {\partial w}}\left( {{x_j},w_j^{\left( {k - 1} \right)}} \right).\left( {w_j^{\left( k \right)} - w_j^{\left( {k - 1} \right)}} \right)} \right\} = 0,\;j = 1,2, \ldots ,N;\;k = 1,2, \ldots
Let zi = yi − ui, 0 ≤ i ≤ N, where yi the solution of (22) and ui the solution of (1)–(2) at mesh point xi. We now estimate the approximate error zi, which satisfies the following discrete problem
l{z_i} = {R_{i\;}};\;{z_0} = {z_N} = 0
where the truncation error Ri is in the equation (23).
Definition 2
(Discrete Maximum Principle). Suppose that a mesh function vi satisfies v0 ≥ 0 and vN ≥ 0. Then ℓvi ≥ 0, for all 0 ≤ i ≤ N − 1 implies that vi ≥ 0 for all 0 ≤ i ≤ N.
Theorem 3
Let f (x, u) ∈ C2 [0, l]. Then the following estimate holds.
{\left\| {y - u} \right\|_{C\left( {{{\bar w}_h}} \right)}} \le C{h^2}
Proof
If we find an estimate for Ri the proof is complete. For function f (x, u), using the Taylor series expansion, we obtain:
\matrix{ {f\left( {x,u} \right) - f\left( {{x_i},{u_i}} \right) = \left( {x - {x_i}} \right)\left\{ {{{\partial f\left( {{x_i},{u_i}} \right)} \over {\partial x}} + {{\partial f\left( {{x_i},{u_i}} \right)} \over {\partial u}}{{du\left( {{x_i}} \right)} \over {dx}}} \right\} + } \cr { + {{{{(x - {x_i})}^2}} \over {2!}}\left\{ {{{{\partial ^2}f({\xi _i},u({\xi _i}))} \over {\partial {x^2}}} + {{{\partial ^2}f({\xi _i},u({\xi _i}))} \over {\partial x\partial u}}{{du\left( {{\xi _i}} \right)} \over {dx}} + {{{\partial ^2}f({\xi _i},u({\xi _i}))} \over {\partial {u^2}}}{{({{du\left( {{\xi _i}} \right)} \over {dx}})}^2}} \right\}}}
If we write this expression instead of Ri we obtain the following relation:
\matrix{ {\left( {x - {x_i}} \right)\left\{ {{{\partial f\left( {{x_i},{u_i}} \right)} \over {\partial x}} + {{\partial f\left( {{x_i},{u_i}} \right)} \over {\partial u}}{{du\left( {{x_i}} \right)} \over {dx}}} \right\} + } \cr { + {{{{(x - {x_i})}^2}} \over {2!}}({{{\partial ^2}f({\xi _i},u({\xi _i}))} \over {\partial {x^2}}} + \;{{{\partial ^2}f({\xi _i},u({\xi _i}))} \over {\partial x\partial u}}{{du\left( {{\xi _i}} \right)} \over {dx}} + {{{\partial ^2}f({\xi _i},u({\xi _i}))} \over {\partial {u^2}}}{{({{du\left( {{\xi _i}} \right)} \over {dx}})}^2}){\varphi _i}\left( x \right)dx}}
Considering the equivalents of u′ and u″, using of discrete maximum principle we have
{\left\| R \right\|_{C\left( {{{\bar w}_h}} \right)}} \le C{h^2}
This estimate conclude the proof of Theorem 4.
Numerical example
To illustrate the applicability of the method proposed in this article, we applied it to an example. Connsider the following semilinear problem:
Example 4. − ɛu″ − exp(−(x2 + u)) = 0, x ∈ [0, 1], u(0) = 0, u(1) = 1
For numerical approximation of solution, we have shown that the method uniform convergens according to the ɛ-parameter. Since the exact solution to this problem could not be found, we used the following the double mesh prensiple for calculate of the maximum absolute errors.
E_\varepsilon ^N = \mathop {\max }\limits_{0 \le i \le N} \left| {u_i^N - u_{2i}^{2N}} \right|
For any N, the ɛ-uniform maximum absolute error is calculated by
{E^N} = \mathop {\max }\limits_\varepsilon {E_\varepsilon }.
The numerical rate of convergence and ɛ-uniform convergence rate for example has been calculated by the following formulas;
R_\varepsilon ^N = {{\log \left| {E_\varepsilon ^N/E_\varepsilon ^{2N}} \right|} \over {log 2}},\,\,{R^N} = {{\log \left| {{E^N}/{E^{2N}}} \right|} \over {log 2}}
The maximum point wise errors and the rates of convergence of the problem in example is presented in Table 1. (
u_i^0 = {x^2}
, 1 ≤ i ≤ N − 1, for arbitrary inital function)
Maximum errors and rates of convergence wh
ɛ
N = 32
N = 64
N = 128
N = 256
N = 512
R_\varepsilon ^N
2−5
2.3e–4
5.9e–5
1.5e–5
3.7e–6
9.2e–7
2.009
2−6
3.8e–4
9.5e–5
2.4e–5
6.2e–6
1.5e–6
2.014
2−7
6.4e–4
1.6e–4
4.0e–5
1.0e–5
2.5e–6
2.019
2−8
1.1e–3
2.4e–4
6.9e–5
1.7e–5
4.3e–6
2.020
2−9
1.9e–3
4.9e–4
1.2e–4
3.1e–5
7.6e–6
2.015
RN
2.013
2.003
2.002
1.998
2.020
Conclusion
In this paper we have presented a Numerov’s scheme to solve a class of singularly perturbed semilinear reaction-diffusion problem. We have introduced computational technique based Numerov’s scheme on a uniform mesh. Uniform convergence of the method is demonstrated with respect to the parameter of perturbation. The accuracy of the uniform convergence was supported by a numerical example.
