Uniformly convergent high order finite element solutions of a singularly perturbed reaction–diffusion equation using mesh equidistribution

Applied Numerical Mathematics 39 (2001) 31–45 www.elsevier.com/locate/apnum

Uniformly convergent high order finite element solutions of a singularly perturbed reaction–diffusion equation using mesh equidistribution G. Beckett ∗ , J.A. Mackenzie Department of Mathematics, University of Strathclyde, Livingstone Tower, 26 Richmond Street, Glasgow, Scotland G1 1XH, UK

Abstract We study the numerical approximation of a singularly perturbed reaction–diffusion equation using a pth order Galerkin finite element method on a non-uniform grid. The grid is constructed by equidistributing a strictly positive monitor function which is a linear combination of a constant floor and a power of the second derivative of a representation of the boundary layers—obtained using a suitable decomposition of the analytical solution. By the appropriate selection of the monitor function parameters we prove that the numerical solution is insensitive to the size of the singular perturbation parameter and achieves the optimal rate of convergence with respect to the mesh density.  2001 IMACS. Published by Elsevier Science B.V. All rights reserved. Keywords: Uniform convergence; Adaptivity; Equidistribution; Singular perturbation; Reaction–diffusion; Finite element

1. Introduction This paper is concerned with the pth order Galerkin finite element discretization of the model singularly perturbed boundary value problem


Lε u := −ε 2 u (x) + b(x)2 u(x) = f (x), u(0) = 0, u(1) = 0.

x ∈ Ω := (0, 1),

(1)

Here, the constant of diffusion ε satisfies 0 < ε  β := infx∈Ω {b(x)}, and we assume that b(x) and f (x) are sufficiently smooth. Such problems are important in areas such as mathematical biology, the modeling of chemical reactions, and in solid mechanics for plate and shell modeling. It is widely appreciated that a standard numerical discretization of (1) is likely to be unreliable, due to the possible presence of steep * Corresponding author.

E-mail address: [email protected] (G. Beckett). 0168-9274/01/$ – see front matter  2001 IMACS. Published by Elsevier Science B.V. All rights reserved. PII: S 0 1 6 8 - 9 2 7 4 ( 0 1 ) 0 0 0 4 9 - 6

32

G. Beckett, J.A. Mackenzie / Applied Numerical Mathematics 39 (2001) 31–45

layers at the boundaries of the domain. These layers are of width O(ε) and so, as ε is reduced, any inaccuracies are likely to be accentuated. To overcome the above-mentioned difficulties, it is desirable to modify the discretization to ensure that the accuracy of the solution is insensitive to the size of ε—a property often referred to as robustness or uniform convergence. We define a numerical scheme to be kth order uniformly convergent in the norm  ·  on a mesh of N intervals if the error satisfies   u(x) − uN (x) < CN −k ,

where uN (x) is the numerical solution and C is a constant that is independent of N and ε. Strategies for achieving uniform convergence may be broken into two categories. The first approach involves obtaining sufficient a priori information about the form of the boundary layers and then tailoring the discretization accordingly. Many such schemes have been proposed in the literature, see, for example, Doolan et al. [12], O’Riordan and Stynes [19], and Roos [24]. These methods aim to produce accurate solutions for a low number of degrees of freedom. However, they are computationally expensive and difficult to extend to nonlinear and higher-dimensional problems. The second category involves the appropriate selection of a non-uniform grid—often referred to as mesh refinement—allowing a standard numerical discretization to be implemented. For the finite element method, there are three types of mesh refinement strategy; h-adaptivity involves adding points to (or removing points from) the domain, r-adaptivity involves moving a fixed number of points around the domain, and p-adaptivity keeps the mesh fixed and instead increasing the order of the local finite element basis functions. The necessary information required to implement these strategies should be deduced from the numerical solution, thus obviating the need for extensive a priori knowledge. A popular way to implement an r-adaptive scheme is to require that the L1 measure of a solutionbased monitor function is distributed equally between mesh cells. This is commonly referred to as mesh equidistribution. The concept is well establish in one dimension, see de Boor [11], Pereyra and Sewell [20], Russell and Christiansen [25], and White [27]. More recently, an extension of the concept to higher dimensions has been suggested by Huang and Russell [13]. Early error analysis in the literature was based on an asymptotic expansion of the truncation error with no consideration of uniform convergence, see Pereyra and Sewell [20] and Ascher et al. [1]. More recently, detailed error analysis by Qiu and Sloan [22], Qiu et al. [23], Mackenzie [16], and Linß [15] has shown that, for a monitor function based on the solution derivative, uniformly convergent numerical solutions of a convection-diffusion problem may be obtained using a low order finite difference scheme. In Beckett and Mackenzie [4, 5], we have constructed and analyzed a monitor function—called the BM monitor function—based on the second derivative of the solution, that can be applied to give uniformly convergent finite difference discretizations of both convection-diffusion and reaction–diffusion problems. The BM monitor function has also been successfully used within an adaptive moving mesh method to solve singularly perturbed evolutionary partial differential equations, see Beckett et al. [6]. The uniform convergence analysis in these references assumes that the grid is given a priori, though in practice the grid is generated in an adaptive manner based on the numerical solution: numerical experiments are used to justify the validity of this assumption. In this paper, we consider an r-adaptive refinement strategy using the BM monitor function with a pth order Galerkin finite element discretization of (1). We prove uniform convergence for the scheme and show how, by relating the monitor function to the polynomial order p, one may obtain the optimal rate of convergence in the L2 and energy norm.

