discontinuous Galerkin methods for singularly perturbed problems of convection-diffusion type

Applied Numerical Mathematics 76 (2014) 48–59

Contents lists available at ScienceDirect

Applied Numerical Mathematics www.elsevier.com/locate/apnum

Higher order uniformly convergent continuous/discontinuous Galerkin methods for singularly perturbed problems of convection-diffusion type ✩ Peng Zhu ∗ , Shenglan Xie School of Mathematics, Physics and Information, Jiaxing University, Jiaxing, Zhejiang 314001, China

a r t i c l e

i n f o

a b s t r a c t

Article history: Received 5 July 2012 Received in revised form 16 June 2013 Accepted 7 October 2013 Available online 15 October 2013

In this paper, we propose and analyze a higher order continuous/discontinuous Galerkin methods for solving singularly perturbed convection-diffusion problems. Based on piecewise polynomial approximations of degree k  1, a uniform convergence rate O ( N −k lnk N ) in associated norm is established on Shishkin mesh, where N is the number of elements. Numerical experiments complement the theoretical results. © 2013 IMACS. Published by Elsevier B.V. All rights reserved.

Keywords: Convection-diffusion equation Local discontinuous Galerkin method Finite element method Shishkin mesh Uniform convergence

1. Introduction In this paper we are interested in the construction and validation of high-order finite element approximations to problems of type



−ε u  + bu  + cu = f

in Ω = (0, 1),

u (0) = u (1) = 0,

(1.1)

where 0 < ε  1 is a small positive parameter, and b, c , f are sufficiently smooth functions with the following properties

b(x)  b0 > 0,

c (x)  0,

c (x) −

1  ¯ b (x)  c 0 > 0, ∀x ∈ Ω, 2

(1.2)

for some constants b0 and c 0 . This assumption guarantees that (1.1) has a unique solution in H 2 (Ω) ∩ H 01 (Ω) for all f ∈ L 2 (Ω) [16]. Typically the solution of (1.1) has an exponential boundary layer at x = 1. Problem (1.1) is a simple model problem that helps understanding the behavior of numerical methods in presence of layers in more complex problems like the Navier–Stokes equations in fluid dynamics or convection-diffusion equations in chemical reaction processes. The smallness of ε causes global unphysical oscillations if standard discretization schemes on general meshes are applied. Therefore, layer-adapted meshes have to be used to obtain efficient discretizations. Early ideas of layer-adapted meshes go ✩ Supported by Zhejiang Provincial Natural Science Foundation of China (Grant No. LQ12A01014) and Science Research Foundation of Jiaxing University of China (Grant No. 70510017). Corresponding author. E-mail addresses: [email protected] (P. Zhu), [email protected] (S. Xie).

*

0168-9274/$36.00 © 2013 IMACS. Published by Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.apnum.2013.10.001

P. Zhu, S. Xie / Applied Numerical Mathematics 76 (2014) 48–59

49

back to Bakhvalov [2]. The piecewise uniform Shishkin meshes [10] were proposed originally for finite difference schemes. Due to the instability of the standard Galerkin method even on layer-adapted meshes, see [6] for numerical results, stabilized discretizations have to be considered. For more information on this topic we refer to the recent book by Roos et al. [9]. When ε approaches zero, the problem (1.1) changes from an elliptic equation to a hyperbolic one. Inspired by the great success of the discontinuous Galerkin (DG) method in solving hyperbolic equations, Cockburn and Shu [5] firstly presented a local discontinuous Galerkin (LDG) method for solving time-dependent convection-diffusion problem. Then Celiker and Cockburn [4] analyzed the superconvergent properties of the LDG method for one-dimensional steady convection-diffusion equations. Recently Xie et al. [12,13] and Zhang et al. [17] investigated the uniform superconvergent properties of the LDG method for one-dimensional singularly perturbed convection-diffusion equation. On the other hand, nonsymmetric discontinuous Galerkin method with interior penalty (the NIPG method), originally designed for elliptic equations, is analyzed by Zarin and Roos [15] for two-dimensional convection-diffusion problems with parabolic layers. A disadvantage of DG method is that the method produces more degrees of freedom than the continuous finite element method (CFEM). The combination of DG and CFEM is an efficient way to overcome this disadvantage. With this motivation, Roos and Zarin [8], Zarin [14] analyzed the uniformly convergence of NIPG/CFEM coupled method with bilinear element on Shishkin mesh for two dimensional convection-diffusion problems with exponentially layers or characteristic layers. In their work, the domain was split into the coarse and the fine part, then the CFEM with bilinear element was adopted in the fine part where the mesh size is comparable with  , and the NIPG method with bilinear element is used in the coarse part for its stabilization. A coupled LDG/CFEM approach has also been studied by Perugia and Schötzau [7] for the modeling of elliptic problems arising in electromagnetics. Zhu et al. [18,19] extended the LDG/CFEM coupled method with linear element to one-dimensional convection-diffusion problems and analyzed its convergence on quasi-uniform mesh [18] and its uniformly convergence on Shishkin mesh [19] respectively. The works in [8,14,19] suggest that the CFEM and DG coupled approach is efficient and robust for solving singularly perturbed problem. But only bilinear element case in the NIPG/CFEM coupled method and linear element case in the LDG/CFEM coupled method were analyzed on Shishkin mesh. In the present work, we will investigate the uniform convergence properties of the LDG/CFEM coupled methods with higher order elements for one-dimensional singularly perturbed problems with convection-diffusion type on Shishkin meshes. In the sequel, with C we shall denote a generic positive constant independent of the perturbation parameter  and mesh size. 2. The layer-adapted mesh and the LDG/CFEM coupled method 2.1. The layer-adapted mesh Let N be an even integer. Denote by fine. This parameter is given by



