- Email: [email protected]
A modified graded mesh and higher order finite element approximation for singular perturbation problems
Journal of Computational Physics 395 (2019) 275–285
Contents lists available at ScienceDirect
Journal of Computational Physics www.elsevier.com/locate/jcp
A modified graded mesh and higher order finite element approximation for singular perturbation problems Aditya Kaushik a,c , Anil K. Vashishth b , Vijayant Kumar b , Manju Sharma d,∗ a
Department of Applied Mathematics, Delhi Technological University, Delhi – 110042, India Department of Mathematics, Kurukshetra University, Kurukshetra – 136119, India University Institute of Engineering and Technology, Panjab University, Chandigarh – 160014, India d Department of Mathematics, KVA DAV College, Karnal – 132001, India b c
a r t i c l e
i n f o
Article history: Received 20 February 2018 Received in revised form 31 March 2019 Accepted 11 April 2019 Available online 12 June 2019 Keywords: Singular perturbation Boundary value problem Layer-adapted meshes Finite element method
a b s t r a c t A singularly perturbed convection diffusion problem is solved numerically using finite element method based on higher order polynomials. More precisely, we introduce a modified graded mesh generated using some implicitly defined functions. Higher order parameter uniform convergence is obtained in -weighted energy norm. Moreover, the error estimates obtained are optimal in the sense that they are free from logarithmic factor. A number of test examples are taken into account and a rigorous comparative analysis is presented. Moreover, we compare the proposed method with others found in the literature. © 2019 Elsevier Inc. All rights reserved.
1. Introduction It is well known that many phenomena in biology [14,15], chemistry [21], engineering [30,10], physics [23] can be described by boundary value problems associated with various types of ordinary and partial differential equations or systems. When we associate a mathematical model with a phenomenon, we generally try to capture what is essential, retaining the important quantities and omitting the negligible ones which involve small parameters. The model that would be obtained by maintaining the small parameters is called the perturbed model, whereas the degenerate simplified model is called unperturbed or reduced model [2]. Of course, the unperturbed model is to be preferred, because it is simpler. What matters is that it should describe faithfully enough the respective phenomenon, which means that its solution must be sufficiently close enough to the solution of the corresponding perturbed model [2]. This fact holds in the case of regular perturbations. On the other hand, in the case of singular perturbations, things get more complicated. If we refer to an initial boundary value problem, the solution of the unperturbed problem does not satisfy in general all the original boundary conditions and/or initial conditions as some of the derivatives may disappear by neglecting the small parameters [2,28]. Thus, some discrepancy may appear between the solution of the perturbed model and that of the corresponding reduced model [31]. Numerical discretization of singularly perturbed differential equation is burdened with difficulties; for example, finite difference method (FDM) and finite element method (FEM) necessitate a mesh to sustain approximations [19,28]. In fact, the discrete solutions obtained using standard Galerkin or cantered finite differences methods demonstrate oscillatory behaviour for small discretization parameters. However, layer adapted meshes seem promising in the discretization of such equations [19]. This leads to growing interest in the generation of layer-adapted meshes. Layer adapted meshes have first been pro-
*
Corresponding author. E-mail addresses: [email protected] (A. Kaushik), [email protected] (A.K. Vashishth), [email protected] (M. Sharma).
https://doi.org/10.1016/j.jcp.2019.04.073 0021-9991/© 2019 Elsevier Inc. All rights reserved.
276
A. Kaushik et al. / Journal of Computational Physics 395 (2019) 275–285
posed by Bakhvalov [1,24,5] in the context of reaction-diffusion problems. Later on, special meshes for convection-diffusion problems has been extensively investigated. That includes Shishkin’s piecewise equidistant meshes [19,17], grids that are based on equidistribution [16,3,22], Gartland type meshes [11], Bakhvalov-Shishkin grids [18,37] and Vulanovic´ improved Shishkin grids [34,35]. Further development leads to the layer-adapted meshes defined using a recursive formula. Examples are the Gartland-Shishkin meshes [27] or the graded meshes analysed in [9,8,6,4,36] and others. The first parameter uniform convergence result, for linear finite elements on a Shishkin mesh, in the -weighted H 1 -norm was proved in [33]. The given interval is divided into two subintervals that includes a small subinterval [0, τ ] and a large subinterval [τ , 1]. Here, the choice of τ = τ0 ln N is important. However, the mesh is equidistant in both subintervals but, fine in [0, τ ] and coarse in [τ , 1] with a mesh size of order O ( N −1 ). Moreover, the resulting error estimate u − u N ≤ C N −1 ln N is not optimal due to the presence of logarithmic factor. Whereas in [26], the authors adopted a different strategy and introduced a class of meshes, where in the fine subinterval [0, τ ] the mesh is not necessarily equidistant but graded based on a mesh generating function. Then, for certain mesh generating functions the optimal estimate u − u N ≤ C N −1 follows. While for linear elements the influence of the logarithmic factor is not that large, the expressions ( N −1 lnN ) p and N − p differ significantly for larger p. Therefore, for higher order finite element methods it is important to use a graded mesh and not a piecewise equidistant mesh [25]. Consequently, this paper is aimed to introduce a modified graded mesh. This mesh is generated recursively using some implicitly defined functions. The fundamental idea is to incorporate the bounds on the convection term into the smooth grid presented in [25,32]. This term not only makes the mesh exponentially smaller but, it also adapts the layer more efficiently. Higher order parameter uniform convergence is obtained in -weighted energy norm. Moreover, similar to the Vulanovic´ and Bakhvalov-Shishkin meshes, the error estimates are optimal in the sense that they are free from logarithmic factor. A number of test examples are taken into account and a rigorous comparative analysis is presented. Test problems are solved using higher order finite element method over modified graded mesh presented in this paper. The results so obtained are compared with the results obtained over Shishkin mesh, mesh with a rational and optimal mesh characterizing function, Vulanovic´ mesh, Bakhvalov-Shishkin mesh, Gartland-Shishkin mesh (see [25]) and with modified Bakhvalov mesh [29]. The results obtained are also compared with B-spline methods [13,20]. 2. Statement of the problem and weak formulation Consider the following singularly perturbed boundary value problem
Lu := u (x) + a(x)u (x) + b(x)u (x) = f (x) in = (0, 1), u (0) = u (1) = 0
(2.1)
where 0 < 1 is a small positive parameter, a, b and f are sufficiently smooth functions. The solution of (2.1) typically has exponential boundary layers in the neighbourhood of x = 0. However, a simple transformation of variable x → 1 − x leads to the problem with layers along x = 1. In general, it may also exhibit interior layers, but for our analysis we shall assume that only exponential boundary and corner layers are present. Moreover, being interested in the case of dominating convection, it is standard to assume N ≤ C . Also assume that
a(x) ≥ α > 0 and
− b(x) + a (x)/2 > 0 for all x ∈ .
(2.2)
These assumptions guaranties the coercivity of the corresponding bilinear form. Consequently, the existence of a unique solution in H 01 () follows for all f ∈ L 2 (). Moreover, the solution of (2.1) admits a decomposition into regular and singular component [28, pp. 23] that reads
u = uR + uS
(2.3)
so that
(k) (k) u R (x) ≤ C and u S (x) ≤ C −k e −α x/
for all
x∈
(2.4)
and for prescribed non-negative integer k. The weak formulation of (2.1) can be interpreted as to determine u ∈ H 01 () such that
A (u , v ) = F ( v ) for all v ∈ H 01 (), (2.5) 2 where A (u , v ) = u , v − au , v − (bu , v ) and F ( v ) = − ( f , v ). Here (·, ·) is the standard scalar product in L -space. A standard norm associated with the bilinear form (2.5) is the energy norm . defined as
2 u 2 = u 0 + u 20 ;
u ∈ H 01 ()
where .0 is the standard L 2 -norm.
(2.6)
A. Kaushik et al. / Journal of Computational Physics 395 (2019) 275–285
277
3. Mesh generation N
Let = {t i : t i = i / N , i = 0, 1, 2 . . . N } be a uniform mesh. We try to find a non-uniform mesh = {xi }iN=0 on [0, 1]. Then there is an invertible monotone map : [0, 1] → [0, 1] such that (xi ) = t i . Thus, the non-uniform mesh
N
is the
N
preimage of under , and −1 might be observed as a deformation of into . We will assume that the non-negative function is continuously differentiable so that (x) = ψ(x) > 0 with (0) = 0 and (1) = 1. The functions and ψ might be perceived as distribution and density function. This terminology is used because we seek a suitable distribution of the grid points, whose density is inversely proportional to the local mesh width. Note that
dt = d = ψ dx so that dx = dt /ψ. In view of mean value theorem we may formally write
x i − x i −1 ≈
t i − t i −1 1 = . (xi −1 ) N ψ(xi −1 )
It is easy to follow that, the higher the density ψ(x), smaller is the mesh width. This implies that the mesh points are distributed densely. The solution of singular perturbation problem (2.1) is known to exhibit layer behaviour at x = 0. This suggests the choice of density function as
ψ(x) = μe −α x/τ + σ ,
x∈
(3.1)
where the constants μ, τ and σ are assumed to be positive and e −α x/τ is termed as layer correction function. Here, the normalization parameter, such that
μ is
x
(x) =
ψ(ξ )dξ
(3.2)
0
satisfies the condition (1) = 1. Intergrating and eliminating
(x) = σ x + (1 − σ )
1 − e −α x/τ 1 − e −α /τ
μ, we obtain the distribution function that reads (3.3)
.
Consequently, the corresponding density function reads
ψ(x) = σ +
α (1 − σ )e−α x/τ . τ (1 − e−α /τ )
(3.4)
Here, τ and σ are chosen such that τ > 0 and σ ∈ [0, 1]. The parameter τ regulates the thickness of boundary layer whereas σ standardize the density of mesh points. At σ = 0 all mesh points will effectively go to the layer region, while at σ = 1 the mesh is uniform. If we take σ = 0.5, approximately half of the mesh points will go to the layer region and the remaining half will go outside layer region. For numerical illustration we will take σ = 0.3 so that 70% of the mesh points will remain in the layer region. To construct the mesh we solve equation (xi ) = t i for xi using Newton’s method j +1
xi
j
j
= xi −
( x i ) − t i j
(3.5)
.
ψ(xi )
We start at x0 = 0 and successively compute xi for increasing value of i. The required mesh is then obtained iteratively with each iteration starting with x0i = xi −1 as the initial approximation. If h i denotes the successive mesh width then
hi =
τ λ − αNλ χi −1
:=
N
(3.6)
i
N (τ λσ + α (1 − σ )e ) α where λ = 1 − e − τ and χi is given by the recurrence relation
χ 0 = 0 , χ1 =
1
τ λσ + α (1 − σ )
and
χ i = χ i −1 +
1
τ λσ + α (1 − σ )e
for i = 2, 3, . . . , N − 1. An algorithm for mesh generation is as follows:
−α λ N
χ i −1
278
A. Kaushik et al. / Journal of Computational Physics 395 (2019) 275–285
Read ε , σ , τ , N, iter, and tol; Initialize mesh = 0; for i = 1 : N do x(i , 1) = mesh j = 2; while j ≤ iter do j −1
j
xi = xi if
j −1
−
( x i
) − ti
j −1
ψ(xi
j −1 |(xi ) − t i |
;
)
≤ tol then
mesh = x(i , j ); break;
end j = j + 1; end end mesh.
Algorithm 1: Mesh generation. 4. Interpolation error Let i = (xi −1 , xi ) be the ith element which is divided into p sub-elements as
x i −1 , j = x i −1 + j
hi p
;
j = 0, 1, . . . , p and i = 1, 2, . . . , N .
Suppose V pN be the finite element space composed of piecewise continuous Lagrange polynomials of degree p over . Further, φ0i , . . . , φ pi denotes Lagrange nodal shape functions associated with each i . The shape functions reads
φ ij (x) =
x − xi −1,k
0≤k≤ p
xi −1, j − xi −1,k
;
j = 0, 1, . . . , p and i = 1, 2, . . . , N .
k = j
The functions φ1i , . . . , φ pi −1 are the bubble functions [19] and the vertex functions are defined as
φi (x) =
φ pi (x) φ0i +1 (x)
if x ∈ i , if x ∈ i +1
for
i = 1, 2, . . . , N − 1.
These functions constitutes the basis of finite element space so that dim V pN = pN − 1. The shape functions vanishes at the boundary and therefore the space is contained in H 01 (). Let u I be the V pN -interpolant of u defined as
u I (x) =
N −1
u i φi (x) +
i =1
p −1 N
u ij φ ij (x),
x∈
(4.1)
i =1 j =1
where the coefficients u i and u ij are to be determined. The interpolant u I admits a decomposition into regular and singular component [19] given by
u I = u IR + u IS .
(4.2)
For each i , the Lagrange interpolation [12] leads to
u ( p +1 ) R and i ( p + 1)! i i ( p +1 ) (k) I u u − uS ≤ g (k) S S ( p + 1 )! i i (k) I u − uR R
≤ g (k)
(4.3a)
(4.3b)
i
where g (x) =
p
x − xi −1, j for x ∈ i and 0 ≤ k ≤ p. Hence, it follows that
j =0
p +1
g i ≤ h i
p and g ≤ ( p + 1)h i . i
(4.4)
A. Kaushik et al. / Journal of Computational Physics 395 (2019) 275–285
279
Consequently, (2.4), (4.3a) and (4.4) results in
(k) I u − uR R
≤ i
C
( p +1−k)
( p + 1 − k)!
hi
; k = 0, 1 .
( p +1)
The second inequality of (2.4) asserts that u S
i
(4.5)
vanishes for i that lies outside layer region. However, within the
layer region, inequalities (2.4), (4.3b) and (4.4) yields
(k) I u − uS S
≤ i
C −( p +1)
( p +1−k)
( p + 1 − k)!
hi
; k = 0, 1.
(4.6)
The next lemma provides the interpolation error bound in -weighted energy norm. Lemma 4.1. Let u I be the V pN -interpolant of the solution u and ε ≤ C 1 N −1 . Then
1 2 A1 B1 I + 2 N −p , u − u < C ε ( p !)2 N (( p + 1)!)2 for some constant C 1 with
A1 = B1 =
C1 N σ 2p
2p +2
+ 3+
1
σ 2 ( p +1 )
+
C1
N 2p +2
3C 1 N
2p +1 and (1 − σ ) (1−σ ) N
2p +3
+
(4.7)
C1
N 2p +3
(4.8)
2p +3 . (1 − σ ) (1−σ ) N
(4.9)
Proof. We may write from (2.6) that
2 2 I 2 I I u − u = ε u − u + u − u ε
0
0
:= I + J .
(4.10)
Now from (2.3) and (4.2) it follows that
2 I I ≤ ε u R (x) − u R (x) dx + ε
2 I u (x) − u S (x) dx S
+2ε u IR (x) − u R (x) u IS (x) − u S (x) dx.
(4.11)
Combining (3.6), (4.5) and (4.6) to discover
2p +1 2 I C2 C 1 (1−σ ) N C1 ε u R (x) − u R (x) dx ≤ + 2p N −(2p +1) . (1 − σ )C 1 N ( p !)2 σ
(4.12) Similarly, we obtain
2p +1 2 I C 2 (1 − σ ) (1−σ ) N ε u S (x) − u S (x) dx < N −2p ( p !)2
(4.13)
and
I 2p +1 C 2 (1 − σ ) I 2ε u R (x) − u R (x) u S (x) − u S (x) dx < 2
(1−σ ) N . ( p !)2 N 2p
(4.14)
Using (4.12), (4.13) and (4.14) into (4.11) to obtain
I <
C2
( p !)2
C1 N σ 2p
+ 3+
2p +2
C1
N 2p +2
(1 − σ ) (1−σ ) N
2p +1
N −2p .
(4.15)
280
A. Kaushik et al. / Journal of Computational Physics 395 (2019) 275–285
Again from (2.3) and (4.2) it follows that
2 2 I J ≤ u R (x) − u R (x) dx + u IS (x) − u S (x) dx
+2 u IR (x) − u R (x) u IS (x) − u S (x) dx.
(4.16)
Combining (3.6), (4.5) and (4.6) to obtain
2 I u R (x) − u R (x) dx ≤
C2
(( p + 1)!)2
(1 − σ )
C 1 (1−σ ) N
2p +3
N
+
1
σ 2 ( p +1 )
N −2 ( p +1 ) .
(4.17)
As in our earlier calculation, we compute
2 2p +3 −(2p +3) C 2 C 1 (1 − σ ) I
(1−σ ) N N u S (x) − u S (x) dx < 2 (( p + 1)!)
(4.18)
and that
2
I u R (x) − u R (x) u IS (x) − u S (x) dx < 2
C 2 C 1 (1 − σ )
(( p + 1)!)2 N 2p +3
2p +3
.
(4.19)
N −2 ( p +1 ) .
(4.20)
(1−σ ) N
Using (4.17), (4.18) and (4.19) into (4.16) to obtain
J<
C2
(( p + 1)!)2
1
σ 2 ( p +1 )
+
3C 1 N
2p +3
+
C1
N 2p +3
(1 − σ ) (1−σ ) N
Consequently from (4.15), (4.20) and (4.10) the required result follows.
2p +3
2
5. Discretization error N The discrete counterpart of the weak form (2.5) can be interpreted as to determine u N p ∈ V p such that
N A (u N p , v ) = F ( v ), for all v ∈ V p .
(5.1)
From (2.5) and (5.1) it follows that N A (u − u N p , v ) = 0, for all v ∈ V p .
(5.2)
The next two lemma provides the necessary results required to estimate the discretization error in -weighted energy norm. Lemma 5.1. Let u I be the V pN -interpolant of the exact solution u and u N p be the solution of (5.1). Then
C b 12 −( p +1) B 1 u I − u N A u I − u , u I − u Np < p N ( p + 1)!
(5.3)
for some constant C . N Proof. Set u I − u N p = v p . Then,
uI − u vN dx + a(x) u I − u vN A u I − u , v Np ≤ p p dx I + b(x) u − u v Np dx
:= K + L + M . Now from (2.3) and (4.2) it follows that
(5.4)
A. Kaushik et al. / Journal of Computational Physics 395 (2019) 275–285
K ≤ u IR − u R vN dx + u IS − u S vN dx . p p Since
vN p
(0) =
vN p
281
(5.5)
(1) = 0, we may formally compute
N C p u IR − u R vN dx ≤ hi vN dx = 0. p p i =1 p !
(5.6)
i
A similar calculation suggests that second term in (5.5) vanishes. Observe also that
u I − u dx = 0
L ≤ a(x) v N p
(5.7)
because u − u (0) = u − u (1) = 0. Using Cauchy-Schwarz inequality to obtain I
I
⎛
⎞1 ⎛ ⎞1 2 2 2 2 I N ⎝ ⎠ ⎝ ⎠ M ≤ b u − u dx v p dx
≤ b u I − u v Np . 0
(5.8)
Hence, the required result follows from (5.4)-(5.8) and (4.20).
2
The next lemma provides the -uniform error estimate in energy norm between the interpolant of exact solution and the solution of discrete weak form (5.1). Lemma 5.2. Let u I be the V pN -interpolant of the exact solution u and u N p be the solution of (5.1). Then
C b 12 −( p +1) I B N u − u Np < ε ( p + 1)! 1
(5.9)
for some constant C . Proof. Condition (2.2) implies the coercivity of the bilinear form. Also from (5.2) it is easy to follow that
2 I u − u Np ≤ A u I − u , u I − u Np .
(5.10)
ε
Consequently, Lemma (5.1) and (5.10) leads to the required result.
2
The next theorem which follows from triangle inequality, Lemma (4.1) and Lemma (5.2) presents the bound.
ε -uniform error
Theorem 5.3. Let u and u N p be the solutions of (2.1) and (5.1) respectively. Then
1 1 2 b A1 B1 N −1 2 + 2 + B N N −p u − u p < C ε ( p + 1)! 1 ( p !)2 N (( p + 1)!)2
(5.11)
for some constant C . 6. Numerical results and discussion In this section we take into account a number of test problems and present a rigorous comparative analysis. Test problems are solved using higher order finite element method over recursive modified graded mesh presented in this paper. The results so obtained are compared with the results obtained using higher order finite element method over Shishkin mesh, mesh with a rational and optimal mesh characterizing function, Bakhvalov-Shishkin mesh [19], Vulanovic´ mesh, GartlandShishkin mesh [25] and with modified Bakhvalov mesh [29]. In addition, results obtained are compared with B-spline method with artificial viscosity [13] and quintic B-spline method [20].
282
A. Kaushik et al. / Journal of Computational Physics 395 (2019) 275–285
Table 1 EN p , for Example 6.1 using Galerkin FEM with linear elements when
= 2−25 .
N
Bakhvalov-Shishkin mesh [19,25]
Gartland-Shishkin mesh [25]
Modified graded mesh
20 41 81 162 322 642
1.162(−2) 8.508(−3) 4.306(−3) 2.140(−3) 1.080(−3) 5.425(−4)
1.083(−2) 4.523(−3) 2.123(−3) 1.048(−3) 5.234(−4) 2.618(−4)
4.773(−4) 1.215(−4) 3.181(−5) 8.020(−6) 2.209(−6) 9.827(−7)
Table 2 EN p , and O rd( N ) for Example 6.1 using Galerkin FEM with quadratic elements when N
Shishkin mesh [25]
Vulanovic´ mesh [25]
Bakhvalov-Shishkin mesh [25]
Modified graded mesh
32
2.534(−3) 1.454 9.251(−4) 1.548 3.164(−4) 1.612 1.035(−4) 1.659 3.277(−5) 1.696 1.012(−5) 1.725 3.060(−6) –
3.818(−4) 1.902 1.022(−4) 1.935 2.674(−5) 1.946 6.939(−6) 1.953 1.792(−6) 1.960 4.605(−7) 1.967 1.178(−7) –
3.607(−4) 1.963 9.253(−5) 1.995 2.322(−5) 1.999 5.807(−6) 1.998 1.454(−6) 1.998 3.640(−7) 1.999 9.107(−8) –
1.249e−4 2.110 2.894e−5 2.435 5.351e−6 2.785 7.766e−7 2.936 1.015e−7 2.982 1.284e−8 2.995 1.610e−9 –
64 128 256 512 1024 2048
Suppose E N p , and O rd( N ) denotes the error in the convergence respectively. Moreover, we define
N EN p , = u − u p
= 2−12 .
and O rd( N ) =
-weighted energy norm and the corresponding numerical rate of
2N ln( E N p , ) − ln( E p , )
ln(2)
.
N Further, e N p , = u − u p N defines the maximum absolute error. In [19], results are presented for linear elements over Shishkin mesh, Bakhvalov-Shishkin mesh and for a graded mesh having a rational mesh generating function. Whereas, in [25] Galerkin FEM with quadratic and cubic elements are brought to attention. They consider five different S-type meshes; the original Shishkin mesh, a mesh with a rational mesh characterizing function, the Vulanovic´ mesh, the Bakhvalov-Shishkin mesh and the mesh with optimal mesh characterizing function ψopt . Moreover, a modified Bakhvalov mesh is the subject matter of [29]. In this section we present and compare results obtained over modified graded mesh with other adaptive/graded meshes developed over the last decades. Numerical results for Example 6.1 are conferred in Tables 1–3 and for Example 6.2 in Table 4. For linear elements, EN p , is tabulated and compared with the error obtained over Bakhvalov-Shishkin and Gartland-Shishkin meshes in the
form of Table 1. Moreover, for quadratic elements E N p , is tabulated and compared with the error obtained over Shishkin mesh, Vulanovic´ mesh, Bakhvalov-Shishkin mesh and Gartland-Shishkin mesh in the form of Table 2. For different values of perturbation parameter e N p , is tabulated and compared with an upwind scheme defined over a modified Bakhvalov mesh [29] and with a B-spline method with artificial viscosity [13]. The details of which can be had from Tables 3–4. Maximum absolute error (e N p , ) obtained using Galerkin FEM over the mesh presented is compared with a number of adaptive/graded meshes developed over the last decades (see Fig. 1). For different values of the method presented is compared with quintic B-spline method [20] (see Fig. 2). Example 6.1. Consider the following singular perturbation problem from [19,29,25]:
u + u − 2u = −e x−1 in = (0, 1), u (0) = u (1) = 0. Example 6.2. Consider the following singular perturbation problem from [7,20]:
u (x) + u (x) − u (x) = 0 in = (0, 1), u (0) = u (1) = 1.
A. Kaushik et al. / Journal of Computational Physics 395 (2019) 275–285
283
Table 3 Comparison of maximum absolute error (e N p , ) for Example 6.1 with an upwind scheme defined over a modified Bakhvalov mesh [29].
eN p ,
32
64
128
256
512
1024
2048
2−10
eN p ,
8.696(−4)
2.167(−4)
5.407(−5)
1.352(−5)
3.379(−6)
8.447(−7)
2.112(−7)
[29]
2.110(−2)
1.108(−2)
5.543(−3)
2.714(−3)
1.356(−3)
6.778(−4)
3.388(−4)
2−12
eN p ,
8.779(−4)
2.192(−4)
5.467(−5)
1.367(−5)
3.417(−6)
8.542(−7)
2.136(−7)
[29]
2.163(−2)
1.104(−2)
5.539(−3)
2.765(−3)
1.379(−3)
6.884(−4)
3.439(−4)
2−14
eN p ,
8.799(−4)
2.200(−4)
5.483(−5)
1.371(−5)
3.426(−6)
8.563(−7)
2.155(−7)
[29]
2.189(−2)
1.114(−2)
5.580(−3)
2.782(−3)
1.387(−3)
6.921(−4)
3.455(−4)
2−16
eN p ,
8.803(−4)
2.202(−4)
5.485(−5)
1.371(−5)
3.430(−6)
8.566(−7)
2.172(−7)
[29]
2.205(−2)
1.120(−2)
5.560(−3)
2.789(−3)
1.390(−3)
6.932(−4)
3.461(−4)
2−18
eN p ,
8.805(−4)
2.203(−4)
5.485(−5)
1.371(−5)
3.427(−6)
8.611(−7)
1.700(−7)
[29]
2.215(−2)
1.123(−2)
5.611(−3)
2.793(−3)
1.391(−3)
6.938(−4)
3.463(−4)
eN p ,
8.805(−4)
2.203(−4)
5.484(−5)
1.371(−5)
3.417(−6)
8.724(−7)
2.860(−7)
[29]
2.223(−2)
1.126(−2)
5.618(−3)
2.796(−3)
1.392(−3)
6.940(−4)
3.464(−4)
eN p ,
8.805(−4)
2.203(−4)
5.485(−5)
1.371(−5)
3.399(−6)
8.747(−7)
5.527(−7)
[29]
2.223(−2)
1.128(−2)
5.624(−3)
2.797(−3)
1.392(−3)
6.941(−4)
3.464(−4)
eN p ,
8.805(−4)
2.203(−4)
5.482(−5)
1.369(−5)
3.198(−6)
1.232(−6)
3.077(−7)
[29]
2.223(−2)
1.128(−2)
5.628(−3)
2.799(−3)
1.392(−3)
6.942(−4)
3.464(−4)
512
1024
2−20 2−22 2−24
Table 4 Comparison of maximum absolute error (e N p , ) for Example 6.2 with a B-spline method with artificial viscosity [13].
eN p ,
16
32
2−12
eN p ,
4.182(−3)
[13]
1.085(−2)
eN p ,
4.185(−3)
1.001(−3)
2.483(−4)
6.218(−5)
1.569(−5)
3.922(−6)
9.804(−7)
[13]
1.089(−2)
5.574(−3)
2.798(−3)
1.385(−3)
6.718(−4)
3.139(−4)
1.347(−4)
eN p ,
4.186(−3)
1.000(−3)
2.484(−4)
6.165(−5)
1.573(−5)
3.923(−6)
9.809(−7)
[13]
1.091(−2)
5.597(−3)
2.821(−3)
1.407(−3)
6.942(−4)
3.363(−4)
1.571(−4)
eN p ,
4.187(−3)
9.994(−4)
2.485(−4)
6.127(−5)
1.582(−5)
3.924(−6)
9.811(−7)
[13]
1.092(−2)
5.608(−3)
2.832(−3)
1.418(−3)
7.054(−4)
3.476(−4)
1.683(−4)
eN p ,
4.187(−3)
9.990(−4)
2.486(−4)
6.103(−5)
1.592(−5)
3.925(−6)
9.812(−7)
[13]
1.093(−2)
5.614(−3)
2.838(−3)
1.424(−3)
7.110(−4)
3.532(−4)
1.739(−4)
eN p ,
4.188(−3)
9.986(−4)
2.487(−4)
6.077(−5)
1.607(−5)
3.894(−6)
9.755(−7)
[13]
1.093(−2)
5.619(−3)
2.843(−3)
1.429(−3)
7.163(−4)
3.584(−4)
1.792(−4)
eN p ,
4.188(−3)
9.986(−4)
2.487(−4)
6.075(−5)
1.609(−5)
3.886(−6)
9.716(−7)
[13]
1.093(−2)
5.619(−3)
2.843(−3)
1.429(−3)
7.166(−4)
3.588(−4)
1.795(−4)
2−13 2−14 2−15 2−16 2−20 2−25
64
128
256
1.004(−3)
2.483(−4)
6.262(−5)
1.567(−5)
3.918(−6)
9.796(−7)
5.530(−3)
2.754(−3)
1.340(−3)
6.269(−4)
2.692(−4)
9.640(−5)
Fig. 1. Error (E N p , ) plot for Example 6.1, using Galerkin FEM with quadratic elements, over a mesh with a rational and optimal mesh characterizing function [25] and the modified graded mesh presented.
284
A. Kaushik et al. / Journal of Computational Physics 395 (2019) 275–285
Fig. 2. Comparison of maximum absolute error (e N p , ) for Example 6.2 with quintic B-spline method [20].
A comparison through numerical experiments reveals the dominance of the modified graded mesh over the state of the art meshes. Moreover, the results obtained are better than the results of B-spline method with artificial viscosity [13] and quintic B-spline method [20]. 7. Conclusion A singularly perturbed convection-diffusion problem is solved numerically using FEM based on higher order polynomials. A modified graded mesh is generated using some defined functions. Higher order convergence is obtained in implicitly -weighted energy norm. A convergence with O N − p is established. The error estimates obtained are optimal in the sense that they are free from logarithmic factor. Test examples are taken into account and a rigorous comparative analysis is presented. Moreover, we compare the proposed method and graded mesh generated with others found in the literature. References [1] N.S. Bakhvalov, On the optimization of methods for solving boundary value problems with boundary layers, Ž. Vyˇcisl. Mat. Mat. Fiz. 9 (1969) 841–859. ¸ Singularly Perturbed Boundary Value Problems, Birkhäusar, Basel, 2007. [2] L. Barbu, G. Morosanu, [3] M.G. Beckett, J.A. Mackenzie, On a uniformly accurate finite difference approximation of a singularly perturbed reaction diffusion problem using grid equidistribution, J. Comput. Appl. Math. 131 (2001) 381–405. [4] M. Brdar, H. Zarin, On graded meshes for a two parameter singularly perturbed problem, Appl. Math. Comput. 282 (2016) 97–107. [5] M. Brdar, H. Zarin, A singularly perturbed problem with two parameters on a Bakhvalov type mesh, J. Comput. Appl. Math. 292 (2016) 307–319. [6] M. Brdar, H. Zarin, L. Teofanov, A singularly perturbed problem with two parameters in two dimensions on graded meshes, Comput. Math. Appl. 72 (2016) 2582–2603. [7] C. Clavero, J.C. Jorge, F. Lisbona, G.I. Shishkin, An alternating direction scheme on a nonuniform mesh for reaction diffusion parabolic problem, IMA J. Numer. Anal. 20 (2000) 263–280. [8] R.G. Durán, A.L. Lombardi, Finite element approximation of convection diffusion problems using graded meshes, Appl. Numer. Math. 56 (2006) 1314–1325. [9] R.G. Durán, A.L. Lombardi, M.I. Prieto, Superconvergence for finite element approximation of a convection diffusion equation using graded meshes, IMA J. Numer. Anal. 32 (2012) 511–533. [10] Z. Gajic, M.T. Lim, Optimal Control of Singularly Perturbed Linear Systems and Applications, Marcel-Dekker, New York, 2001. [11] E.C. Gartland, Graded mesh difference schemes for singularly perturbed two point boundary value problems, Math. Comput. 51 (1988) 631–657. [12] G.W. Howell, Derivative error bounds for Lagrange interpolation: an extension of Cauchy’s bound for the error of Lagrange interpolation, J. Approx. Theory 67 (1991) 164–173. [13] M.K. Kadalbajoo, P. Arora, B-spline with artificial viscosity for solving singularly perturbed boundary value problems, Math. Comput. Model. 52 (2010) 654–666. [14] A. Kaushik, Singular perturbation analysis of bistable differential equation arising in the nerve pulse propagation, Nonlinear Anal. 9 (2008) 2106–2127. [15] A. Kaushik, M.D. Sharma, Numerical analysis of a mathematical model for propagation of an electrical pulse in a neuron, Numer. Methods Partial Differ. Equ. 27 (2008) 1–18. [16] N. Kopteva, N. Madden, M. Stynes, Grid equidistribution for reaction diffusion problems in one dimension, Numer. Algorithms 40 (2005) 305–322. [17] N. Kopteva, E. O’Riordan, Shishkin meshes in the numerical solution of singularly perturbed differential equations, Int. J. Numer. Anal. Model. 7 (2010) 393–415. [18] T. Linß, Analysis of a Galerkin finite element method on a Bakhvalov-Shishkin mesh for a linear convection diffusion problem, IMA J. Numer. Anal. 20 (2000) 621–632. [19] T. Linß, Layer Adapted Meshes for Reaction Convection Diffusion Problems, Springer, Berlin, 2010.
A. Kaushik et al. / Journal of Computational Physics 395 (2019) 275–285
285
[20] R.K. Lodhi, H.K. Mishra, Quintic B-spline method for solving second order linear and nonlinear singularly perturbed two point boundary value problems, J. Comput. Appl. Math. 319 (2017) 170–187. [21] J.J.H. Miller, Singular Perturbation Problems in Chemical Physics, John Wiley & Sons, New York, 1997. [22] J. Mohapatra, S. Natesan, Uniformly convergent numerical method for singularly perturbed differential difference equation using grid equidistribution, Int. J. Numer. Methods Biomed. Eng. 27 (2011) 1427–2445. [23] R.E. O’Malley, Introduction to Singular Perturbations, Academic Press, New York, 1974. [24] H.-G. Roos, Error estimates for linear finite elements on Bakhvalov type meshes, Appl. Math. 51 (2006) 63–72. [25] H.-G. Roos, Graded meshes for higher order FEM, J. Comput. Math. 33 (2015) 1–16. [26] H.-G. Roos, T. Linß, Sufficient conditions for uniform convergence on layer adapted grids, Computing 63 (1999) 27–45. [27] H.-G. Roos, T. Skalický, A comparison of the finite element method on Shishkin and Garland type meshes for convection diffusion problems, CWI Q. 10 (1997) 277–300. [28] H.-G. Roos, M. Stynes, L. Tobiska, Robust Numerical Methods for Singularly Perturbed Differential Equations, Springer, Berlin, 2008. [29] H.-G. Roos, L. Teofanov, Z. Uzelac, A modified Bakhvalov mesh, Appl. Math. Lett. 31 (2014) 7–11. [30] H. Schlichting, K. Gersten, Boundary Layer Theory, Springer-Verlag, Berlin, 2003. [31] K.K. Sharma, A. Kaushik, A solution of the discrepancy occurs due to using the fitted mesh approach rather than to the fitted operator for solving singularly perturbed differential equations, Appl. Math. Comput. 181 (2006) 756–766. [32] G. Söderlind, A.S. Yadaw, The impact of smooth W -grids in the numerical solution of singular perturbation two point boundary value problems, Appl. Math. Comput. 218 (2012) 6045–6055. [33] G. Sun, M. Stynes, Finite element methods for singularly perturbed higher order elliptic two point boundary value problems ii: convection diffusion type, IMA J. Numer. Anal. 15 (1995) 197–219. ´ A priori meshes for singularly perturbed quasilinear two point boundary value problems, IMA J. Numer. Anal. 21 (2001) 349–366. [34] R. Vulanovic, ´ L. Teofanov, A modification of the Shishkin discretization mesh for one dimensional reaction diffusion problems, Appl. Math. Comput. [35] R. Vulanovic, 220 (2013) 104–116. [36] C. Xenophontos, S. Franz, L. Ludwig, Finite element approximation of convection diffusion problems using an exponentially graded mesh, Comput. Math. Appl. 72 (2016) 1532–1540. [37] Q. Zheng, X. Zheng, Y.F. Liu, Uniform second order hybrid schemes on Bakhvalov-Shishkin mesh for quasilinear convection diffusion problems, Adv. Mater. Res. 871 (2014) 135–140.
Contents lists available at ScienceDirect
Journal of Computational Physics www.elsevier.com/locate/jcp
A modified graded mesh and higher order finite element approximation for singular perturbation problems Aditya Kaushik a,c , Anil K. Vashishth b , Vijayant Kumar b , Manju Sharma d,∗ a
Department of Applied Mathematics, Delhi Technological University, Delhi – 110042, India Department of Mathematics, Kurukshetra University, Kurukshetra – 136119, India University Institute of Engineering and Technology, Panjab University, Chandigarh – 160014, India d Department of Mathematics, KVA DAV College, Karnal – 132001, India b c
a r t i c l e
i n f o
Article history: Received 20 February 2018 Received in revised form 31 March 2019 Accepted 11 April 2019 Available online 12 June 2019 Keywords: Singular perturbation Boundary value problem Layer-adapted meshes Finite element method
a b s t r a c t A singularly perturbed convection diffusion problem is solved numerically using finite element method based on higher order polynomials. More precisely, we introduce a modified graded mesh generated using some implicitly defined functions. Higher order parameter uniform convergence is obtained in -weighted energy norm. Moreover, the error estimates obtained are optimal in the sense that they are free from logarithmic factor. A number of test examples are taken into account and a rigorous comparative analysis is presented. Moreover, we compare the proposed method with others found in the literature. © 2019 Elsevier Inc. All rights reserved.
1. Introduction It is well known that many phenomena in biology [14,15], chemistry [21], engineering [30,10], physics [23] can be described by boundary value problems associated with various types of ordinary and partial differential equations or systems. When we associate a mathematical model with a phenomenon, we generally try to capture what is essential, retaining the important quantities and omitting the negligible ones which involve small parameters. The model that would be obtained by maintaining the small parameters is called the perturbed model, whereas the degenerate simplified model is called unperturbed or reduced model [2]. Of course, the unperturbed model is to be preferred, because it is simpler. What matters is that it should describe faithfully enough the respective phenomenon, which means that its solution must be sufficiently close enough to the solution of the corresponding perturbed model [2]. This fact holds in the case of regular perturbations. On the other hand, in the case of singular perturbations, things get more complicated. If we refer to an initial boundary value problem, the solution of the unperturbed problem does not satisfy in general all the original boundary conditions and/or initial conditions as some of the derivatives may disappear by neglecting the small parameters [2,28]. Thus, some discrepancy may appear between the solution of the perturbed model and that of the corresponding reduced model [31]. Numerical discretization of singularly perturbed differential equation is burdened with difficulties; for example, finite difference method (FDM) and finite element method (FEM) necessitate a mesh to sustain approximations [19,28]. In fact, the discrete solutions obtained using standard Galerkin or cantered finite differences methods demonstrate oscillatory behaviour for small discretization parameters. However, layer adapted meshes seem promising in the discretization of such equations [19]. This leads to growing interest in the generation of layer-adapted meshes. Layer adapted meshes have first been pro-
*
Corresponding author. E-mail addresses: [email protected] (A. Kaushik), [email protected] (A.K. Vashishth), [email protected] (M. Sharma).
https://doi.org/10.1016/j.jcp.2019.04.073 0021-9991/© 2019 Elsevier Inc. All rights reserved.
276
A. Kaushik et al. / Journal of Computational Physics 395 (2019) 275–285
posed by Bakhvalov [1,24,5] in the context of reaction-diffusion problems. Later on, special meshes for convection-diffusion problems has been extensively investigated. That includes Shishkin’s piecewise equidistant meshes [19,17], grids that are based on equidistribution [16,3,22], Gartland type meshes [11], Bakhvalov-Shishkin grids [18,37] and Vulanovic´ improved Shishkin grids [34,35]. Further development leads to the layer-adapted meshes defined using a recursive formula. Examples are the Gartland-Shishkin meshes [27] or the graded meshes analysed in [9,8,6,4,36] and others. The first parameter uniform convergence result, for linear finite elements on a Shishkin mesh, in the -weighted H 1 -norm was proved in [33]. The given interval is divided into two subintervals that includes a small subinterval [0, τ ] and a large subinterval [τ , 1]. Here, the choice of τ = τ0 ln N is important. However, the mesh is equidistant in both subintervals but, fine in [0, τ ] and coarse in [τ , 1] with a mesh size of order O ( N −1 ). Moreover, the resulting error estimate u − u N ≤ C N −1 ln N is not optimal due to the presence of logarithmic factor. Whereas in [26], the authors adopted a different strategy and introduced a class of meshes, where in the fine subinterval [0, τ ] the mesh is not necessarily equidistant but graded based on a mesh generating function. Then, for certain mesh generating functions the optimal estimate u − u N ≤ C N −1 follows. While for linear elements the influence of the logarithmic factor is not that large, the expressions ( N −1 lnN ) p and N − p differ significantly for larger p. Therefore, for higher order finite element methods it is important to use a graded mesh and not a piecewise equidistant mesh [25]. Consequently, this paper is aimed to introduce a modified graded mesh. This mesh is generated recursively using some implicitly defined functions. The fundamental idea is to incorporate the bounds on the convection term into the smooth grid presented in [25,32]. This term not only makes the mesh exponentially smaller but, it also adapts the layer more efficiently. Higher order parameter uniform convergence is obtained in -weighted energy norm. Moreover, similar to the Vulanovic´ and Bakhvalov-Shishkin meshes, the error estimates are optimal in the sense that they are free from logarithmic factor. A number of test examples are taken into account and a rigorous comparative analysis is presented. Test problems are solved using higher order finite element method over modified graded mesh presented in this paper. The results so obtained are compared with the results obtained over Shishkin mesh, mesh with a rational and optimal mesh characterizing function, Vulanovic´ mesh, Bakhvalov-Shishkin mesh, Gartland-Shishkin mesh (see [25]) and with modified Bakhvalov mesh [29]. The results obtained are also compared with B-spline methods [13,20]. 2. Statement of the problem and weak formulation Consider the following singularly perturbed boundary value problem
Lu := u (x) + a(x)u (x) + b(x)u (x) = f (x) in = (0, 1), u (0) = u (1) = 0
(2.1)
where 0 < 1 is a small positive parameter, a, b and f are sufficiently smooth functions. The solution of (2.1) typically has exponential boundary layers in the neighbourhood of x = 0. However, a simple transformation of variable x → 1 − x leads to the problem with layers along x = 1. In general, it may also exhibit interior layers, but for our analysis we shall assume that only exponential boundary and corner layers are present. Moreover, being interested in the case of dominating convection, it is standard to assume N ≤ C . Also assume that
a(x) ≥ α > 0 and
− b(x) + a (x)/2 > 0 for all x ∈ .
(2.2)
These assumptions guaranties the coercivity of the corresponding bilinear form. Consequently, the existence of a unique solution in H 01 () follows for all f ∈ L 2 (). Moreover, the solution of (2.1) admits a decomposition into regular and singular component [28, pp. 23] that reads
u = uR + uS
(2.3)
so that
(k) (k) u R (x) ≤ C and u S (x) ≤ C −k e −α x/
for all
x∈
(2.4)
and for prescribed non-negative integer k. The weak formulation of (2.1) can be interpreted as to determine u ∈ H 01 () such that
A (u , v ) = F ( v ) for all v ∈ H 01 (), (2.5) 2 where A (u , v ) = u , v − au , v − (bu , v ) and F ( v ) = − ( f , v ). Here (·, ·) is the standard scalar product in L -space. A standard norm associated with the bilinear form (2.5) is the energy norm . defined as
2 u 2 = u 0 + u 20 ;
u ∈ H 01 ()
where .0 is the standard L 2 -norm.
(2.6)
A. Kaushik et al. / Journal of Computational Physics 395 (2019) 275–285
277
3. Mesh generation N
Let = {t i : t i = i / N , i = 0, 1, 2 . . . N } be a uniform mesh. We try to find a non-uniform mesh = {xi }iN=0 on [0, 1]. Then there is an invertible monotone map : [0, 1] → [0, 1] such that (xi ) = t i . Thus, the non-uniform mesh
N
is the
N
preimage of under , and −1 might be observed as a deformation of into . We will assume that the non-negative function is continuously differentiable so that (x) = ψ(x) > 0 with (0) = 0 and (1) = 1. The functions and ψ might be perceived as distribution and density function. This terminology is used because we seek a suitable distribution of the grid points, whose density is inversely proportional to the local mesh width. Note that
dt = d = ψ dx so that dx = dt /ψ. In view of mean value theorem we may formally write
x i − x i −1 ≈
t i − t i −1 1 = . (xi −1 ) N ψ(xi −1 )
It is easy to follow that, the higher the density ψ(x), smaller is the mesh width. This implies that the mesh points are distributed densely. The solution of singular perturbation problem (2.1) is known to exhibit layer behaviour at x = 0. This suggests the choice of density function as
ψ(x) = μe −α x/τ + σ ,
x∈
(3.1)
where the constants μ, τ and σ are assumed to be positive and e −α x/τ is termed as layer correction function. Here, the normalization parameter, such that
μ is
x
(x) =
ψ(ξ )dξ
(3.2)
0
satisfies the condition (1) = 1. Intergrating and eliminating
(x) = σ x + (1 − σ )
1 − e −α x/τ 1 − e −α /τ
μ, we obtain the distribution function that reads (3.3)
.
Consequently, the corresponding density function reads
ψ(x) = σ +
α (1 − σ )e−α x/τ . τ (1 − e−α /τ )
(3.4)
Here, τ and σ are chosen such that τ > 0 and σ ∈ [0, 1]. The parameter τ regulates the thickness of boundary layer whereas σ standardize the density of mesh points. At σ = 0 all mesh points will effectively go to the layer region, while at σ = 1 the mesh is uniform. If we take σ = 0.5, approximately half of the mesh points will go to the layer region and the remaining half will go outside layer region. For numerical illustration we will take σ = 0.3 so that 70% of the mesh points will remain in the layer region. To construct the mesh we solve equation (xi ) = t i for xi using Newton’s method j +1
xi
j
j
= xi −
( x i ) − t i j
(3.5)
.
ψ(xi )
We start at x0 = 0 and successively compute xi for increasing value of i. The required mesh is then obtained iteratively with each iteration starting with x0i = xi −1 as the initial approximation. If h i denotes the successive mesh width then
hi =
τ λ − αNλ χi −1
:=
N
(3.6)
i
N (τ λσ + α (1 − σ )e ) α where λ = 1 − e − τ and χi is given by the recurrence relation
χ 0 = 0 , χ1 =
1
τ λσ + α (1 − σ )
and
χ i = χ i −1 +
1
τ λσ + α (1 − σ )e
for i = 2, 3, . . . , N − 1. An algorithm for mesh generation is as follows:
−α λ N
χ i −1
278
A. Kaushik et al. / Journal of Computational Physics 395 (2019) 275–285
Read ε , σ , τ , N, iter, and tol; Initialize mesh = 0; for i = 1 : N do x(i , 1) = mesh j = 2; while j ≤ iter do j −1
j
xi = xi if
j −1
−
( x i
) − ti
j −1
ψ(xi
j −1 |(xi ) − t i |
;
)
≤ tol then
mesh = x(i , j ); break;
end j = j + 1; end end mesh.
Algorithm 1: Mesh generation. 4. Interpolation error Let i = (xi −1 , xi ) be the ith element which is divided into p sub-elements as
x i −1 , j = x i −1 + j
hi p
;
j = 0, 1, . . . , p and i = 1, 2, . . . , N .
Suppose V pN be the finite element space composed of piecewise continuous Lagrange polynomials of degree p over . Further, φ0i , . . . , φ pi denotes Lagrange nodal shape functions associated with each i . The shape functions reads
φ ij (x) =
x − xi −1,k
0≤k≤ p
xi −1, j − xi −1,k
;
j = 0, 1, . . . , p and i = 1, 2, . . . , N .
k = j
The functions φ1i , . . . , φ pi −1 are the bubble functions [19] and the vertex functions are defined as
φi (x) =
φ pi (x) φ0i +1 (x)
if x ∈ i , if x ∈ i +1
for
i = 1, 2, . . . , N − 1.
These functions constitutes the basis of finite element space so that dim V pN = pN − 1. The shape functions vanishes at the boundary and therefore the space is contained in H 01 (). Let u I be the V pN -interpolant of u defined as
u I (x) =
N −1
u i φi (x) +
i =1
p −1 N
u ij φ ij (x),
x∈
(4.1)
i =1 j =1
where the coefficients u i and u ij are to be determined. The interpolant u I admits a decomposition into regular and singular component [19] given by
u I = u IR + u IS .
(4.2)
For each i , the Lagrange interpolation [12] leads to
u ( p +1 ) R and i ( p + 1)! i i ( p +1 ) (k) I u u − uS ≤ g (k) S S ( p + 1 )! i i (k) I u − uR R
≤ g (k)
(4.3a)
(4.3b)
i
where g (x) =
p
x − xi −1, j for x ∈ i and 0 ≤ k ≤ p. Hence, it follows that
j =0
p +1
g i ≤ h i
p and g ≤ ( p + 1)h i . i
(4.4)
A. Kaushik et al. / Journal of Computational Physics 395 (2019) 275–285
279
Consequently, (2.4), (4.3a) and (4.4) results in
(k) I u − uR R
≤ i
C
( p +1−k)
( p + 1 − k)!
hi
; k = 0, 1 .
( p +1)
The second inequality of (2.4) asserts that u S
i
(4.5)
vanishes for i that lies outside layer region. However, within the
layer region, inequalities (2.4), (4.3b) and (4.4) yields
(k) I u − uS S
≤ i
C −( p +1)
( p +1−k)
( p + 1 − k)!
hi
; k = 0, 1.
(4.6)
The next lemma provides the interpolation error bound in -weighted energy norm. Lemma 4.1. Let u I be the V pN -interpolant of the solution u and ε ≤ C 1 N −1 . Then
1 2 A1 B1 I + 2 N −p , u − u < C ε ( p !)2 N (( p + 1)!)2 for some constant C 1 with
A1 = B1 =
C1 N σ 2p
2p +2
+ 3+
1
σ 2 ( p +1 )
+
C1
N 2p +2
3C 1 N
2p +1 and (1 − σ ) (1−σ ) N
2p +3
+
(4.7)
C1
N 2p +3
(4.8)
2p +3 . (1 − σ ) (1−σ ) N
(4.9)
Proof. We may write from (2.6) that
2 2 I 2 I I u − u = ε u − u + u − u ε
0
0
:= I + J .
(4.10)
Now from (2.3) and (4.2) it follows that
2 I I ≤ ε u R (x) − u R (x) dx + ε
2 I u (x) − u S (x) dx S
+2ε u IR (x) − u R (x) u IS (x) − u S (x) dx.
(4.11)
Combining (3.6), (4.5) and (4.6) to discover
2p +1 2 I C2 C 1 (1−σ ) N C1 ε u R (x) − u R (x) dx ≤ + 2p N −(2p +1) . (1 − σ )C 1 N ( p !)2 σ
(4.12) Similarly, we obtain
2p +1 2 I C 2 (1 − σ ) (1−σ ) N ε u S (x) − u S (x) dx < N −2p ( p !)2
(4.13)
and
I 2p +1 C 2 (1 − σ ) I 2ε u R (x) − u R (x) u S (x) − u S (x) dx < 2
(1−σ ) N . ( p !)2 N 2p
(4.14)
Using (4.12), (4.13) and (4.14) into (4.11) to obtain
I <
C2
( p !)2
C1 N σ 2p
+ 3+
2p +2
C1
N 2p +2
(1 − σ ) (1−σ ) N
2p +1
N −2p .
(4.15)
280
A. Kaushik et al. / Journal of Computational Physics 395 (2019) 275–285
Again from (2.3) and (4.2) it follows that
2 2 I J ≤ u R (x) − u R (x) dx + u IS (x) − u S (x) dx
+2 u IR (x) − u R (x) u IS (x) − u S (x) dx.
(4.16)
Combining (3.6), (4.5) and (4.6) to obtain
2 I u R (x) − u R (x) dx ≤
C2
(( p + 1)!)2
(1 − σ )
C 1 (1−σ ) N
2p +3
N
+
1
σ 2 ( p +1 )
N −2 ( p +1 ) .
(4.17)
As in our earlier calculation, we compute
2 2p +3 −(2p +3) C 2 C 1 (1 − σ ) I
(1−σ ) N N u S (x) − u S (x) dx < 2 (( p + 1)!)
(4.18)
and that
2
I u R (x) − u R (x) u IS (x) − u S (x) dx < 2
C 2 C 1 (1 − σ )
(( p + 1)!)2 N 2p +3
2p +3
.
(4.19)
N −2 ( p +1 ) .
(4.20)
(1−σ ) N
Using (4.17), (4.18) and (4.19) into (4.16) to obtain
J<
C2
(( p + 1)!)2
1
σ 2 ( p +1 )
+
3C 1 N
2p +3
+
C1
N 2p +3
(1 − σ ) (1−σ ) N
Consequently from (4.15), (4.20) and (4.10) the required result follows.
2p +3
2
5. Discretization error N The discrete counterpart of the weak form (2.5) can be interpreted as to determine u N p ∈ V p such that
N A (u N p , v ) = F ( v ), for all v ∈ V p .
(5.1)
From (2.5) and (5.1) it follows that N A (u − u N p , v ) = 0, for all v ∈ V p .
(5.2)
The next two lemma provides the necessary results required to estimate the discretization error in -weighted energy norm. Lemma 5.1. Let u I be the V pN -interpolant of the exact solution u and u N p be the solution of (5.1). Then
C b 12 −( p +1) B 1 u I − u N A u I − u , u I − u Np < p N ( p + 1)!
(5.3)
for some constant C . N Proof. Set u I − u N p = v p . Then,
uI − u vN dx + a(x) u I − u vN A u I − u , v Np ≤ p p dx I + b(x) u − u v Np dx
:= K + L + M . Now from (2.3) and (4.2) it follows that
(5.4)
A. Kaushik et al. / Journal of Computational Physics 395 (2019) 275–285
K ≤ u IR − u R vN dx + u IS − u S vN dx . p p Since
vN p
(0) =
vN p
281
(5.5)
(1) = 0, we may formally compute
N C p u IR − u R vN dx ≤ hi vN dx = 0. p p i =1 p !
(5.6)
i
A similar calculation suggests that second term in (5.5) vanishes. Observe also that
u I − u dx = 0
L ≤ a(x) v N p
(5.7)
because u − u (0) = u − u (1) = 0. Using Cauchy-Schwarz inequality to obtain I
I
⎛
⎞1 ⎛ ⎞1 2 2 2 2 I N ⎝ ⎠ ⎝ ⎠ M ≤ b u − u dx v p dx
≤ b u I − u v Np . 0
(5.8)
Hence, the required result follows from (5.4)-(5.8) and (4.20).
2
The next lemma provides the -uniform error estimate in energy norm between the interpolant of exact solution and the solution of discrete weak form (5.1). Lemma 5.2. Let u I be the V pN -interpolant of the exact solution u and u N p be the solution of (5.1). Then
C b 12 −( p +1) I B N u − u Np < ε ( p + 1)! 1
(5.9)
for some constant C . Proof. Condition (2.2) implies the coercivity of the bilinear form. Also from (5.2) it is easy to follow that
2 I u − u Np ≤ A u I − u , u I − u Np .
(5.10)
ε
Consequently, Lemma (5.1) and (5.10) leads to the required result.
2
The next theorem which follows from triangle inequality, Lemma (4.1) and Lemma (5.2) presents the bound.
ε -uniform error
Theorem 5.3. Let u and u N p be the solutions of (2.1) and (5.1) respectively. Then
1 1 2 b A1 B1 N −1 2 + 2 + B N N −p u − u p < C ε ( p + 1)! 1 ( p !)2 N (( p + 1)!)2
(5.11)
for some constant C . 6. Numerical results and discussion In this section we take into account a number of test problems and present a rigorous comparative analysis. Test problems are solved using higher order finite element method over recursive modified graded mesh presented in this paper. The results so obtained are compared with the results obtained using higher order finite element method over Shishkin mesh, mesh with a rational and optimal mesh characterizing function, Bakhvalov-Shishkin mesh [19], Vulanovic´ mesh, GartlandShishkin mesh [25] and with modified Bakhvalov mesh [29]. In addition, results obtained are compared with B-spline method with artificial viscosity [13] and quintic B-spline method [20].
282
A. Kaushik et al. / Journal of Computational Physics 395 (2019) 275–285
Table 1 EN p , for Example 6.1 using Galerkin FEM with linear elements when
= 2−25 .
N
Bakhvalov-Shishkin mesh [19,25]
Gartland-Shishkin mesh [25]
Modified graded mesh
20 41 81 162 322 642
1.162(−2) 8.508(−3) 4.306(−3) 2.140(−3) 1.080(−3) 5.425(−4)
1.083(−2) 4.523(−3) 2.123(−3) 1.048(−3) 5.234(−4) 2.618(−4)
4.773(−4) 1.215(−4) 3.181(−5) 8.020(−6) 2.209(−6) 9.827(−7)
Table 2 EN p , and O rd( N ) for Example 6.1 using Galerkin FEM with quadratic elements when N
Shishkin mesh [25]
Vulanovic´ mesh [25]
Bakhvalov-Shishkin mesh [25]
Modified graded mesh
32
2.534(−3) 1.454 9.251(−4) 1.548 3.164(−4) 1.612 1.035(−4) 1.659 3.277(−5) 1.696 1.012(−5) 1.725 3.060(−6) –
3.818(−4) 1.902 1.022(−4) 1.935 2.674(−5) 1.946 6.939(−6) 1.953 1.792(−6) 1.960 4.605(−7) 1.967 1.178(−7) –
3.607(−4) 1.963 9.253(−5) 1.995 2.322(−5) 1.999 5.807(−6) 1.998 1.454(−6) 1.998 3.640(−7) 1.999 9.107(−8) –
1.249e−4 2.110 2.894e−5 2.435 5.351e−6 2.785 7.766e−7 2.936 1.015e−7 2.982 1.284e−8 2.995 1.610e−9 –
64 128 256 512 1024 2048
Suppose E N p , and O rd( N ) denotes the error in the convergence respectively. Moreover, we define
N EN p , = u − u p
= 2−12 .
and O rd( N ) =
-weighted energy norm and the corresponding numerical rate of
2N ln( E N p , ) − ln( E p , )
ln(2)
.
N Further, e N p , = u − u p N defines the maximum absolute error. In [19], results are presented for linear elements over Shishkin mesh, Bakhvalov-Shishkin mesh and for a graded mesh having a rational mesh generating function. Whereas, in [25] Galerkin FEM with quadratic and cubic elements are brought to attention. They consider five different S-type meshes; the original Shishkin mesh, a mesh with a rational mesh characterizing function, the Vulanovic´ mesh, the Bakhvalov-Shishkin mesh and the mesh with optimal mesh characterizing function ψopt . Moreover, a modified Bakhvalov mesh is the subject matter of [29]. In this section we present and compare results obtained over modified graded mesh with other adaptive/graded meshes developed over the last decades. Numerical results for Example 6.1 are conferred in Tables 1–3 and for Example 6.2 in Table 4. For linear elements, EN p , is tabulated and compared with the error obtained over Bakhvalov-Shishkin and Gartland-Shishkin meshes in the
form of Table 1. Moreover, for quadratic elements E N p , is tabulated and compared with the error obtained over Shishkin mesh, Vulanovic´ mesh, Bakhvalov-Shishkin mesh and Gartland-Shishkin mesh in the form of Table 2. For different values of perturbation parameter e N p , is tabulated and compared with an upwind scheme defined over a modified Bakhvalov mesh [29] and with a B-spline method with artificial viscosity [13]. The details of which can be had from Tables 3–4. Maximum absolute error (e N p , ) obtained using Galerkin FEM over the mesh presented is compared with a number of adaptive/graded meshes developed over the last decades (see Fig. 1). For different values of the method presented is compared with quintic B-spline method [20] (see Fig. 2). Example 6.1. Consider the following singular perturbation problem from [19,29,25]:
u + u − 2u = −e x−1 in = (0, 1), u (0) = u (1) = 0. Example 6.2. Consider the following singular perturbation problem from [7,20]:
u (x) + u (x) − u (x) = 0 in = (0, 1), u (0) = u (1) = 1.
A. Kaushik et al. / Journal of Computational Physics 395 (2019) 275–285
283
Table 3 Comparison of maximum absolute error (e N p , ) for Example 6.1 with an upwind scheme defined over a modified Bakhvalov mesh [29].
eN p ,
32
64
128
256
512
1024
2048
2−10
eN p ,
8.696(−4)
2.167(−4)
5.407(−5)
1.352(−5)
3.379(−6)
8.447(−7)
2.112(−7)
[29]
2.110(−2)
1.108(−2)
5.543(−3)
2.714(−3)
1.356(−3)
6.778(−4)
3.388(−4)
2−12
eN p ,
8.779(−4)
2.192(−4)
5.467(−5)
1.367(−5)
3.417(−6)
8.542(−7)
2.136(−7)
[29]
2.163(−2)
1.104(−2)
5.539(−3)
2.765(−3)
1.379(−3)
6.884(−4)
3.439(−4)
2−14
eN p ,
8.799(−4)
2.200(−4)
5.483(−5)
1.371(−5)
3.426(−6)
8.563(−7)
2.155(−7)
[29]
2.189(−2)
1.114(−2)
5.580(−3)
2.782(−3)
1.387(−3)
6.921(−4)
3.455(−4)
2−16
eN p ,
8.803(−4)
2.202(−4)
5.485(−5)
1.371(−5)
3.430(−6)
8.566(−7)
2.172(−7)
[29]
2.205(−2)
1.120(−2)
5.560(−3)
2.789(−3)
1.390(−3)
6.932(−4)
3.461(−4)
2−18
eN p ,
8.805(−4)
2.203(−4)
5.485(−5)
1.371(−5)
3.427(−6)
8.611(−7)
1.700(−7)
[29]
2.215(−2)
1.123(−2)
5.611(−3)
2.793(−3)
1.391(−3)
6.938(−4)
3.463(−4)
eN p ,
8.805(−4)
2.203(−4)
5.484(−5)
1.371(−5)
3.417(−6)
8.724(−7)
2.860(−7)
[29]
2.223(−2)
1.126(−2)
5.618(−3)
2.796(−3)
1.392(−3)
6.940(−4)
3.464(−4)
eN p ,
8.805(−4)
2.203(−4)
5.485(−5)
1.371(−5)
3.399(−6)
8.747(−7)
5.527(−7)
[29]
2.223(−2)
1.128(−2)
5.624(−3)
2.797(−3)
1.392(−3)
6.941(−4)
3.464(−4)
eN p ,
8.805(−4)
2.203(−4)
5.482(−5)
1.369(−5)
3.198(−6)
1.232(−6)
3.077(−7)
[29]
2.223(−2)
1.128(−2)
5.628(−3)
2.799(−3)
1.392(−3)
6.942(−4)
3.464(−4)
512
1024
2−20 2−22 2−24
Table 4 Comparison of maximum absolute error (e N p , ) for Example 6.2 with a B-spline method with artificial viscosity [13].
eN p ,
16
32
2−12
eN p ,
4.182(−3)
[13]
1.085(−2)
eN p ,
4.185(−3)
1.001(−3)
2.483(−4)
6.218(−5)
1.569(−5)
3.922(−6)
9.804(−7)
[13]
1.089(−2)
5.574(−3)
2.798(−3)
1.385(−3)
6.718(−4)
3.139(−4)
1.347(−4)
eN p ,
4.186(−3)
1.000(−3)
2.484(−4)
6.165(−5)
1.573(−5)
3.923(−6)
9.809(−7)
[13]
1.091(−2)
5.597(−3)
2.821(−3)
1.407(−3)
6.942(−4)
3.363(−4)
1.571(−4)
eN p ,
4.187(−3)
9.994(−4)
2.485(−4)
6.127(−5)
1.582(−5)
3.924(−6)
9.811(−7)
[13]
1.092(−2)
5.608(−3)
2.832(−3)
1.418(−3)
7.054(−4)
3.476(−4)
1.683(−4)
eN p ,
4.187(−3)
9.990(−4)
2.486(−4)
6.103(−5)
1.592(−5)
3.925(−6)
9.812(−7)
[13]
1.093(−2)
5.614(−3)
2.838(−3)
1.424(−3)
7.110(−4)
3.532(−4)
1.739(−4)
eN p ,
4.188(−3)
9.986(−4)
2.487(−4)
6.077(−5)
1.607(−5)
3.894(−6)
9.755(−7)
[13]
1.093(−2)
5.619(−3)
2.843(−3)
1.429(−3)
7.163(−4)
3.584(−4)
1.792(−4)
eN p ,
4.188(−3)
9.986(−4)
2.487(−4)
6.075(−5)
1.609(−5)
3.886(−6)
9.716(−7)
[13]
1.093(−2)
5.619(−3)
2.843(−3)
1.429(−3)
7.166(−4)
3.588(−4)
1.795(−4)
2−13 2−14 2−15 2−16 2−20 2−25
64
128
256
1.004(−3)
2.483(−4)
6.262(−5)
1.567(−5)
3.918(−6)
9.796(−7)
5.530(−3)
2.754(−3)
1.340(−3)
6.269(−4)
2.692(−4)
9.640(−5)
Fig. 1. Error (E N p , ) plot for Example 6.1, using Galerkin FEM with quadratic elements, over a mesh with a rational and optimal mesh characterizing function [25] and the modified graded mesh presented.
284
A. Kaushik et al. / Journal of Computational Physics 395 (2019) 275–285
Fig. 2. Comparison of maximum absolute error (e N p , ) for Example 6.2 with quintic B-spline method [20].
A comparison through numerical experiments reveals the dominance of the modified graded mesh over the state of the art meshes. Moreover, the results obtained are better than the results of B-spline method with artificial viscosity [13] and quintic B-spline method [20]. 7. Conclusion A singularly perturbed convection-diffusion problem is solved numerically using FEM based on higher order polynomials. A modified graded mesh is generated using some defined functions. Higher order convergence is obtained in implicitly -weighted energy norm. A convergence with O N − p is established. The error estimates obtained are optimal in the sense that they are free from logarithmic factor. Test examples are taken into account and a rigorous comparative analysis is presented. Moreover, we compare the proposed method and graded mesh generated with others found in the literature. References [1] N.S. Bakhvalov, On the optimization of methods for solving boundary value problems with boundary layers, Ž. Vyˇcisl. Mat. Mat. Fiz. 9 (1969) 841–859. ¸ Singularly Perturbed Boundary Value Problems, Birkhäusar, Basel, 2007. [2] L. Barbu, G. Morosanu, [3] M.G. Beckett, J.A. Mackenzie, On a uniformly accurate finite difference approximation of a singularly perturbed reaction diffusion problem using grid equidistribution, J. Comput. Appl. Math. 131 (2001) 381–405. [4] M. Brdar, H. Zarin, On graded meshes for a two parameter singularly perturbed problem, Appl. Math. Comput. 282 (2016) 97–107. [5] M. Brdar, H. Zarin, A singularly perturbed problem with two parameters on a Bakhvalov type mesh, J. Comput. Appl. Math. 292 (2016) 307–319. [6] M. Brdar, H. Zarin, L. Teofanov, A singularly perturbed problem with two parameters in two dimensions on graded meshes, Comput. Math. Appl. 72 (2016) 2582–2603. [7] C. Clavero, J.C. Jorge, F. Lisbona, G.I. Shishkin, An alternating direction scheme on a nonuniform mesh for reaction diffusion parabolic problem, IMA J. Numer. Anal. 20 (2000) 263–280. [8] R.G. Durán, A.L. Lombardi, Finite element approximation of convection diffusion problems using graded meshes, Appl. Numer. Math. 56 (2006) 1314–1325. [9] R.G. Durán, A.L. Lombardi, M.I. Prieto, Superconvergence for finite element approximation of a convection diffusion equation using graded meshes, IMA J. Numer. Anal. 32 (2012) 511–533. [10] Z. Gajic, M.T. Lim, Optimal Control of Singularly Perturbed Linear Systems and Applications, Marcel-Dekker, New York, 2001. [11] E.C. Gartland, Graded mesh difference schemes for singularly perturbed two point boundary value problems, Math. Comput. 51 (1988) 631–657. [12] G.W. Howell, Derivative error bounds for Lagrange interpolation: an extension of Cauchy’s bound for the error of Lagrange interpolation, J. Approx. Theory 67 (1991) 164–173. [13] M.K. Kadalbajoo, P. Arora, B-spline with artificial viscosity for solving singularly perturbed boundary value problems, Math. Comput. Model. 52 (2010) 654–666. [14] A. Kaushik, Singular perturbation analysis of bistable differential equation arising in the nerve pulse propagation, Nonlinear Anal. 9 (2008) 2106–2127. [15] A. Kaushik, M.D. Sharma, Numerical analysis of a mathematical model for propagation of an electrical pulse in a neuron, Numer. Methods Partial Differ. Equ. 27 (2008) 1–18. [16] N. Kopteva, N. Madden, M. Stynes, Grid equidistribution for reaction diffusion problems in one dimension, Numer. Algorithms 40 (2005) 305–322. [17] N. Kopteva, E. O’Riordan, Shishkin meshes in the numerical solution of singularly perturbed differential equations, Int. J. Numer. Anal. Model. 7 (2010) 393–415. [18] T. Linß, Analysis of a Galerkin finite element method on a Bakhvalov-Shishkin mesh for a linear convection diffusion problem, IMA J. Numer. Anal. 20 (2000) 621–632. [19] T. Linß, Layer Adapted Meshes for Reaction Convection Diffusion Problems, Springer, Berlin, 2010.
A. Kaushik et al. / Journal of Computational Physics 395 (2019) 275–285
285
[20] R.K. Lodhi, H.K. Mishra, Quintic B-spline method for solving second order linear and nonlinear singularly perturbed two point boundary value problems, J. Comput. Appl. Math. 319 (2017) 170–187. [21] J.J.H. Miller, Singular Perturbation Problems in Chemical Physics, John Wiley & Sons, New York, 1997. [22] J. Mohapatra, S. Natesan, Uniformly convergent numerical method for singularly perturbed differential difference equation using grid equidistribution, Int. J. Numer. Methods Biomed. Eng. 27 (2011) 1427–2445. [23] R.E. O’Malley, Introduction to Singular Perturbations, Academic Press, New York, 1974. [24] H.-G. Roos, Error estimates for linear finite elements on Bakhvalov type meshes, Appl. Math. 51 (2006) 63–72. [25] H.-G. Roos, Graded meshes for higher order FEM, J. Comput. Math. 33 (2015) 1–16. [26] H.-G. Roos, T. Linß, Sufficient conditions for uniform convergence on layer adapted grids, Computing 63 (1999) 27–45. [27] H.-G. Roos, T. Skalický, A comparison of the finite element method on Shishkin and Garland type meshes for convection diffusion problems, CWI Q. 10 (1997) 277–300. [28] H.-G. Roos, M. Stynes, L. Tobiska, Robust Numerical Methods for Singularly Perturbed Differential Equations, Springer, Berlin, 2008. [29] H.-G. Roos, L. Teofanov, Z. Uzelac, A modified Bakhvalov mesh, Appl. Math. Lett. 31 (2014) 7–11. [30] H. Schlichting, K. Gersten, Boundary Layer Theory, Springer-Verlag, Berlin, 2003. [31] K.K. Sharma, A. Kaushik, A solution of the discrepancy occurs due to using the fitted mesh approach rather than to the fitted operator for solving singularly perturbed differential equations, Appl. Math. Comput. 181 (2006) 756–766. [32] G. Söderlind, A.S. Yadaw, The impact of smooth W -grids in the numerical solution of singular perturbation two point boundary value problems, Appl. Math. Comput. 218 (2012) 6045–6055. [33] G. Sun, M. Stynes, Finite element methods for singularly perturbed higher order elliptic two point boundary value problems ii: convection diffusion type, IMA J. Numer. Anal. 15 (1995) 197–219. ´ A priori meshes for singularly perturbed quasilinear two point boundary value problems, IMA J. Numer. Anal. 21 (2001) 349–366. [34] R. Vulanovic, ´ L. Teofanov, A modification of the Shishkin discretization mesh for one dimensional reaction diffusion problems, Appl. Math. Comput. [35] R. Vulanovic, 220 (2013) 104–116. [36] C. Xenophontos, S. Franz, L. Ludwig, Finite element approximation of convection diffusion problems using an exponentially graded mesh, Comput. Math. Appl. 72 (2016) 1532–1540. [37] Q. Zheng, X. Zheng, Y.F. Liu, Uniform second order hybrid schemes on Bakhvalov-Shishkin mesh for quasilinear convection diffusion problems, Adv. Mater. Res. 871 (2014) 135–140.