G. Beckett, J.A. Mackenzie / Applied Numerical Mathematics 39 (2001) 31–45

33

Throughout this paper we will use standard finite element notation: H p (Ω) will represent the Sobolev p space of functions on Ω with L2 integrable derivative of order up to and including p. The space H0 (Ω) will represent the restriction of H p (Ω) to functions with compact support on Ω. In addition, C will denote a generic constant that is independent of ε and N , that can take different values at different places, even in the same argument.

2. Solution decomposition Lemma 1. For any prescribed order r, the solution u(x) of (1) has the form u(x) = v(x) + w(x), where the smooth component v(x) satisfies Lε v(x) = f (x),

x ∈ Ω,

v(0) = −w(0),

v(1) = −w(1),

and the singular component satisfies Lε w(x) = 0,

x ∈ Ω,

w(0) = −v(0),

w(1) = −v(1).

Furthermore, for 0
k
r the following bounds hold: and

 (k)  v (x) < C

(2)

 (k)    w (x) < Cε −k e−(b(0)/ε)x + e−(b(1)/ε)(1−x) .

(3)

Proof. See Lemma 2 of Beckett and Mackenzie [5].

✷

3. Discretization and non-uniform grids 3.1. Finite element discretization The weak formulation of the boundary value problem (1) is to find u ∈ H01 (Ω) such that for all z ∈ H01 (Ω),

B(u, z) = (f, z)Ω ,

(4)

where (· , ·)Ω is the L2 inner product on Ω and the bilinear form B : H01 (Ω) × H01 (Ω) → R is given by 

B(y, z) = ε 2 y  , z



Ω



+ b(x)2 y, z



Ω,

with associated energy norm zE,Ω = B(z, z)1/2. It is easy to show that B(u, z) is continuous and uniformly coercive on H01 (Ω) × H01 (Ω). The mapping z → (f, z)Ω is clearly a bounded linear functional and hence, from the Lax–Milgram lemma, we deduce that (4) has a unique solution.

34

G. Beckett, J.A. Mackenzie / Applied Numerical Mathematics 39 (2001) 31–45

We consider the Galerkin finite element discretization of (4) on a non-uniform partition ∆N := {0 = x0 < x1 < · · · < xN = 1}. Set hj = xj − xj −1 and denote by Ωj the mesh element (xj −1 , xj ). Let Πp (J ) be the space of polynomials of degree p on the interval J . Then, define the space of piecewise polynomials on ∆N by 

S p (∆N ) = z ∈ H 1 (Ω): zΩj ∈ Πp (Ωj ), j = 1, 2, . . . , N



and p

S0 (∆N ) = S p (∆N ) ∩ H01 (Ω). p

The Galerkin finite element solution uN ∈ S0 (∆N ) is such that p

B(uN , z) = (f, z),

for all z ∈ S0 (∆N ).

(5)

The standard error estimate for the Galerkin solution on an arbitrary mesh is given by the next theorem. Theorem 2. Let u be the solution of (1) and uN be the solution of (5) on an arbitrary mesh ∆N . Then u − uN E,Ω <

inf p

z∈S0 (∆N )

u − zE,Ω .

Proof. This is a standard result.

✷

It is easy to show that the above theorem will not lead to ε-uniform convergence on arbitrary grids. However, we will prove that for a suitably chosen mesh, it is possible to achieve ε-uniform convergence. 3.2. Grid equidistribution A commonly used technique in adaptive grid generation is mesh equidistribution: see for example de Boor [11] or Pereyra and Sewell [20]. A grid is said to be equidistributed if xj

1

1 M(s) ds = N

xj−1

M(s) ds,

j = 1, . . . , N,

(6)

0

L+ 1 (Ω)

:= {z ∈ L1 (Ω): essinf(z) > 0} is called a monitor function. Eq. (6) is often where M(x) ∈ referred to as the global equidistribution principle. It gives rise to an invertible mapping x(ξ ), from a computational domain (0, 1) onto the physical domain Ω, which satisfies x(ξ ) 

M(s) ds = ξ