τ = min

1

,

σ

2 b0

τ the transition parameter which indicate where the mesh changes from coarse to



ε ln N .

The parameter σ is typically chosen to equal the formal order of the method or to accommodate the error analysis. Notice that ε  1, here and below we take τ = bσ ε ln N. Moreover, we suppose that ε  N −1 which is realistic for this type of 0 problems. The domain Ω is divided into two subdomains Ω1 = [0, 1 − τ ] and Ω2 = [1 − τ , 1]. Let T N = { I j = (x j −1 , x j ): j = 1, . . . , N } be a partition of the domain Ω . We choose



xj =

2(1 − τ ) j / N ,

j = 0, 1 , . . . , N /2 ,

1 − 2τ ( N − j )/ N ,

j = N /2 + 1 , . . . , N ,

for Shishkin meshes. N /2 Set T N1 = { I j } j =1 and T N2 = { I j } Nj= N /2+1 . We denote by u (x+ ) and u (x−j ) the values of u at x j , from the right cell and the j left cell of x j , respectively. 2.2. The weak formulation of the LDG/CFEM coupled method We discretize problem (1.1) by using the LDG method in Ω1 , and CFEM in Ω2 . Let u i = u |Ωi , i = 1, 2, and q = (u 1 ) in Ω1 . Rewrite problem (1.1) as the following equivalent transmission problem:

⎧   q − u1 = 0 ⎪ ⎪ ⎪ ⎪  1  ⎪  1 ⎪ ⎪ ⎨ −εq + b u + cu = f  2   2  −ε u + b u + cu 2 = f ⎪ ⎪ ⎪ ⎪ ⎪ u 1 (1 − τ ) = u 2 (1 − τ ), ⎪ ⎪   ⎩ q(1 − τ ) = u 2 (1 − τ ),

in Ω1 , in Ω1 , in Ω2 ,

(2.1)

50

P. Zhu, S. Xie / Applied Numerical Mathematics 76 (2014) 48–59

with boundary conditions

u 1 (0) = u 2 (1) = 0.

(2.2)

Let us now denote by Pk ( I j ) the space of polynomials of degree at most k on I j , and define the finite element space 1,k

VN

2,k

and V N

as follows,







V N = v 1 ∈ L 2 (Ω1 ): v 1  I ∈ Pk ( I j ), ∀ I j ∈ T N1 , 1,k

j







V N = v 2 ∈ H 1 (Ω2 ): v 2 (1) = 0, v 2  I ∈ Pk ( I j ), ∀ I j ∈ T N2 . 2,k

j

2,k

1,k

The space V N is a standard conforming finite element space, whereas the functions in V N are completely discontinuous across interelement boundaries. 1,k 1,k 2,k We will search for approximate solutions (q N , u 1N , u 2N ) of (2.1) and (2.2) in the finite element space V N × V N × V N that satisfy (2.1) and (2.2) in a weak sense. The LDG/CFEM coupled method (see more details in [18,19]) for problem (2.1) and (2.2) is defined as: 1,k 1,k 2,k Find (q N , u 1N , u 2N ) ∈ V N × V N × V N such that





Ij

Ij











Ij



c − b u 1N v 1 dx −

εq N − b u 1N v 1 dx + Ij

1

+ εqˆ − bu˜ (x j −1 ) v x+j −1 = 





+ uˆ 1 (x j −1 ) w x+j −1 = 0, u 1N w  dx − uˆ 1 (x j ) w x− j

q N w dx +

 1







εqˆ − bu˜ 1 (x j ) v 1 x−j

(2.3)



f v 1 dx,

(2.4)

Ij 1,k

1,k

for all test function ( w , v 1 ) ∈ V N × V N

 







ε u 2N − bu 2N v 2 dx +

Ω2

and for all I j ∈ T N1 , and





c − b u 2N v 2 dx +





Ω2 2,k VN ,

εqˆ − bu˜ 1 (1 − τ ) v 2 (1 − τ ) =

f v 2 dx,

(2.5)

Ω2

ˆ1

˜1

for all test function v ∈ where u , u and qˆ are the numerical fluxes, which approximate the traces of u 1N and q N on the boundary of the elements of T N1 . To complete the specification of the method, it only remains to define the numerical fluxes. 2

2.2.1. The numerical fluxes We use the following notation to describe the numerical fluxes at the interior nodes. The average and jump of the trace of smooth function v at the interior node x j are given by





v (x j ) =

v (x+ ) + v (x−j ) j 2

and













− v x−j , v (x j ) = v x+ j

respectively. We now define the numerical fluxes uˆ 1 and qˆ by

⎧ if j = 0, ⎪ ⎨0 1 1 1 uˆ (x j ) = {u N (x j )} − β[u N (x j )] if j = 1, . . . , N /2 − 1, ⎪ ⎩ 2 u N (x j ) if j = N /2,

and

⎧ q (x+ ) + α u 1N (x+ ) if j = 0, ⎪ ⎪ j ⎨ N j qˆ (x j ) = {q N (x j )} + β[q N (x j )] + α [u 1N (x j )] if j = 1, . . . , N /2 − 1, ⎪ ⎪ ⎩ q (x− ) + α (u 2 (x ) − u 1 (x− )) if j = N /2. N j N N j j

(2.6)

(2.7)

Here the scalar α = α (x) and β = β(x) are auxiliary parameters. Their purpose is to enhance the stability and accuracy properties of the LDG method (see [3] and [5]). The numerical flux associated with the convection is the classical upwinding flux, namely,



u˜ (x j ) = 1

0 u 1N (x− ) j

if j = 0, if j = 1, . . . , N /2.

(2.8)

P. Zhu, S. Xie / Applied Numerical Mathematics 76 (2014) 48–59

51

3. Stability and error analysis of the coupled method This section is devoted to the existence and uniqueness of the solution of the coupled method (2.3)–(2.5) with numerical fluxes (2.6)–(2.8), and its corresponding error analysis. Firstly, we rewrite our method in the primal form by eliminating q following Arnold et al. [1]. And then we get stability of the CFEM–LDG coupled method, if the stabilization parameter α is taken of order O (1/ H ). Under this condition, we obtain the higher order uniform convergence of the coupled method. 3.1. Primal formulation We denote by H 1 (T N1 ) the space of functions on Ω1 whose restriction to each element K belongs to the Sobolev space H ( K ). A straightforward computation shows that for all φ ∈ H 1 (T N1 ) and ψ ∈ H 1 (T N1 ), 1

φψ  dx = −

Ω1

φ  ψ dx −

  



  φ(x j ) ψ(x j ) + φ(x j ) ψ(x j ) j =1

Ω1

    + φ (1 − τ )− ψ (1 − τ )− − φ 0+ ψ 0+ .

Summing up (2.3) for all I j ∈

N /2−1





q N − u 1N



T N1 ,

and combining with (2.6) and (3.1), we obtain

N /2−1

w dx =

(3.1)

 



u 1N (x j )







w (x j ) + β w (x j )

  + u 1N 0+ w 0+

j =1

Ω1

   + u 2N (1 − τ ) − u 1N (1 − τ )− w (1 − τ )− . 1,k

(3.2)

2,k

k Let us introduce the space V N := { v = ( v 1 , v 2 ): v 1 ∈ V N , v 2 ∈ V N } and V k ( N ) = [ H 2 (Ω) ∩ H 01 (Ω)] + V Nk , where H 01 (Ω) =

{ v ∈ H (Ω): v |Ω1 = v , v |Ω2 = 1

1

v , v (0+ ) 2

1

satisfying

N /2−1

L1 ( v )r dx =

 

= 0, v (1) = 0}. For v ∈ V k ( N ), we define L1 ( v ) as the unique element in V N1,k



v 1 (x j )

2







r (x j ) + β r (x j )

  + v 1 0+ r 0+

j =1

Ω1

   + v 2 (1 − τ ) − v 1 (1 − τ )− r (1 − τ )−

(3.3)

1,k VN .

for all r ∈ As a result, (3.2) can be rewritten as



q N = u 1N



+ L1 (u N ) in V N1,k .

(3.4)

Summing up (2.4) for all I j ∈ T N1 and adding together with (2.5), we get







 

εq N − b u 1N v 1 dx +

Ω1







N /2−1

ε u 2N − b u 2N v 2 dx +

 





    c − b u N v dx + εqˆ − b u˜ (1 − τ ) v 2 (1 − τ ) − v 1 (1 − τ )− + εqˆ − b u˜ 1 (0) v 1 0+ =

+







j =1

Ω2



εqˆ − bu˜ 1 (x j ) v 1 (x j )



1



Ω

f v dx. Ω

Inserting (2.7) and (2.8) into the equation above, and recalling the definition of L1 (·), we obtain



εq N v

1 

+ L1 ( v ) dx +

Ω1

ε



   u 2N v 2 dx −

Ω2 N /2−1

+









bu N v dx + Ω









c − b u N v dx

Ω



εα u 1N (x j ) v 1 (x j ) + εα u 2N (1 − τ ) − u 1N (1 − τ )−





v 2 (1 − τ ) − v 1 (1 − τ )−



j =1 N /2−1      1     b(x j )u 1N x− v (x j ) − b(1 − τ )u 1N (1 − τ )− v 2 (1 − τ ) − v 1 (1 − τ )− + εα u 1N 0+ v 1 0+ − j j =1

=

f v dx. Ω

(3.5)

52

P. Zhu, S. Xie / Applied Numerical Mathematics 76 (2014) 48–59 1,k

Similar to the definition of L1 ( v ), for v ∈ V k ( N ), we define L2 ( v ) as the unique element in V N

N /2−1



bL2 ( v )r dx =











satisfying



+ b(1 − τ ) v 2 (1 − τ ) − v 1 (1 − τ )− r (1 − τ )− b(x j ) v 1 (x j ) r x− j



(3.6)

j =1

Ω1 1,k

for all r ∈ V N . As a consequence, (3.5) can be rewritten as









εq N v 1 + L1 ( v ) dx + Ω1

ε u 2N

 



v 2 dx −

Ω2

+





Ω

Ω

c − b u N v dx −

N /2−1