0

1

M(s) ds.

(7)

0

If the equidistribution principle is used in an a posteriori adaptive algorithm, the monitor function will depend on the finite element solution uN . The coupled problem of finding the mesh and uN is then highly nonlinear, even though the original problem is linear. To make some headway with the analysis, we will consider an a priori layer adapted mesh based on the monitor function 

1/m

M = α + wN (x)

,

(8)

G. Beckett, J.A. Mackenzie / Applied Numerical Mathematics 39 (2001) 31–45

35

where wN is the singular component of uN , analogous to the decomposition in Lemma 2 and m is a positive integer. The parameter α is a positive constant that is independent of N , chosen to prevent mesh starvation and to improve the robustness of the adaptive grid procedure, see Pryce [21]. Often this floor is chosen in an ad hoc way. However, in [4,5] we show how to choose α so as to ensure uniform convergence for finite difference approximations of convection–diffusion and reaction–diffusion problems. 3.3. Grid structure We now consider the structure of a grid based on equidistribution of (8) and suggest how to choose α. In order to make headway with the analysis, we replace the approximation wN (x) in (8) by the leading order term of w  (x), as detailed in [5], to give        (x) = f (0)1/m ε −2/m e−(b(0)/mε)x + f (1)ε −2/m e−(b(1)/mε)(1−x). w (x)1/m ≈ W N

(9)

Hence, 1







(x) dx = mε (m−2)/m f (0)1/m b(0)−1 1 − e−b(0)/(mε) W

0



1/m

+ f (1)



b(1)−1 1 − e−b(1)/(mε)





=: Φ.

(10)

Using the assumption (9) and substituting (8) into (7) results in the mapping



    α α x(ξ ) + λ0 1 − e−b(0)/(mε)x(ξ ) + λ1 e−b(1)/(mε)(1−x(ξ )) − e−b(1)/(mε) = ξ +1 , Φ Φ

(11)

where λ0 = and

|f (0)|1/m b(0)−1 Ψ 

1/m

Ψ := f (0)

λ1 = 

|f (1)|1/m b(1)−1 Ψ 



1/m

b(0)−1 1 − e−b(0)/(mε) + f (1)

(12) 



b(1)−1 1 − e−b(1)/(mε) .

Notice immediately that 0
λ0 , λ1
1 + O(e−1/ε ). N A non-uniform grid {xj }N j =0 on the physical domain corresponds to an equispaced grid {ξj = j/N}j =0 on the computational domain. This identification gives the mesh



    j α α xj + λ0 1 − e−(b(0)/(mε))xj + λ1 e−(b(1)/(mε))(1−xj ) − e−b(1)/(mε) = +1 . Φ N Φ

(13)

The location of the grid point xj is given implicitly by the solution of the nonlinear algebraic equation (13). However, an important insight into the distribution of the mesh points is given by the following lemma. Lemma 3. Assuming that the non-uniform grid (13) is generated with α = Φ then ε ε log(N) < xkl +1 and xkr −1 < 1 − m log(N) < xkr , xkl < m b(0) b(1)

36

G. Beckett, J.A. Mackenzie / Applied Numerical Mathematics 39 (2001) 31–45

where





  ε 1 N log(N) + λ1 Ne−b(1)/(mε) N b(1)/b(0) − 1 λ0 (N − 1) + m kl = 2 b(0)

and





, 

  1 ε N log(N) + λ0 Ne−b(0)/(mε) N b(0)/b(1) − 1 λ1 (N − 1) + m kr = N − 2 b(1) and [·] denotes the integer part of. Furthermore, we have

e−(b(0)/(mε))xj < CN −1 ,

+ 1,

j  kl − 1

and e−(b(1)/(mε))(1−xj ) < CN −1 ,

j
kr + 1.

Proof. See Beckett and Mackenzie [5].

✷

Throughout the rest of this paper we will assume that α = Φ and that |wN |1/m in (8) may be approximated by (9) to give (x). M(x) =α+W

(14)

We will also assume that λ0 ε λ1 ε N log(N) < and m N log(N) < . (15) m b(0) 2 b(1) 2 Note that this is not an overly restrictive assumption if |v0 (0)| does not differ excessively from |v0 (1)| and N log(N) is small in relation to ε −1 . In fact, (15) is exactly the regime for which we are interested in using an adaptive algorithm to solve (1). Finally, we are only interested in analyzing what goes on when N −1  ε, otherwise a uniform grid could be used and a classical error analysis performed. Lemma 4. For j = 1, 2, . . . , N , the grid spacing satisfies hj < CN −1 . Proof. This is a direct consequence of the global equidistribution principle (6) and the choice of α = Φ. ✷ Lemma 5. The mesh spacing inside the boundary layers satisfies ε ε N −1 < hj < Cm , j = 1, 2, . . . , kl , Cm b(0) b(0) and ε ε N −1 < hj < Cm , j = kr + 1, . . . , N. Cm b(1) b(1) Proof. The upper bounds are proved in Beckett and Mackenzie [5]. For the lower bounds we only consider the lefthand boundary layer as the righthand layer may be treated in an analogous manner. When j = 1, 2, . . . , kl we have hj > h1 = x1 and so using (13), 