bL2 ( v )u 1N dx +

+ εα

u 2N (1 −

τ

) − u 1N









εα u 1N (x j ) v 1 (x j )

j =1

Ω1



bu N v  dx



−

(1 − τ )



   v (1 − τ ) − v (1 − τ )− + εα u 1N 0+ v 1 0+ = 2

1

f v dx.

(3.7)

Ω k Inserting (3.4) into (3.7), we obtain the so-called primal form of our coupled method which reads: find u N ∈ V N such that

A N (u N , v ) := B N (u N , v ) + C N (u N , v ) + S N (u N , v ) = F N ( v ) ∀ v ∈ V Nk with





Ω



bw v  + L2 ( v ) dx +

CN ( w , v ) = − Ω

FN ( v ) =

f v dx, Ω





c − b w v dx,

Ω

N /2−1



SN ( w , v ) =



ε w  + L1 ( w ) v  + L1 ( v ) dx,

BN ( w , v ) =

(3.8)













εα w 1 (x j ) v 1 (x j ) + εα w 1 0+ v 1 0+



j =1

    + εα w 2 (1 − τ ) − w 1 (1 − τ )− v 2 (1 − τ ) − v 1 (1 − τ )− . Here, L1 (·) and L2 (·) have been defined in L 2 (Ω) by a trivial extension. From the following lemma, the primal formulation is consistent. Lemma 3.1. Let u be the exact solution of the problem (2.1) and (2.2). Then the primal form (3.8) has the Galerkin orthogonality property k for all v ∈ V N .

A N ( u − u N , v ) = 0,

(3.9)

Proof. Since u is the exact solution, we have [u 1 (x j )] = 0, j = 1, . . . , N /2 − 1, u 1 (0) = 0, u 1 (1 − τ ) = u 2 (1 − τ ), (u 1 ) (1 − τ ) = (u 2 ) (1 − τ ). Consequently, S N (u , v ) = 0, L1 (u ) = 0. Then for all v ∈ V Nk , we have

B N (u , v ) =





ε u  v  + L1 ( v ) dx

Ω



ε u1

=

 



v 1 dx +

Ω1





εL1 ( v ) u 1 dx +

Ω1



ε u2

 



v 2 dx.

(3.10)

Ω2

Taking φ = ε (u 1 ) and ψ = v 1 in (3.1), we have



ε u1 Ω1

 



v 1 dx = − Ω1





ε u 1 v 1 dx − 

+ε u

 1 

N /2−1



ε







u 1 (x j )



v 1 (x j ) +







u 1 (x j )



v 1 (x j )

j =1

     (1 − τ )− v 1 (1 − τ )− − ε u 1 0+ v 1 0+ .

(3.11)

P. Zhu, S. Xie / Applied Numerical Mathematics 76 (2014) 48–59

53

Inserting (3.11) into (3.10), integrating by parts the third term in the right-hand side of (3.10), and combining with the definition of L1 (·), we obtain

−ε u  v dx +

B N (u , v ) =

N /2−1







 











ε u 1 (x j ) β v 1 (x j ) − v 1 (x j ) + ε u 2 (1) v 2 (1).

j =1

Ω

k By the definition of V N , we get v 2 (1) = 0. This, together with [u  ] = 0, yields

−ε u  v dx,

B N (u , v ) =

k for all v ∈ V N .

Ω

Now we consider the term C N (u , v ). It can be rewritten as

bu 1

C N (u , v ) = −



v1



+ L2 ( v ) dx −

Ω1





bu 2 v 2 dx +

Ω2





c − b uv dx.

(3.12)

Ω

Taking φ = −bu 1 and ψ = v 1 in (3.1), we have





bu 1 v 1 dx =

− Ω1

 



b u1

N /2−1   



  + b u 1 v 1 dx + b (x j ) u 1 (x j ) v 1 (x j ) + u 1 (x j ) v 1 (x j )

Ω1 1



−

− b(1 − τ )u (1 − τ )



j =1

v

1



  (1 − τ )− + b(0)u 1 0+ v 1 0+ .

(3.13)

Inserting (3.13) into (3.12), integrating by parts the second term in the right-hand side of (3.12), and combining with the definition of L2 (·), we obtain



C N (u , v ) =

N /2−1





bu + cu v dx +









b(x j ) v 1 x+ u 1 (x j ) . j

j =1

Ω

Recalling [u 1 (x j )] = 0, j = 1, . . . , N /2 − 1, and −ε u  + bu  + cu = f , we have

A N (u , v ) =

∀ v ∈ V Nk .

f v dx, Ω

2

In view of (3.8), (3.9) obviously holds. 3.2. Stability analysis

To consider the stability of the primal form A N , we define the following norms and seminorms for v ∈ V k ( N ):

| v |2ε =  v 20,Ω + ε | v |21, N + ε | v |2∗ + | v |c2 , N /2   1  2     2  v  , | v |21, N =  v 2 0,Ω + 0, I 2

N /2−1 2

| v |∗ = | v |c2 =





j

j =1



2



α v 1 ( x j ) + α v 1 0+

2

  2 + α v 2 (1 − τ ) − v 1 (1 − τ )− ,

j =1 N /2−1 1 

2



2

b (x j ) v 1 (x j )

j =1

 2   2 1  1 + b(0) v 1 0+ + b(1 − τ ) v 2 (1 − τ ) − v 1 (1 − τ )− , 2

2

where  · 0, D is the usual Sobolev norm defined on region D. Lemma 3.2. If α = O (1/ H ), there exists a constant C 1 > 0, such that

A N (u N , u N )  C 1 |u N |2ε ,

∀u N ∈ V Nk .

(3.14)

54

P. Zhu, S. Xie / Applied Numerical Mathematics 76 (2014) 48–59

k Proof. By direct computation, we obtain, for all u N ∈ V N ,

B N (u N , u N ) =



2

ε u N + L1 (u N ) dx, SN (u N , u N ) = ε|u N |2∗ ,

Ω

 C N (u N , u N ) =



c−

1  2 b u N dx + |u N |c2 . 2

Ω

Adding together the above three equality, and combining with (1.2), we have

A N (u N , u N ) 



2 ε u N + L1 (u N ) dx + c0 

Ω

= ε |u N |21, N + 2ε

u 2N dx + ε |u N |2∗ + |u N |c2 Ω



2

u N L1 (u N ) dx + ε L1 (u N )0,Ω + c 0 u N 20,Ω + ε |u N |2∗ + |u N |c2 .

Ω

Applying the arithmetic-geometric mean inequality, we have, for every θ > 0,

 2   A N (u N , u N )  ε (1 − θ)|u N |21, N + (1 − 1/θ)L1 (u N )0,Ω + c 0 u N 20,Ω + ε |u N |2∗ + |u N |c2 . According to [7] (p. 422), when the coefficient

  L1 (u )

0,Ω

 C |u |∗ ,

α = O( H −1 ), there exists a constant C > 0, such that

∀u ∈ V k ( N ).

Inserting (3.16) into (3.15), for any θ satisfying



(3.15)

(3.16) C2 C 2 +1

< θ < 1, we easily get,

 A N (u N , u N )  ε (1 − θ)|u N |21, N + ε C 2 (1 − 1/θ) + 1 |u N |2∗ + c 0 u N 20,Ω + |u N |c2   γ ε |u N |21, N + ε |u N |2∗ + c 0 u N 20,Ω + |u N |c2  min{γ , c 0 }|u N |2ε , where

γ = min{1 − θ, C 2 (1 − 1/θ) + 1}. Taking C 1 = min{γ , c0 }, the proof is complete. 2

From Lemma 3.2, we easily get

|u N |ε  C  f 0,Ω , which implies the uniqueness of the solution to (3.8). Further, since (3.8) is a linear problem over the finite-dimensional k space V N , the existence of the solution follows from its uniqueness. Consequently, by (3.4), we get the existence and uniqueness of the solution to the problem (2.3)–(2.5) with numerical fluxes (2.6)–(2.8). Remark 1. In fact, following [7] or [12], for any α  0, we can prove the existence and uniqueness of the solution to the problem (2.3)–(2.5) with numerical fluxes (2.6)–(2.8). In this paper, we are only interested in the special case α = O ( H −1 ). 3.3. Error analysis We are now going to provide a ε -uniform estimate for the error u − u N in the norm (3.14). The error analysis presented in this paper relies on a priori estimate of the exact solution of (1.1) and a special interpolation which was firstly introduced in [11]. Lemma 3.3. (See [9].) Let q be some positive integer. Consider the boundary value problem (1.1) with the assumption of (1.2). Its exact solution u can be composed as u = S + E, where the smooth part S and the layer part E satisfies

−ε S  + b S  + c S = f , −ε E  + bE  + c E = 0, and

 (l)   S (x)  C ,

   (l)   E (x)  C  −l exp − b0 (1 − x) for 0  l  q.



(3.17)

P. Zhu, S. Xie / Applied Numerical Mathematics 76 (2014) 48–59

55

Next we introduce a special interpolant in [11] that will be useful later. On each element K = [x j −1 , x j ], we define k + 1 nodal functionals Nl by

N 0 ( v ) = v (x j −1 ), Nl ( v ) =

N k ( v ) = v (x j ),

x j

1

(x j − x j −1 )l

(x − x j −1 )l−1 v (x) dx,

l = 1, . . . , k − 1.

x j −1

Now a local interpolation v I | K ∈ Pk ( K ) is defined by Nl ( v I − v ) = 0, l = 0, . . . , k, which can be extended to a continuous k global interpolation v I ∈ V N . Obviously, if k = 1, this special interpolation just is the Lagrange linear interpolation. Lemma 3.4. Let the exact solution u = S + E of the problem (1.1) can be decomposed into a smooth and layer part, respectively, S I and E I their interpolants on a Shishkin mesh (assuming ε ln N  b0 /(2σ )). Then, we have u I = S I + E I and the estimates

 k+1 u − u I L ∞ (Ω2 )  C N −1 ln N ,  −(k+1) − u − u I L ∞ (Ω1 )  C N +N σ ,   ( S − S I )(l)  2  C N l−(k+1) , l = 0, . . . , k, L (Ω)  k+1  E − E I L 2 (Ω2 )  C ε 1/2 N −1 ln N ,    1/2 −σ −1    N E I L 2 (Ω ) +  E I  L 2 (Ω1 )  C ε N + N −(1/2+σ ) ,

(3.18a) (3.18b) (3.18c) (3.18d) (3.18e)