x1 + λ0 1 − e−b(0)x1/(mε) + λ1 e−(b(1)/(mε))(1−x1 ) − e−b(1)/(mε) = 2N −1 .

G. Beckett, J.A. Mackenzie / Applied Numerical Mathematics 39 (2001) 31–45

37

In addition, from (15) x1 < xkl <

λ0 mε log(N) < , b(0) 2N

and thus, x1 < N −1 . Finally, (e−(b(1)/(mε))(1−x1 ) − e−b(1)/(mε)) = O(e−1/ε ) and so 

  b(0) x1 > − log 1 − (λ0 N)−1 mε   mε mε −1 log 1 − (λ0 N)−1 > C N . ⇒ x1 > − b(0) b(0)

λ0 1 − e−(b(0)/(mε))x1 ) > N −1 ⇒

✷

4. Error analysis Lemma 6. Let k be an integer satisfying k
p + 1 and let z ∈ C(Ω) such that zΩj ∈ C k (Ωj ). Then there exists a piecewise polynomial interpolant of z, ip z ∈ S p (∆N ), and a constant C which is independent of hj and z, such that z(xj ) = ip z(xj ), j = 0, 1, . . . , N , and k−q

|z − ip z|q,∞,Ωj < Chj

|z|k,∞,Ωj

for q = 0, 1, . . . , k. Proof. See Theorem 3.1.5 of Ciarlet [9].

✷

Theorem 7. Let u be the solution of (1). Then on the equidistributed mesh (13) we have


|u − ip u|0,∞,Ωj < C and

N −(p+1) , N −m ,




ε|u − ip u|1,∞,Ωj < C

N −p , N −m ,

m  p + 1, m
(16)

m  p, m
(17)

for j = 1, 2, . . . , N . Proof. Using the decomposition of the exact solution we have ip u = ip v + ip w. We now bound the error in both terms separately. Firstly, for the smooth component using Lemmas 1, 4 and 6 with r  p + 1 and k = p + 1 we have p+1−q

|v − ip v|q,∞,Ωj < Chj

|v|p+1,∞,Ωj < CN −(p+1−q) ,

For the singular component we similarly have for q = 0, 1 that p+1−q

|w − ip w|q,∞,Ωj < Chj

|w|p+1,∞,Ωj .

q = 0, 1.

38

G. Beckett, J.A. Mackenzie / Applied Numerical Mathematics 39 (2001) 31–45

For j
kl − 1, we have from Lemmas 1 and 5 that p+1

|w − ip w|0,∞,Ωj < Chj < < < <

p+1 

|w|p+1,∞,Ωj = Chj



w (p+1)(η)

  p+1 Chj ε −(p+1) e−(b(0)/ε)η + e−(b(1)/ε)(1−η)  p+1  −(p+1) −(b(0)/ε)η Chj ε e +1  p+1  −(p+1) −(b(0)/ε)xj (b(0)/ε)hj Chj ε e e +1  p+1  −(p+1) −(b(0)/ε)xj Chj ε e +1 .

for some η ∈ Ωj

Therefore, using the global equidistribution principle (6), when m  p + 1 we have p+1  −(p+1) −(b(0)/ε)xj

|w − ip w|0,∞,Ωj < Chj

ε

e



 xj


(2−m)(p+1)/m



+1

p+1

ds M(s)



+N

−(p+1)

xj−1

< CN −(p+1) , and when m < p + 1 we have p+1  −(p+1) −(b(0)/ε)xj

|w − ip w|0,∞,Ωj < Chj

ε

e



(p+1)−m 2−(p+1) hj ε




+1

 xj

ds M(s)

m



+N

−(p+1)

xj−1

< CN

−m

.

For the error in interpolating the derivative, we can show using a similar argument that


ε|w − ip w|1,∞,Ωj <

CN −p , CN −m ,

m  p, m < p,

j
kl − 1.

Similar arguments establish the same results for j  kr + 2. Away from the boundary layer region, when kl
j
kr + 1, we have from Lemmas 3 and 6 with k = 0 |w − ip w|0,∞,Ωj < C|w|0,∞,Ωj 

< C e−(b(0)/ε)η + e−(b(1)/ε)(1−η)


e−(b(0)/ε)(xkl −1 /m)

m





+ e−(b(1)/ε)(xkr +1 /m)