1

 E L ∞ (Ω1 ) + ε −1/2  E L 2 (Ω1 )  C N −σ ,   E  2  C ε −1/2 N −σ . L (Ω )

(3.18f) (3.18g)

1

Proof. (3.18a)–(3.18f) can be obtained directly from Lemma 12 of [11]. We only have to check (3.18g). By Lemma 3.3, we obtain

  2 E  2

L (Ω1

1

−τ

  2  E (x) dx 

= ) 0

=



1 2b0 ε

Thus,  E   L 2 (Ω1 )  √ 1

2b0

1

−τ



ε−2 exp −

b0 (1 − x)

 dx

ε

0

    2b0 τ 2b0 − exp −  exp −

ε

ε



1 2b0 ε

exp −

2b0 τ



ε

=

1 2b0 ε

N −2 σ .

ε−1/2 N −σ . 2

The following statement is the direct consequence of Lemma 3.4. Theorem 3.1. Under the conditions of Lemma 3.4 and taking σ = k + 1, we have η = u − u I satisfies

 k |η|ε  C N −1 ln N

(3.19)

for Shishkin mesh. Proof. Since u − u I is continuous in Ω , we have |η|∗ = 0, |η|c = 0. Then, |η|2ε = η20,Ω + ε |η|21, N . By (3.18c), (3.18d), (3.18e) and (3.18f) of Lemma 3.4, we easily conclude

u − u I 0,Ω   S − S I 0,Ω +  E − E I 0,Ω2 +  E 0,Ω1 +  E I 0,Ω1

    C N −(k+1) 1 + ε 1/2 (ln N )k+1 + N −σ ε 1/2 + N −1/2 .

(3.20)

On Ω2 , by (3.17), we have

| E − E I |1,Ω2 =

1/2

 

|E −

E I |21, K

K ⊂Ω2

 C

2τ N

k   K ⊂Ω2 K



 

 C

K ⊂Ω2



ε−2(k+1) exp −

2τ

2k

 (k+1) 2 E 

0, K

N

2β(1 − x)

ε



1/2 dx

1/2

56

P. Zhu, S. Xie / Applied Numerical Mathematics 76 (2014) 48–59

 k  C ε N −1 ln N · ε −(k+1) ·



 exp −

2β(1 − x)



ε

1/2 dx

Ω2

 k  C ε −1/2 N −1 ln N .

(3.21)

On Ω1 , due to (3.18e) and (3.18g) we get

     | E − E I |1,Ω1   E   L 2 (Ω ) +  E I  L 2 (Ω )  C ε −1/2 + ε 1/2 N + N 1/2 N −σ . 1

1

(3.22)

Combined with (3.18c), (3.21) and (3.22), we obtain



ε1/2 |u − u I |1,Ω  ε1/2 | S − S I |1,Ω + | E − E I |1,Ω2 + | E − E I |1,Ω1



    C ε 1/2 + (ln N )k N −k + 1 + ε N + ε 1/2 N 1/2 N −σ .

Taking

(3.23)

σ = k + 1, by (3.20) and (3.23), we easily get (3.19). 2

Now we turn to estimate |ξ |ε . Theorem 3.2. Under the conditions of Lemma 3.4 and taking σ = k + 1. Assuming α = O (1/ H ), then ξ = u I − u N satisfies

 k |ξ |ε  C N −1 ln N .

(3.24)

Proof. By Lemma 3.1 and Lemma 3.2, we first obtain

C 1 |ξ |2ε  A N (ξ, ξ ) = −A N (η, ξ ) = −B N (η, ξ ) − C N (η, ξ ) − S N (η, ξ ). By the definition of u I , we have

(3.25)

η(x j ) = 0, j = 0, 1, . . . , N. Consequently,





εη ξ  + L1 (ξ ) dx,

B N (η, ξ ) = Ω

C N (η, ξ ) = −

bηξ  dx +

Ω



c − b



ηξ dx ≡ I 1 + I 2 ,

Ω

S N (η, ξ ) = 0.

(3.26)

Then, by (3.16) and (3.19), we have

        k B N (η, ξ )  C ε η  2 |ξ |1, N + |ξ |∗  C ε 1/2 η  L 2 (Ω) |ξ |ε  C N −1 ln N |ξ |ε . L (Ω)

(3.27)

The first term in the right-hand side of (3.26) can be estimated with,

     | I 1 |  C ηL ∞ (Ω1 ) ξ   L 1 (Ω ) + ηL ∞ (Ω2 ) ξ   L 1 (Ω ) . 1

2

On Ω1 we use an inverse inequality to estimate

  ξ 

L 1 (Ω1 )

 C N ξ L 1 (Ω1 )  C N |ξ |ε ,

while on Ω2 we use Cauchy–Schwarz inequality to estimate

  ξ 

L 1 (Ω2 )



√   τ ξ 

0,Ω2

 C (ln N )1/2 |ξ |ε .

These two bounds and (3.18a), (3.18b) of Lemma 3.4 yield

   k | I 1 |  C N −k + N −(k+1) (ln N )k+3/2 |ξ |ε  C N −1 ln N |ξ |ε , where we have used the fact (ln N )3/2  C N. Taking σ = k + 1 in (3.20), we have

η0,Ω  C N −(k+1/2) . Hence, the second term in the right-hand side of (3.26) can be easily estimated with

| I 2 |  C η0,Ω ξ 0,Ω  C N −(k+1/2) |ξ |ε .

(3.28)

P. Zhu, S. Xie / Applied Numerical Mathematics 76 (2014) 48–59

57

This, combined with (3.28), yields

   C N (η, ξ )  C N −1 ln N k |ξ |ε .

(3.29)

2

Collecting (3.25), (3.27) and (3.29), we have (3.24).

The combination of Theorem 3.1 and Theorem 3.2 leads to our main results directly, i.e.,

Theorem 3.3. Let u and u N be the solutions of the continuous problem (2.1) and the discrete problem (3.8), respectively. Taking

σ = k + 1 and assuming α = O(1/ H ), then  k |u − u N |ε  C N −1 ln N .

(3.30)

Corollary 3.1. Let (q N , u N ) be the solution obtained by the coupled method (2.3)–(2.5) with numerical fluxes (2.6)–(2.8). Under the assumption of Theorem 3.3, we have

   (q − q N , u − u N )  C N −1 ln N k , A

(3.31)

N

where |(·, ·)|AN is a problem-related norm defined by

    (r , v )2 =  v 2 + ε r 2 + ε  v 2  2 + ε | v |2 + | v |2 . c ∗ 0,Ω 0,Ω1 A 0,Ω N

2

(3.32)

Proof. From (3.4), we have q − q N = (u 1 − u 1N ) − L1 (u N ). Since L1 (u ) = 0, we obtain q − q N = (u 1 − u 1N ) + L1 (u − u N ). In terms of (3.16), we conclude |(q − q N , u − u N )|AN  C |u − u N |ε , which implies the conclusion. 2

4. Numerical experiments

In this section, we numerically verify the sharpness of our theoretical findings. In our numerical experiments, we take

α = 1/ H , β = 1/2 in (2.6) and (2.7), σ = k + 1 in the transition parameter of Shishkin mesh. Example. We solve the model problem (1.1) with b = c = 1, and









f (x) = (cos x) 1 + e (x−1)/ε + (1 + ε )(sin x) 1 − e (x−1)/ε . The exact solution is u (x) = (sin x) (1 − e (x−1)/ε ), which exhibits a boundary layer with the width O (ε ln 1ε ) at the outflow boundary x = 1. We display the history of convergence of our coupled method in Table 1. A Shishkin mesh with N elements is called mesh N. Let err( N ) denote the error of the approximate solution computed on the mesh N. Then the approximate order of convergence, i.e., order(2N ), is defined by

order (2N ) :=

ln(err( N )/err (2N )) ln(2 ln( N )/ ln(2N ))

.

From Table 1 we observe that the numerical results for the problem-related norm (3.32) of the error agree with those predicted in Corollary 3.1.

58

P. Zhu, S. Xie / Applied Numerical Mathematics 76 (2014) 48–59

Table 1 History of convergence of the coupled method, under the norm (3.32), Shishkin mesh, k

N

ε = 1.0e−03

σ = k + 1.

ε = 1.0e−04

ε = 1.0e−05

error

order

error

order

error

order

1

16 32 64 128 256 512

1.175746e−01 7.408547e−02 4.460330e−02 2.605187e−02 1.489342e−02 8.378784e−03

– 0.98 0.99 1.00 1.00 1.00

1.174759e−01 7.402205e−02 4.456482e−02 2.602932e−02 1.488050e−02 8.371513e−03

– 0.98 0.99 1.00 1.00 1.00

1.174661e−01 7.401573e−02 4.456098e−02 2.602707e−02 1.487921e−02 8.370788e−03

– 0.98 0.99 1.00 1.00 1.00

2

16 32 64 128 256 512

2.250080e−02 9.138153e−03 3.344589e−03 1.145197e−03 3.747559e−04 1.186607e−04

– 1.92 1.97 1.99 2.00 2.00

2.247004e−02 9.125263e−03 3.339808e−03 1.143552e−03 3.742164e−04 1.184898e−04

– 1.92 1.97 1.99 2.00 2.00

2.246696e−02 9.123973e−03 3.339329e−03 1.143387e−03 3.741624e−04 1.184727e−04

– 1.92 1.97 1.99 2.00 2.00

3

16 32 64 128 256 512

4.421321e−03 1.161643e−03 2.590585e−04 5.205066e−05 9.753553e−06 1.738388e−06

– 2.84 2.94 2.98 2.99 3.00

4.413009e−03 1.159369e−03 2.585423e−04 5.194622e−05 9.733933e−06 1.734886e−06

– 2.84 2.94 2.98 2.99 3.00

4.412177e−03 1.159142e−03 2.584906e−04 5.193575e−05 9.731966e−06 1.734533e−06

– 2.84 2.94 2.98 2.99 3.00

1

16 32 64 128 256 512

1.174651e−01 7.401510e−02 4.456060e−02 2.602684e−02 1.487908e−02 8.370715e−03

– 0.98 0.99 1.00 1.00 1.00

1.174650e−01 7.401503e−02 4.456056e−02 2.602682e−02 1.487907e−02 8.370708e−03

– 0.98 0.99 1.00 1.00 1.00

1.174650e−01 7.401503e−02 4.456055e−02 2.602682e−02 1.487907e−02 8.370707e−03

– 0.98 0.99 1.00 1.00 1.00

2

16 32 64 128 256 512

2.246666e−02 9.123844e−03 3.339281e−03 1.143370e−03 3.741570e−04 1.184710e−04