< CN −m . For the error in the derivative we have for kl
j
kr + 1 ε 2 |w − ip w|1,∞,Ωj < Cε 2 |w|1,∞,Ωj < CN −m . This completes the proof.

✷

m 

G. Beckett, J.A. Mackenzie / Applied Numerical Mathematics 39 (2001) 31–45

39

In terms of the energy norm the next result follows immediately from Theorem 7. Corollary 8. Let u be the solution of problem (1). Let ip u be the piecewise polynomial interpolant from p S0 (∆N ) to u on the equidistributed mesh (13), as defined in Lemma 6. Then


u − ip uE,Ω < C

N −p , N −m ,

m  p, m < p.

(18)

We now state a uniform convergence result in the weighted energy norm. p

Theorem 9. Let uN ∈ S0 (∆N ) be the Galerkin finite element solution of (1) on the equidistributed mesh (13). Then


u − uN E,Ω < C

N −p , N −m ,

m  p, m < p.

(19) ✷

Proof. The result follows immediately from Theorem 2 and Corollary 8.

From the interpolation estimate given in Theorem 7, one may expect to obtain an improved rate of convergence in the L2 norm. This is usually obtained from the energy norm estimate using the Aubin– Nitche technique. Here, this will not produce sharp estimates due to the requirement that the error is decomposed into a smooth and a singular component. However, using a technique similar to that of Li and Navon [14], the higher rate of convergence can be proved for a slightly more stringent restriction on N using the following lemma. Lemma 10. The interpolating polynomial satisfies 



ip u − uN E < C 1 + εN + ε 1/2 N log(N)1/2 |u − ip u|0,∞ . Proof. Using the Galerkin orthogonality property ip u − uN 2E,Ω = B(ip u − uN , ip u − uN ) = B(ip u − u, ip u − uN ) 







= ε 2 (ip u − u) , χ  + b(x)2 (ip u − u), χ ,

(20)

where χ = ip u − uN . The first term on the right-hand side of (20) may be written as 



ε (ip u − u) , χ 2





 N         = ε (ip u − u) (εχ) dx    i=1 Ω

i

  N       = ε (ip u − u)(εχ) dx    i=1 Ω

i


ε|ip u − u|0,∞


 S1



εχ dx +

 S2



εχ dx +







εχ dx , R

where S1 = (0, xkl ), S2 = (xkr +1 , 1), and R = (xkl , xkr +1 ). We consider each of these integrals in turn. For  the boundary layer components, we use the standard inverse estimate χ  L2 (Ωj ) < Ch−1 j χ L2 (Ωj ) , see

40

G. Beckett, J.A. Mackenzie / Applied Numerical Mathematics 39 (2001) 31–45

Brenner and Scott [7]. Therefore, using Lemma 5 we have 





εχ  dx
meas(S1 )1/2 εχ  L2 (S1 )

S1



= meas(S1 )

1/2

kl  j =1



< Cmeas(S1 )

1/2

1/2

εχ  2L2 (Ωj )

kl  j =1

< <





1/2

εχ  2 h−2 j L2 (Ωj ) 



εχ   Cε log(N)1/2 h−1 1 L (S )   2 1 −1/2 1/2   Cε N log(N) εχ L2 (S1 ) . 1/2

Similarly, for the right-hand boundary layer we have 





εχ  dx < Cε −1/2 N log(N)1/2 εχ  L2 (S2 ) .

S2

Away from the boundary layers, we have 





εχ  dx < CN εχ  L2 (R) .

R

Combining these estimates together gives 









ε 2 (ip u − u) , χ  < C|ip u − u|0,∞,Ω εN + ε 1/2 N log(N)1/2 εχ  L2 (Ω) .

(21)

Furthermore, for the second component of (20), we may easily obtain 



b(x)2 (ip u − u), χ
C|ip u − uN |0,∞,Ω b(x)χL2 (Ω) .

Combining (21) and (22) gives 

(22)



ip u − uN E,Ω < C 1 + εN + ε 1/2 N log(N)1/2 |u − ip u|0,∞,Ω , ✷

and the proof is complete. p

Theorem 11. Let uN ∈ S0 (∆N ) be the Galerkin finite element solution of (1) on the equidistributed mesh (13). Then, for ε 1/2 N log(N)1/2 < C, we have


u − uN L2 (Ω) < C

N −(p+1) , N −m ,

m  p + 1, m < p + 1.

Proof. Using Theorem 7 and Lemma 10, we have u − uN L2 (Ω)
u − ip uL2 (Ω) + ip u − uN L2 (Ω)
u − ip uL2 (Ω) + ip u − uN E,Ω
CN −(p+1) , which completes the proof.

✷

(23)

G. Beckett, J.A. Mackenzie / Applied Numerical Mathematics 39 (2001) 31–45

41

The restriction on N imposed by Theorem 11 is stricter than a normal adaptivity criterion, though the numerical results presented in the next section suggest that the improved rate of convergence holds for more general N . Remark 12. Theorems 9 and 11 prove that if m  p + 1 then the Galerkin approximation converges at the optimal rate with respect to N .

5. Adaptive mesh generation The generation of the Galerkin finite element solution of (1) requires the completion of two tasks; firstly the adaptive mesh given by (13) has to be determined, and then the finite element solution has to be computed on this mesh. To provide a benchmark for truly adaptive computations, we will confirm the error analysis by solving an example problem using the exactly equidistributing grid for the monitor function (14)—additional examples, including a problem with a non-constant reaction coefficient, are presented in Beckett [3]. The problem we will consider here is










−ε 2 u (x) + u(x) = 2 1 + ε 2 π 2 cos(π x) + 1 + 16π 2 ε 2 cos(4π x), u(0) = 0, u(1) = 0. This problem has the exact solution u(x) = A− e−x/ε + A+ e−(1−x)/ε + 2 cos(π x) + cos(4π x),

(24)

(25)

where 3 + e−1/ε 1 + 3e−1/ε + and A = . 1 − e−2/ε 1 − e−2/ε The solution has two boundary layers of equal thickness O(ε) and of height 3 and −1, respectively. In Tables 1 and 2, we present numerical results in the L∞ , L2 , and energy norms using the exactly equidistributing grid, piecewise cubic basis functions, and with single perturbation parameter ε = 10−4 and ε = 10−8 , respectively. In both cases, we see that the error converges at the expected rate in each norm: fourth order for the L∞ and L2 norms, and third order for the energy norm. For the case of ε = 10−8 , we see that the energy norm actually converges at the higher than expected rate of O(N −4 ). As will be explained in more detail below, this is due to the dominance of the L2 component of the error. Comparing analogous results in Tables 1 and 2, we see ε independent convergence in all three norms. In fact, for the energy norm we see better than uniform convergence. For small values of N , the energy norm is dominated by the L2 component and, as asserted by the error analysis of the previous section, this converges at fourth order with N . In contrast, the H1 component of the error only converges at third order. For ε = 10−4 , the H1 component of the error becomes appreciable as N increases and the convergence rate drops to third order. However, for the case when ε = 10−8 —because of the scaling of the energy norm—the H1 component is two orders of magnitude smaller and does not become significant over the range of N considered here. In practice, generation of the adaptive grid requires that the global equidistribution principle (6) and monitor function (8) be discretized. The resulting set of equations could then be coupled with the finite element discretization (5) to give a nonlinear system for ∆N and uN . However, this system is both large A− = −

42

G. Beckett, J.A. Mackenzie / Applied Numerical Mathematics 39 (2001) 31–45

and restrictive: it is likely to be expensive to solve and dictates that the grid and numerical solution be evaluated with the same level of accuracy. Thus, a popular alternative is to decouple the calculation of the grid from the finite element method, and solve the two sets of equations in an iterative manner. Firstly the finite element solution (5) is computed on a fixed grid: this is used to compute the monitor function (8) which can be exactly equidistributed to give a new grid. This procedure is repeated until some measure of convergence is achieved. Details of the discretization of the monitor function and the solution of the decoupled system of equations is given in [3]: what follows is a brief summary. Table 1 Problem (24) using exactly equidistributing grid with ε = 10−4 , p = 3 and m = 4 N

eL∞

Rate

eL2

Rate

eE

Rate

16

2.475 × 10−3

–

8.811 × 10−4

–

8.886 × 10−4

–

32

2.234 × 10−4

3.47

5.398 × 10−5

4.03

5.590 × 10−5

3.99

64

1.414 × 10−5

3.98

3.301 × 10−6

4.03

3.756 × 10−6

3.90

128

8.834 × 10−7

4.00

2.070 × 10−7

4.00

3.048 × 10−7

3.62

256

5.563 × 10−8

3.99

1.292 × 10−8

4.00

3.081 × 10−8

3.31

Table 2 Problem (24) using exact equidistribution with ε = 10−8 , p = 3 and m = 4 N

eL∞

Rate

eL2

Rate

eE

Rate

16

2.458 × 10−3

–

8.808 × 10−4

–

8.808 × 10−4

–

32

2.492 × 10−4

3.30

5.424 × 10−5

4.02

5.424 × 10−5

4.02

64

1.415 × 10−5

4.14

3.305 × 10−6

4.03

3.305 × 10−6

4.03

128

8.865 × 10−7

4.00

2.069 × 10−7

4.00

2.069 × 10−7

4.00

256

5.674 × 10−8

3.97

1.297 × 10−8

4.00

1.298 × 10−8

4.00