– 1.92 1.97 1.99 2.00 2.00

2.246662e−02 9.123831e−03 3.339276e−03 1.143369e−03 3.741565e−04 1.184708e−04

– 1.92 1.97 1.99 2.00 2.00

2.246662e−02 9.123830e−03 3.339276e−03 1.143369e−03 3.741565e−04 1.184708e−04

– 1.92 1.97 1.99 2.00 2.00

3

16 32 64 128 256 512

4.412093e−03 1.159119e−03 2.584854e−04 5.193471e−05 9.731752e−06 1.734492e−06

– 2.84 2.94 2.98 2.99 3.00

4.412085e−03 1.159116e−03 2.584849e−04 5.193439e−05 9.731696e−06 1.734704e−06

– 2.84 2.94 2.98 2.99 3.00

4.412085e−03 1.159118e−03 2.584877e−04 5.193465e−05 9.733665e−06 1.735449e−06

– 2.84 2.94 2.98 2.99 3.00

ε = 1.0e−06

ε = 1.0e−07

ε = 1.0e−08

5. Conclusions In this paper, we introduce a higher order coupled LDG–CFEM method for solving singularly perturbed convectiondiffusion problems, and analyze its stability and uniform convergence property on the Shishkin mesh. For kth order element, a rate O (( N −1 ln N )k ) in an associated norm is established. Our numerical experiments indicate the sharpness of this error estimate. References [1] D.N. Arnold, F. Brezzi, B. Cockburn, L.D. Marini, Unified analysis of discontinuous Galerkin methods for elliptic problems, SIAM J. Numer. Anal. 39 (2002) 1749–1779. [2] N.S. Bakhvalov, The optimization of methods of solving boundary value problems with a boundary layer, USSR Comput. Math. Math. Phys. 9 (1969) 139–166. [3] P. Castillo, B. Cockburn, I. Perugia, D. Schötzau, An a priori error analysis of the local discontinuous Galerkin method for elliptic problems, SIAM J. Numer. Anal. 38 (2000) 1676–1706. [4] F. Celiker, B. Cockburn, Superconvergence of the numerical traces of discontinuous Galerkin and Hybridized methods for convection-diffusion problems in one space dimension, Math. Comput. 76 (2007) 67–96. [5] B. Cockburn, C.-W. Shu, The local discontinuous Galerkin finite element method for convection-diffusion systems, SIAM J. Numer. Anal. 35 (1998) 2440–2463. [6] T. Linß, M. Stynes, Asymptotic analysis and Shishkin-type decomposition for an elliptic convection-diffusion problem, J. Math. Anal. Appl. 261 (2001) 604–632. [7] I. Perugia, D. Schötzau, On the coupling of local discontinuous Galerkin and conforming finite element methods, J. Sci. Comput. 16 (2001) 411–433. [8] H.-G. Roos, H. Zarin, A supercloseness result for the discontinuous Galerkin stabilization of convection-diffusion problems on Shishkin meshes, Numer. Methods Partial Diff. Equ. 23 (6) (2007) 1560–1576.

P. Zhu, S. Xie / Applied Numerical Mathematics 76 (2014) 48–59

59

[9] H.-G. Roos, M. Stynes, L. Tobiska, Robust Numerical Methods for Singularly Perturbed Differential Equations, Springer Ser. Comput. Math., vol. 24, Springer-Verlag, Berlin, Heidelberg, 2008. [10] G.I. Shishkin, Grid approximation of singularly perturbed elliptic and parabolic equations, Second doctoral thesis, Keldysh Institute, Moscow, 1990 (in Russian). [11] L. Tobiska, Analysis of a new stabilized higher order finite element method for advection-diffusion equations, Comput. Methods Appl. Mech. Eng. 196 (2006) 538–550. [12] Z.Q. Xie, Z. Zhang, Superconvergence of DG method for one-dimensional singularly perturbed problems, J. Comput. Math. 25 (2007) 185–200. [13] Z.Q. Xie, Z.Z. Zhang, Z. Zhang, A numerical study of uniform superconvergence for solving singularly perturbed problems, J. Comput. Math. 27 (2009) 280–298. [14] H. Zarin, Continuous-discontinuous finite element method for convection-diffusion problems with characteristic layers, J. Comput. Appl. Math. 231 (2009) 626–636. [15] H. Zarin, H.-G. Roos, Interior penalty discontinuous approximations of convection-diffusion problems with parabolic layers, Numer. Math. 100 (2005) 735–759. [16] Z. Zhang, Finite element superconvergence on Shishkin mesh for 2-D convection-diffusion problems, Math. Comput. 245 (2003) 1147–1177. [17] Z.Z. Zhang, Z.Q. Xie, Z. Zhang, Superconvergence of discontinuous Galerkin methods for convection-diffusion problems, J. Sci. Comput. 41 (2009) 70–93. [18] P. Zhu, Z.Q. Xie, S.Z. Zhou, A coupled continuous-discontinuous FEM approach for convection-diffusion equations, Acta Math. Sci. 31B (2011) 601–612. [19] P. Zhu, Z.Q. Xie, S.Z. Zhou, A uniformly convergent continuous-discontinuous Galerkin method for singularly perturbed problems of convection-diffusion type, Appl. Math. Comput. 217 (2011) 4781–4790.