Table 3 Problem (24) using recovered equidistributing grid with ε = 10−4 , p = 3, m = 4, NINIT = 32 and mINIT = 2.0 N

eL∞

Rate

eL2

Rate

eE

Rate

16

2.475 × 10−3

–

8.811 × 10−4

–

8.886 × 10−4

–

32

2.234 × 10−4

3.47

5.398 × 10−5

4.03

5.590 × 10−5

3.99

64

1.414 × 10−5

3.98

3.301 × 10−6

4.03

3.756 × 10−6

3.90

128

8.834 × 10−7

4.00

2.070 × 10−7

4.00

3.048 × 10−7

3.62

256

5.563 × 10−8

3.99

1.292 × 10−8

4.00

3.081 × 10−8

3.31

G. Beckett, J.A. Mackenzie / Applied Numerical Mathematics 39 (2001) 31–45

43

Table 4 Problem (24) using recovered equidistributing grid with ε = 10−8 , p = 3, m = 4, NINIT = 32, and mINIT = 2.0 N

eL∞

Rate

eL2

Rate

eE

Rate

16

2.458 × 10−3

–

8.808 × 10−4

–

8.808 × 10−4

–

32

2.492 × 10−4

3.30

5.424 × 10−5

4.02

5.424 × 10−5

4.02

64

1.415 × 10−5

4.14

3.305 × 10−6

4.03

3.305 × 10−6

4.03

128

8.865 × 10−7

4.00

2.069 × 10−7

4.00

2.069 × 10−7

4.00

256

5.674 × 10−8

3.97

1.297 × 10−8

4.00

1.298 × 10−8

4.00

The complete solution {∆N , uN } is generated in two stages. In the first stage the adaptive grid is computed. It is desirable that this be done in a robust and efficient manner. Because of this, it is inappropriate to use a high degree piecewise polynomial numerical solution, since in the initial stages of adaptivity the grid is likely to be some distance from equidistribution with little or no points inside the boundary layers. This leads to oscillations in the solution in the neighbourhood of the boundary layers and hence the approximation to the monitor function will be unreliable. Instead, a piecewise linear finite element solution is generated on a coarse mesh containing NINIT cells, using the decoupled iterative procedure outlined above. This has the additional advantage of reducing the number of unknowns which have to be computed at each iteration. Then, the standard Galerkin orthogonality condition is used in conjunction with our knowledge of the exponential form of the boundary layers to recover a more accurate rendering of the singular component of the solution, as described in Barrat and Morton [2]. From this a more accurate approximation to the grid is generated, on which the final finite element discretization is computed. In Tables 3 and 4 we present results for the recovery procedure outlined above. As can be seen, the results presented are in excellent agreement with those of the exact grid for both ε = 10−4 and ε = 10−8 . 5.1. Guidance on the selection of m The analysis given in this paper provides some guidance on a suitable value for the parameter m: specifically, to achieve the optimal rate of convergence with N , it is necessary to choose m  p + 1. We now investigate further, through numerical experiments, how the solution accuracy is affected by m. In Fig. 1 we plot the error, measured in the energy norm, for the exactly equidistributing grid (13) for a range of values of m. The error decreases rapidly as m increases from 2. It is clear that the error plots for different values of N diverge over the range 2
m
4, confirming that the convergence rate with respect to N increases with m for m
p + 1. For all values of N the minimum error occurs in the range 72
m
4. This agrees with the analysis of Carey and Dinh [8] which shows that, for a piecewise cubic interpolant of (25), the minimum error in the energy norm will be achieved for the grid with m = 72 . As m increases beyond 4, the error slowly begins to rise. This indicates that, for reliability and efficiency, it is better to choose a value of m slightly larger than the optimal rather than slightly smaller.

44

G. Beckett, J.A. Mackenzie / Applied Numerical Mathematics 39 (2001) 31–45

Fig. 1. Effect of varying m on the numerical solution for problem (24) using mesh equidistribution with ε = 10−4 , p = 3 and N = 16, 32, 64, and 128.

6. Conclusions In this paper we have shown that uniform convergence of the pth order Galerkin finite element method may be achieved using an adaptive non-uniform grid obtained by equidistribution. We have proved that the uniform convergence rates are optimal with respect to the number of mesh elements. The error analysis has been carried out in the energy and L2 norms. Numerical experiments and optimal interpolation estimates suggest that these results are also true in the L∞ norm. Future work will concentrate on the extension of these results to two-dimensional problems. Depending upon the nature of the domain, the adaptive mesh will also need to detect and resolve corner layers.

Acknowledgements The authors would like to thank Torsten Linß (Technische Universität Dresden) for a number of useful discussions and suggestions.

References [1] U.M. Ascher, M.M. Mattheij, R.D. Russell, Numerical solution of boundary value problems for ordinary differential equations, Prentice-Hall, Englewood Cliffs, NJ, 1988. [2] J.W. Barrett, K.W. Morton, Optimal finite element solutions to diffusion-convection problems in one dimension, Internat J. Numer. Methods Eng. 15 (1980) 1457–1474.

G. Beckett, J.A. Mackenzie / Applied Numerical Mathematics 39 (2001) 31–45

45

[3] G. Beckett, The robust and efficient numerical solution of singularly perturbed boundary value problems using grid adaptivity, Ph.D. Thesis, Department of Mathematics, University of Strathclyde, 1998. [4] G. Beckett, J.A. Mackenzie, Convergence analysis of finite difference approximations on equidistributed grids to singularly perturbed boundary value problems, Appl. Numer. Math. 35 (2000) 87–109. [5] G. Beckett, J.A. Mackenzie, On a uniformly accurate finite difference approximation of a singularly perturbed reaction–diffusion problem using grid equidistribution, Comput. Appl. Math. 131 (2001) 381–405. [6] G. Beckett, J.A. Mackenzie, A. Ramage, D.M. Sloan, On the numerical solution of one-dimensional PDEs using adaptive methods based on equidistribution, J. Comput. Phys. 167 (2001) 372–392. [7] S.C. Brenner, L.R. Scott, The Mathematical Theory of Finite Element Methods, Springer, Berlin, 1994. [8] G.F. Carey, H.T. Dinh, Grading functions and mesh redistribution, SIAM J. Numer. Anal. 22 (1985) 1028– 1040. [9] P.G. Ciarlet, The Finite Element Method for Elliptic Problems, North-Holland, Amsterdam, 1978. [10] J.M. Coyle, J.E. Flaherty, R. Ludwig, On the stability of mesh equidistribution strategies for time-dependent partial differential equations, J. Comput. Phys. 62 (1986) 26–39. [11] C. De Boor, Good approximation by splines with variables knots II, in: Lecture Notes in Mathematics, Vol. 363, 1974, pp. 12–20. [12] E.P. Doolan, J.J.H. Miller, W.H.A. Schilders, Uniform Numerical Methods for Problems with Initial and Boundary Layers, Boole Press, Dublin, 1980. [13] W. Huang, R.D. Russell, A high dimensional moving mesh strategy, Appl. Numer. Math. 26 (1997) 63–76. [14] J. Li, I.M. Navon, Global uniformly convergent finite element methods for singularly perturbed elliptic boundary value problems: higher-order elements, Comput. Methods Appl. Mech. Eng. 171 (1999) 1–23. [15] T. Linß, Uniform pointwise convergence of finite difference schemes using grid equidistribution, Computing 66 (2001) 27–39. [16] J.A. Mackenzie, Uniform convergence analysis of an upwind finite-difference approximation of a convectiondiffusion boundary value problem on an adaptive grid, IMA J. Numer. Anal. 19 (1999) 233–249. [17] J.M. Melenk, On the robust exponential convergence of hp finite element methods for problems with boundary layers, IMA J. Numer. Anal. 17 (1997) 577–601. [18] J.H.H. Miller, E. O’Riordan, G.I. Shishkin, Fitted Numerical Methods for Singular Perturbation Problems, World Scientific, Singapore, 1995. [19] E. O’Riordan, M. Stynes, A uniformly accurate finite element method for a singularly perturbed onedimensional reaction–diffusion problem, Math. Comp. 47 (1986) 555–570. [20] V. Pereyra, E.G. Sewell, Mesh selection for discrete solution of boundary problems in ordinary differential equations, Numer. Math. 23 (1975) 261–268. [21] J.D. Pryce, On the convergence of iterated remeshing, IMA. J. Numer. Anal. 9 (1989) 315–335. [22] Y. Qiu, D.M. Sloan, Analysis of difference approximations to a singularly perturbed two-point boundary value problem on an adaptively generated grid, J. Comput. Appl. Math. 101 (1999) 1–25. [23] Y. Qiu, D.M. Sloan, T. Tang, Convergence analysis for an adaptive finite difference solution of a variable coefficient singular perturbation problem, J. Comput. Appl. Math. 116 (2000) 121–143. [24] H.-G. Roos, Global uniformly convergent schemes for a singularly perturbed boundary value problem using patch base spline-functions, J. Comput. Appl. Math. 29 (1990) 69–77. [25] R.D. Russell, J. Christiansen, Adaptive mesh selection strategies for solving boundary value problems, SIAM J. Numer. Anal. 15 (1978) 59–80. [26] C. Schwab, M. Suri, The p and hp versions of the finite element method for problems with boundary layers, Math. Comp. 65 (1996) 1403–1429. [27] A.B. White Jr., On selection of equidistributing meshes for two-point boundary problems, SIAM J. Numer. Anal. 16 (1979) 472–502.