Adobe Acrobat Reader | Vanity Fair HDTV 720p AC3 5.1 | adobe acrobat

A C1-conforming hp finite element method for fourth order singularly perturbed boundary value problems

A C1-conforming hp finite element method for fourth order singularly perturbed boundary value problems

JID:APNUM AID:2992 /FLA [m3G; v1.172; Prn:15/02/2016; 13:24] P.1 (1-17) Applied Numerical Mathematics ••• (••••) •••–••• Contents lists available a...

935KB Sizes 0 Downloads 1 Views

JID:APNUM AID:2992 /FLA

[m3G; v1.172; Prn:15/02/2016; 13:24] P.1 (1-17)

Applied Numerical Mathematics ••• (••••) •••–•••

Contents lists available at ScienceDirect

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

A C 1 -conforming hp finite element method for fourth order singularly perturbed boundary value problems Pandelitsa Panaseti a , Antri Zouvani a , Niall Madden b , Christos Xenophontos a,∗ a b

Department of Mathematics and Statistics, University of Cyprus, P.O. Box 20537, 1678 Nicosia, Cyprus Department of Mathematics, Statistics and Applied Mathematics, National University of Ireland, Galway, Ireland

a r t i c l e

i n f o

Article history: Available online xxxx Keywords: hp finite element method Fourth order singularly perturbed boundary value problem Boundary layers Exponential convergence Spectral Boundary Layer Mesh

a b s t r a c t We consider a fourth order singularly perturbed boundary value problem (BVP) in onedimension and the approximation of its solution by the hp version of the Finite Element Method (FEM). The given problem’s boundary conditions are not suitable for writing the BVP as a second order system, hence the approximation will be sought from a finite dimensional subspace of the Sobolev space H 2 . We construct suitable C 1 hierarchical basis functions for the approximation and we show that the hp FEM on the Spectral Boundary Layer Mesh yields a robust approximation that converges exponentially in the energy norm, as the number of degrees of freedom is increased. Numerical examples that validate (and extend) the theory are also presented. © 2016 IMACS. Published by Elsevier B.V. All rights reserved.

1. Introduction: the model problem Singularly Perturbed Problems (SPPs) arise in numerous applications from science and engineering, such as electrical networks, vibration problems, and in the theory of hydrodynamic stability [3,12]. In such problems, the highest derivative in the differential equation is multiplied by a very small positive parameter. This causes the solution to contain boundary layers, which are rapidly varying solution components that have support in a narrow neighborhood of the boundary of the domain. Their numerical approximation must be carefully constructed in order for their effects to be accurately captured [14]. In the context of the FEM, this requires the use of graded meshes which include refinement near the boundary layer region that depends on the singular perturbation parameter. Examples of such meshes include the Bakhvalov [1] and Shishkin [16] meshes, which are used with finite differences and the h version of the FEM, as well as the Spectral Boundary Element Mesh [8] which is used with the p and hp versions of the FEM. Most SPPs that have been studied in the literature concern second order differential operators; notable exceptions are the works [9,10,12,15,17]. In this paper we consider a fourth order elliptic boundary value problem with boundary conditions that are not suitable for writing it as a second order system. The model problem we study here is a simplified version of the well-known Orr–Sommerfeld equation from hydrodynamics (cf. [3,12]), and reads: Find u ∈ C 4 ([0, 1]) such that

ε2 u (4) (x) − (α (x)u  (x)) + β(x)u (x) = f (x) in I = (0, 1) ,

*

Corresponding author. E-mail address: [email protected] (C. Xenophontos).

http://dx.doi.org/10.1016/j.apnum.2016.02.002 0168-9274/© 2016 IMACS. Published by Elsevier B.V. All rights reserved.

(1)

JID:APNUM AID:2992 /FLA

[m3G; v1.172; Prn:15/02/2016; 13:24] P.2 (1-17)

P. Panaseti et al. / Applied Numerical Mathematics ••• (••••) •••–•••

2

where α > 0, β ≥ 0 and f are given (sufficiently smooth) functions and ε ∈ (0, 1] is a given parameter (that can approach 0). The notation u (n) , n ∈ N, denotes the nth derivative of u with respect to x and the interval is taken as (0, 1) for simplicity. Equation (1) is supplemented with the following boundary conditions:

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

(2)

The rest of the paper is organized as follows: In Section 2 we construct a C 1 approximation from a finite dimensional subspace of H 2 , using a hierarchical basis from [13]. In Section 3 we obtain results on the regularity of the solution to (1), (2) that will be used for the approximation, which is presented Section 4. Section 5 contains the results of numerical computations and finally Section 6 contains some conclusions. Throughout the paper we will utilize the notation H k ( I ) to mean the usual Sobolev space of functions defined on I , whose 0, 1, . . . , k (generalized) derivatives belong to L 2 ( I ), with associated norm and seminorm ·k, I , |·|k, I , respectively; the Lebesgue spaces L p ( I ), 1 ≤ p ≤ ∞ (and their norms) are defined in the usual way and ·, · I denotes the usual L 2 ( I ) inner-product. We will also use the spaces





H k0 ( I ) = u ∈ H k ( I ) : u (k) |∂ I = 0 , k = 1, 2, where ∂ I denotes the boundary of I . Finally, the letters c , C (with or without subscripts) will denote generic positive constants independent of the solution, ε , and any discretization parameters. 2. Finite element formulation The variational formulation of (1), (2) reads: Find u ∈ H 02 ( I ) such that

B(u , v ) = F ( v ) ∀ v ∈ H 02 ( I ) ,

(3)

where

B (u , v ) =

1 



ε2 u  (x) v  (x) + α (x)u  (x) v  (x) + β(x)u (x) v (x) dx,

(4)

0

1 F (v ) =

f (x) v (x)dx.

(5)

0

We assume analyticity of the input data, i.e., there exist positive constants C α , C β , C f , γα , γβ , γ f independent of ∀ n = 0, 1, 2, . . .

   (n)  α 

L∞ (I )

    ≤ C α γαn n! , β (n) 

L∞ (I )

    ≤ C β γβn n! ,  f (n) 

L∞ (I )

≤ C f γ fn n! .

ε such that (6)

Coercivity of the bilinear form B (·, ·) holds in the energy norm

|u |2E , I := B(u , u ).

(7)

Existence and uniqueness follow from the Lax–Milgram lemma as usual. The discrete version of (3) reads: Find u N ∈ S N ⊂ H 02 ( I ) such that

B (u N , v ) = F ( v ) ∀ v ∈ S N ,

(8)

and we have

|u − u N | E , I ≤ C inf |u − v | E , I ∀ v ∈ S N . v∈S N

(9)

In order to define the space S N , let p ≥ 3 be an integer and introduce, for t ∈ I ST := (−1, 1), the four nodal basis functions

N 1 (t ) = N 3 (t ) =

1 4 1 4

(1 + t )2 (2 − t ) , N 2 (t ) = (t − 1) (1 + t )2 , N 4 (t ) =

1 4 1 4

(1 − t )2 (2 + t ) , (1 + t ) (1 − t )2 ,

as well as the p − 3 internal basis functions (cf. [13])

 N i (t ) =

2i − 5 2

 t η P i −3 (η)dηdx, i ≥ 5, −1 −1

(10)

JID:APNUM AID:2992 /FLA

[m3G; v1.172; Prn:15/02/2016; 13:24] P.3 (1-17)

P. Panaseti et al. / Applied Numerical Mathematics ••• (••••) •••–•••

3

Fig. 1. Hierarchical basis functions N i (x), i = 1, . . . , 7.

where P k (t ) is the Legendre polynomial of degree k. The properties of the Legendre polynomials allow us to obtain the following formula for N i :

N i (t ) = √



1 4i − 10

1 2i − 7

P i −5 (t ) −

2(2i − 5)

(2i − 7)(2i − 3)

P i −3 (t ) +

1 2i − 3



P i −1 (t ) ,

for i ≥ 5. Fig. 1 shows the plots of the first seven hierarchical basis functions. Now, let  p ( I ) denote the space of polynomials of degree ≤ p on I and note that



 p ( I ST ) = span N 1 (t ), . . . , N p +1 (t ) , p ≥ 3.

Let = {x0 , . . . , xm } be an arbitrary partition of I and set I j = (x j −1 , x j ). We note that the mapping x = Q j (t ) = t )x j −1 +

1 (1 + t )x j 2



1 (1 2



maps I ST to I j . With this in mind, and with ◦ denoting composition of functions, we define the space



1 S p ( ) := w ∈ H 02 ( I ) : w ◦ Q − ∈  p j ( I ST ) , j = 1, . . . , m , j

and choose as our subspace S N ⊂ H 02 ( I ), p

S N ≡ S 0 ( ) = S p ( ) ∩ H 02 ( I ).

(11)

We restrict our attention here to constant polynomial degree p for all elements, i.e. p j = p ∀ j, but clearly more general settings with variable polynomial degrees are possible. The following Spectral Boundary Layer mesh is essentially the minimal mesh that yields exponential convergence, as p → ∞. Definition 1 (Spectral Boundary Layer mesh). For mials by



S (κ , p ) :=

p

S 0 ( ); p S 0 ( );

κ > 0, p ∈ N and 0 < ε ≤ 1, define the spaces S (κ , p ) of piecewise polyno-

= {0, κ p ε , 1 − κ p ε , 1} if κ p ε < = { 0, 1 } if κ p ε ≥

1 2 1 2

The parameter κ is user specified and depends on the problem under consideration as well as the length scales of the boundary layers – we refer to [14] for a more detailed discussion on this issue and we note that in practice the value κ = 1 yields satisfactory results for most problems. We also note that the method we are considering is not a true hp FEM since the location and not the number of the elements changes; a more correct characterization would be a p version FEM on a moving mesh, but in order to be consistent with the bibliography we utilize the term hp FEM for our method. Obviously, additional refinement and/or using a true hp version would yield better results but at the cost of using more degrees of freedom. Remark 1. We will be interested in error estimates in the energy norm defined by (7), since one could not hope to obtain estimates in the H 2 norm that are uniform in ε . The use of more balanced norms has been investigated recently, see [2,11] and is something that we are currently considering for second order [6] as well as fourth order [18] problems.

JID:APNUM AID:2992 /FLA

[m3G; v1.172; Prn:15/02/2016; 13:24] P.4 (1-17)

P. Panaseti et al. / Applied Numerical Mathematics ••• (••••) •••–•••

4

We close this section with the following. Theorem 1. There exist constants C , γ > 0, independent of ε , such that the solution u of (3) satisfies

   (n)  u 

L∞ (I )

≤ C γ n max{nn , ε 1−n } ∀ n = 0, 1, 2, . . .

Proof. This follows from equations (2.8)–(2.12) in [17], in an inductive way.

2

3. Regularity results For the analysis, we will make the simplifying assumption that α (x) = β(x) = 1, in order to avoid technical difficulties that arise when non-constant coefficients are considered. The analysis extends to the variable coefficient case, as was done in [8] for second order systems. x We introduce the stretched variables x˜ = εx , xˆ = 1− ε and make the formal ansatz

u=



ε j {u j (x) + uˆ j (ˆx) + u˜ j (˜x)},

(12)

j =0

where u j , uˆ j , u˜ j are to be determined. This is done by inserting (12) into (1) and equating like powers of ε for the slow (i.e., x) and the fast (i.e., xˆ , x˜ ) variables. This procedure gives information on the unknown functions u j , uˆ j , u˜ j ; in particular, we find the following:

−u 0 + u 0 = f in I ,

(13)



−u 1 + u 1 = 0 in I ,

(14)



(4 ) −u j + u j = −u j −2 , in I , j (4 ) (2 ) (4 )

≥2

(15)

u˜ 0 = 0 , u˜ 1 − u˜ 1 = 0 in I ,

(16)

(2 )

(17)

u˜ j − u˜ j

= −u˜ j −2 , in I , j = 2, 3, . . .

(uˆ j satisfies analogous equations as (16)–(17)). The above problems are supplemented with appropriate boundary conditions in order for (2) to be satisfied; namely,

u j (0) = u j (1) = 0 , j ≥ 0,

(18)

˜

(19)



u˜ j (0) = 0 , u j (0) = −u j −1 (0) , j ≥ 1,

˜

lim u j (x) = 0 , j ≥ 1.

x→∞

(20)

Analogous conditions as (19), (20) are enforced at the right endpoint, involving uˆ j . Now, let M ∈ N and define

w M (x) :=

M

ε j u j (x) , u˜ M (˜x) :=

j =0

M

ε j u˜ j (˜x) , uˆ M (ˆx) :=

j =0

M

ε j uˆ j (ˆx),

(21)

j =0

r M := u − ( w M + u˜ M + uˆ M ).

(22)

This gives the decomposition

u = w M + u˜ M + uˆ M + r M ,

(23)

is the smooth part, u˜ M , uˆ M are the left and right boundary layer parts, respectively, and r M is the remainder.

where w M In what follows, we obtain regularity results for each one of these components; we begin with an auxiliary lemma that will be used for the regularity of u j . Lemma 2. Let g be an analytic function on I = [0, 1] satisfying, for some positive constants γ > 1, C g ,

   (n)  g 

L∞ (I )

≤ C g γ n n! ∀n ∈ N0 ,

(24)

and let v be the solution of the boundary value problem

− v  + v = g in I −

v (0) = g , v (1) = g

+

(25) (26)

JID:APNUM AID:2992 /FLA

[m3G; v1.172; Prn:15/02/2016; 13:24] P.5 (1-17)

P. Panaseti et al. / Applied Numerical Mathematics ••• (••••) •••–•••

5

for some constants g ∓ . Then,

   (n)  v 

L∞ (I )



C˜ {C g + | g − | + | g + |},

n = 0, 1 ˜C γ n−2 (n − 2)!{C g + | g − | + | g + |}, n ≥ 2

(27)

with C˜ > 1 a constant independent of n (but depending on γ ). Proof. The variational formulation of (25)–(26) gives the stability estimate

  v 

0, I

 +  v 0, I ≤  g 0, I + C | g − | + | g + | ,



with C ≤ 2 7/6. Sobolev’s imbedding theorem then gives (27) for n = 0 by choosing C˜ ≥ 1 appropriately. For n = 1, (27) follows again from Sobolev’s imbedding theorem and the fact that

  v  0, I ≤  g 0, I +  v 0, I ≤ 2 g 0, I + C | g − | + | g + | . We suppose now that (27) holds for n and we will show it for n + 2. We have from (25) that v (n+2) = v (n) − g (n) , hence

   (n+2)  v 

L∞ (I )

    ≤  v (n)  ≤ C˜ γ

L∞ (I )

n −2

    +  g (n) 

L∞ (I )

(n − 2)!{C g + | g − | + | g + |} + C g γ n n!

where (24) and the induction hypothesis were used. Thus

   (n+2)  v 

 −

L∞ (I )

+



n

≤ {C g + | g | + | g |}γ n! 1 +



γ2

≤ C˜ γ n n!{C g + | g − | + | g + |}, 

provided C˜ ≥ 1 − γ12

−1

. This completes the proof.

2

Lemma 3. There exist constants C , K 1 , K 2 > 0, independent of ε , such that u j defined by (13)–(15), satisfies

   (n)  u j 

j

L∞ (I )

≤ C K 1n K 2 j !n!, j = 0, 2, 4, . . . ∀ n ∈ N0 .

(28)

Proof. For j = 0, we use Lemma 2 to see that u 0 defined by (13), (18), satisfies

   (n)  u 0 

L∞ (I )

≤ C f γ fn n!,

(29)

since f satisfies (6). The desired result follows (for j = 0) by taking C = C f , K 1 = γ f in (28). For j = 1, the homogeneous boundary conditions (and right hand side) of (14) yield u 1 = 0. In fact, it is not hard to see that u 2 j +1 = 0, j = 0, 1, . . . . To show the desired result for n ≥ 2, we write

u j +2 = u j + r j +2 ,

(30)

where r j +2 satisfies

−r j +2 + r j +2 = u j ,

(31)



(32)



r j +2 (0) = −u j (0) , r j +2 (1) = −u j (1). The use of a Green’s function gives the representation

r j +2 (x) = −

  u j (0)  u j (1)  2− x x 1+ x 1− x e − e e − e − (e 2 − 1) (e 2 − 1)



1

⎧ ⎨

e 2(e 2 − 1) ⎩

2

e

−x

−e

x



1 0

e − y u j ( y )dy

JID:APNUM AID:2992 /FLA

[m3G; v1.172; Prn:15/02/2016; 13:24] P.6 (1-17)

P. Panaseti et al. / Applied Numerical Mathematics ••• (••••) •••–•••

6



+ e x − e −x

1

e y u j ( y )dy + (1 − e 2 )e −x

0

1 + (e 2 − 1)e x

e − y u j ( y )dy

e y u j ( y )dy

0

⎫ ⎬ ⎭

x

.

0

Thus

r j +2 L ∞ ( I ) ≤ C u j L ∞ ( I ) , where the constant C > 0 is independent of u j . The triangle inequality and (30) give

u j +2 L ∞ ( I ) ≤ C u j L ∞ ( I ) ,

(33)

with C > independent of u j . We now proceed as in the proof of Lemma 2 in [4] and we prove the following, stronger assertion: Let G ⊂ C be a complex neighborhood of [0, 1] and let

G δ = { z ∈ G : dist( z, ∂ G ) ≥ δ}. Then there exists a constant K > 0 such that for 0 < δ ≤ 1,

u j L ∞ (G δ ) ≤ C δ − j K j j !u 0 L ∞ (G ) , j = 0, 2, 4, . . .

(34)

where C is independent of j , δ and K . We will show (34) by induction on j, and thereafter the desired result will follow from Cauchy’s Integral Theorem for derivatives. Clearly (34) holds for j = 0, as u 0 may be extended to a holomorphic function (denoted again by u 0 ) due to its analyticity (cf. (29)). So, we assume (34) holds for j and we will show it for j + 2. For any κ ∈ (0, 1), x ∈ G δ , we may apply Cauchy’s Integral Theorem for derivatives, where the integration path is chosen as the circle of radius κ δ about x, to get together with (33),

u j +2 L ∞ (G δ ) ≤ C u j L ∞ (G δ ) ≤ C

2 2πκ δ 2π (κ δ)3

u j L ∞ (G (1−κ )δ ) .

By the induction hypothesis, we obtain

u j +2 L ∞ (G δ ) ≤ 2C (κ δ)−2 j ! K j ((1 − κ )δ)− j u 0 L ∞ (G )  −( j +2) j +2 ∞ ≤ Cδ K ( j + 2)!u 0 L (G ) 2 2 K

Choosing

κ=

1 j +2

2

κ (1 − κ ) j ( j + 1)( j + 2)

 .

yields

1

κ 2 (1 − κ ) j ( j + 1)( j + 2)

 =

j+2 j+1

 j +1

 = 1+

1

j

j+1

≤ e ∀ j ≥ 0,

hence the quantity in brackets above can by bounded by 1 if K 2 > 2e. This completes the proof.

2

We would like to obtain an analogous result to the above lemma, for the functions u˜ j defined by (16)–(17). To accomplish this, we need the following Lemma 4. Let F be an entire function satisfying, for some C F > 0, q ≥ j + 2,

   F  ≤ C F e −Re(z) (q + | z|) j /2 ∀ z ∈ C.

Let g 1 , g 2 ∈ C, g 2 = 0 and let w : (−1, ∞) → C be the solution of

w (4) − w (2) = F , w (0) = g 1 , w  (0) = g 2 , lim w  (x) = 0. x→∞

Then, w can be extended to an entire function (denoted again by w) which satisfies

  | w  ( z)| ≤ C C F (q + | z|) j /2+2 ( j /2 + 2)−1 + | g 2 | e −Re(z) ∀ z ∈ C, with C a constant independent of w and j.

JID:APNUM AID:2992 /FLA

[m3G; v1.172; Prn:15/02/2016; 13:24] P.7 (1-17)

P. Panaseti et al. / Applied Numerical Mathematics ••• (••••) •••–•••

7

Proof. The proof follows that of Lemma 7.3.6 in [5]. First note that v = w  satisfies

v (4) − v (2) = F  , v (0) = g 2 , lim v (x) = 0. x→∞

With the aid of a Green’s function, we have the representation

v ( z) = e −z

∞

y F  ( y )dy +

1 −z e 2

∞

0

∞ −

e − y F  ( y )dy −

0

y F  ( y )dy + z

z

∞

1 −z e 2

z

e y F  ( y )dy

0

F  ( y )dy −

1 2

∞ ez

z

e − y F  ( y )dy + g 2 e −z .

(35)

z

We estimate each integral above, as was done in the proof of Lemma 7.3.6 in [5]. We have for the last integral in (35)

 ∞  ∞       ∞    z    −y  −t  e   e F ( y )dy  =  e F (t + z)dt  ≤ e −t F  (t + z) dt      z

0

0

∞    ≤ e −t C F e −(z+t ) (q + | z| + |t |) j /2  dt 0

∞   −2t  e (q + | z| + |t |) j /2  dt

≤ CFe

−Re( z)

≤ CFe

−Re( z) (q+|z|) −1− j /2

0

e

2

(1 + j /2, q + | z|),

where (·, ·) denotes the incomplete Gamma function (see, e.g., [5] p. 245). There holds

| (α , ξ )| ≤

e −ξ ξ α

|ξ | − α0

, α0 = max{α − 1, 0} , Re(ξ ) ≥ 0 , |ξ | > α0 .

Hence,

 ∞     1+ j /2  z  −y e  ≤ C F e −Re(z) e (q+|z|) 2−1− j /2 e −(q+|z|) (q + | z|) e F ( y ) dy   q + | z | − z − j /2   z

≤ C F e −Re(z)

(q + | z|)1+ j /2 (q + | z|)2+ j /2 ≤ C F e −Re(z) . 2 + j /2 2 + j /2

We then consider the integral

    z 1  −z  y  −Re( z) e  e F ( y )dy  ≤ e C F (q + t | z|) j /2 | z|e −Re(tz) e Re(tz) dt    0

0

≤ C F e −Re(z)

1 j /2 + 1

≤ 2C F e −Re(z) Next we note that the integral

∞ 0

{(q + | z|) j /2+1 − q j /2+1 }

(q + | z|)2+ j /2 . 2 + j /2

e − y F  ( y )dy is the same as the last integral in (35) with z = 0. Thus,

    ∞ 2+ j /2  −z  −y  e  ≤ C F e −Re(z) (q + | z|) e F ( y ) dy .   2 + j /2   0

∞

The remaining terms in (35) are bounded as follows: The last term g 2 e −z , as well as the term e −z 0 y F  ( y )dy, satisfy the desired bound, as they are merely multiples of e −z . Combining the two terms that remain in (35) we have, using again the incomplete function and its properties,

JID:APNUM AID:2992 /FLA

[m3G; v1.172; Prn:15/02/2016; 13:24] P.8 (1-17)

P. Panaseti et al. / Applied Numerical Mathematics ••• (••••) •••–•••

8

∞   ∞            ( z − y ) F  ( y )dy  = − t F  (t + z)dt          z

0

≤ CFe

−Re( z)

∞

|t |e −Re(t ) (q + |t | + | z|) j /2 dt

0

≤ C F e −Re(z)

∞

e −Re(t ) (q + |t | + | z|) j /2+1 dt

0

≤ CFe

−Re( z) (q+|z|)

e

(2 + j /2, q + | z|)

(q + | z|)2+ j /2 ≤ C F e −Re(z) j /2 + 1 /2 ≤ C C F e −Re(z)

(q + | z|)2+ j /2 , j /2 + 2

for a positive constant C , independent of j. Combining all of the estimates above, we get the desired result.

2

Lemma 5. The functions u˜ j given by (16)–(17), are entire and there exist constants C , K , γ > 0, independent of ε , such that for any n = 1, 2, . . .

|u˜ j (˜x)| ≤ C e −|˜x| K n γ j j j . (n)

(36)

Proof. We begin with a claim: there exist constants C , γ > 0, independent of ε , such that

|u˜ j (˜x)| ≤ C

γj j!

( j /2 + 2 + |˜x|) j /2+2 e −|˜x| .

(37)

To show the claim, we note that for u˜ 0 , (37) holds trivially. We also check, by direct calculation, that

u˜ 1 (˜x) = u 0 (0)(e −|˜x| − 1), hence (37) is satisfied for j = 1. Assuming now that (37) holds for j > 1, we will show it for j + 2. The function u˜ j +2 , satisfies (17), (19), (20):

u˜ j +2 − u˜ j +2 = −u˜ j in I , u˜ j +2 (0) = −u j +1 (0), lim u˜ j +2 (x) = 0, (4 )

(2 )

x→∞

with u˜ j satisfying (37), by the induction hypothesis. Lemma 4 then gives the desired bound on u˜ j +2 . To complete the proof, we comment on how (37) leads to (36). Cauchy’s Integral Theorem for derivatives allows us to infer from (37),

|u˜ (jn) (˜x)| ≤ C e −|˜x| ≤ C e −|˜x|

n!

(n

+ 1)n n!

(n

+ 1)n

γj j!

γj j!

( j /2 + 2 + |˜x|) j /2+2 en ( j /2 + 2 + n) j /2+2 en ,

where Lemma 4.8 of [7] was used. Using the observation that ( j /2 + 2 + n) j /2+2 ≤ ( j /2)( j /2+2) (1 + 2(n + 2)/ j ) j /2+2 ≤ ( j /2)( j /2+2) e (n+2) , the proof is completed. 2 Combining Lemmas 3 and 5, we have the following result for the components appearing in the decomposition (23). Theorem 6. Assume (6) holds. Then there exist constants C , K , K 1 , K 2 , q, γ > 0 independent of ε ∈ (0, 1] such that the solution u of (1), (2) can be written as u = w M + u˜ M + uˆ M + r M , with

   (n)   w M  ≤ C K 1n n! ∀ n = 0, 1, 2, . . . , 0, I      (n)   (n)  u˜ M (x) + uˆ M (x) ≤ C K n ε 1−n e −dist(x,∂ I )/ε ∀n ∈ N,   r M  L ∞ (∂ I ) + r M  L ∞ (∂ I ) + |r M | E , I ≤ C e −q/ε ,

provided ε K 2 (1 + γ ) M < 1.

(38) (39) (40)

JID:APNUM AID:2992 /FLA

[m3G; v1.172; Prn:15/02/2016; 13:24] P.9 (1-17)

P. Panaseti et al. / Applied Numerical Mathematics ••• (••••) •••–•••

9

Proof. First we show (38). We have from (21) and Lemma 3,

   (n)  w M 

0, I



M





ε j u (jn) 

0, I

j =0

≤ C K 1n n!

≤C

M

ε j K 1n K 2j j !n!

j =0

M

(ε K 2 j ) j

j =0

∞ ≤ C K 1n n! (ε K 2 M ) j j =0



C K 1n n!,

since the sum is a converging geometric series due to the assumption ε K 2 M < 1. (n) In order to show (39), we utilize Lemma 5; we will only show the procedure for u˜ M (˜x) as the other one is similar. We have

|u˜ (Mn) (˜x)| ≤

M

ε j |u˜ (jn) (˜x)| ≤

j =0

M

ε j e−˜x K n γ j j j

j =1

≤ K n e −˜x

M

−1 M

j =1

j =0

(εγ M ) j ≤ K n e −˜x

(εγ M ) j +1

∞ ≤ C ε K n e −˜x (εγ M ) j j =0 n −˜x

≤ CεK e

,

since the sum above is a converging geometric series under the assumption that It remains to show (40). To this end note that

εγ M < 1. Equation (39) follows.

⎛ ⎞ M M M r M (0) = u (0) − ⎝ ε j u j (0) + ε j u˜ j (0) + ε j uˆ j (1/ε)⎠ j =0

=−

M

j =0

j =0

ε j uˆ j (1/ε),

j =0

since u (0) = u j (0) = u˜ j (0) = 0. Hence, by the analogous lemma to Lemma 5 for uˆ j , we have

|r M (0)| ≤

M

ε j |uˆ j (1/ε)| ≤

j =0

M

ε j C K e−1/ε γ j j !

j =0

≤ C e −1/ε

M

(ε M γ ) j ≤ C e −1/ε

j =0

≤ Ce

−q/ε



(ε M γ ) j

j =0

,

for some positive q ∈ R, independent of

ε and clearly bounded away from 0. Similarly, ⎞ M M M r M (0) = u  (0) − ⎝ ε j u j (0) + ε j u˜ j (0) + ε j uˆ j (1/ε)⎠ ⎛

⎛ = u  (0) − ⎝

j =0

j =0

M

M

ε j u j (0) +

j =0

=

M j =0

ε j (−uˆ j (1/ε))

j =0

j =0

ε j+1 u˜ j+1 (0) +

M j =0



ε j uˆ j (1/ε)⎠

JID:APNUM AID:2992 /FLA

[m3G; v1.172; Prn:15/02/2016; 13:24] P.10 (1-17)

P. Panaseti et al. / Applied Numerical Mathematics ••• (••••) •••–•••

10

since u  (0) = 0, u˜ 0 (0) = 0, u˜ j (0) + u j −1 (0) = 0. Hence,

|r M (0)| ≤

M

ε j |uˆ j (1/ε)| ≤

j =0

M

ε j C K e−q/ε γ j j ! ≤ C e−q/ε .

j =0

In the same way we obtain analogous results for |r M (1)|, |r M (1)|, thus

  r M  L ∞ (∂ I ) + r M  L ∞ (∂ I ) ≤ e −q/ε . 4

d Now, let L := ε 2 dx 4 −

d2 dx2

L (u − w M ) = f −

+ Id, and consider the following: M

ε j Lu j = f −

j =0

M

ε j (ε2 u (j4) − u j + u j ),

j =0

with u j satisfying (13)–(15). After some calculations we find that

L ( u − w M ) = ε M +2 u M , (4 )

hence by Lemma 3 we have

 L (u − w M ) L ∞ ( I ) ≤ C ε M +2 K 2M M ! ≤ C (ε K 2 M ) M +2 .

(41)

We next consider L u˜ M which takes the form (in the variable x˜ ),

L u˜ M =

M

ε j (ε−2 u˜ (j4) − ε−2 u˜ (j2) + u˜ j ),

j =0

with u˜ j satisfying (16)–(17). These equations allow us to obtain

L u˜ M =

M

ε j (−ε−2 u˜ j−2 − ε2 u˜ j ) = ε M +2 u˜ M .

j =2

Thus, by Lemma 5 we have

   L u˜ M 

≤ C ε M +2 e −˜x γ M ( M ) M ≤ C ε 2 e −˜x (εγ M ) M .

L∞ (I )

(42)

Combining the analogous result for L uˆ M with (41) and (42) we have

   Lr M  L ∞ ( I ) =  L (u − w M − u˜ M − uˆ M )∞, I     ≤  L (u − w M )∞, I +  L u˜ M ∞, I +  L uˆ M ∞, I ≤ C ε 2 (ε K 2 (γ + 1) M ) M .

We have shown that r M has exponentially small values at the endpoints of [0, 1] and Lr M is uniformly bounded by an exponentially small quantity on the interval (0, 1). By, e.g., stability, we have the desired result. 2 4. Error estimates Before we prove our main theorem, we present the following approximation result. Lemma 7. There exists a linear operator I p : H 2 ( I ST ) →  p ( I ST ), p ≥ 3, with the following property: (k)

I p u (±1) = u (k) (±1) , k = 0, 1.

 Furthermore, if u ∈ C ∞ I ST , then

  u − I p u 2

0, I ST

    2  u − Ip u 

0, I ST

    2  u − Ip u 

0, I ST

 1 ( p − s)!   ( s +1 )  2 u , ∀ s = 0, 1 , . . . , p ,   0, I ST p 4 ( p + s)!  1 ( p − s)!   ( s +1 )  2 ≤ 2 , ∀ s = 0, 1 , . . . , p , u  0, I ST p ( p + s)!  ( p − s)!   ( s +1 )  2 ≤ , ∀ s = 0, 1 , . . . , p .  u 0, I ST ( p + s)! ≤

(43)

(44) (45) (46)

JID:APNUM AID:2992 /FLA

[m3G; v1.172; Prn:15/02/2016; 13:24] P.11 (1-17)

P. Panaseti et al. / Applied Numerical Mathematics ••• (••••) •••–•••

11

Proof. This result is an extension of Corollary 3.15 in [13] (from H 1 to H 2 ) and our proof follows the same steps – see also Lemma 5.3 in [20]. Let I p u be the Legendre series of u  , truncated after the Legendre polynomial L p −1 . From Lemma 3.9 in [13] we have

   u − I  u 2 p

0, I ST

=

∞ i= p

2 2i + 1

|ai |2 ,

(47)

where

ai =

1

2i + 1 2

u  (ξ ) L i (ξ )dξ.

(48)

−1

Define





I p u (ξ ) =

I p u (η)dη + u  (−1)

(49)

I p u (η)dη + u (−1).

(50)

−1

ξ I p u (ξ ) = −1

Then it is obvious that I p u (−1) = u  (−1) and I p u (−1) = u (−1). Moreover, from (48) we have that

a0 =

1 2

1

1

u  (η)dη =

2

−1

(u  (1) − u  (−1)).

But we also have 

1



I p u (1) − I p u (−1) =

I p u (η)dη = 2a0 ,

−1

from the fact that I p is a series of orthogonal Legendre polynomials. The two equations above give us

u  (1) − I p u (1) − (u  (−1) − I p u (−1)) = 0, from which we see that u  (1) = I p u (1). It remains to show that I p u (1) = u (1) in order to establish (43). We have from (50), and the fact that I p is a series of orthogonal Legendre polynomials,

1 I p u (1) =

I p u (η)dη + u (−1)

−1

=

⎧ 1 ⎨η

−1





−1

1 η

p

=

i =0

⎫ ⎬ I p u (ξ )dξ + u  (−1) dη + u (−1) ⎭

ai

P i (ξ )dξ dη + 2u  (−1) + u (−1)

−1 −1

= 2a0 +

p i =1

= 2a0 −

2a1 3

1 ai −1

1 2i + 1

{ P i +1 (η) − P i −1 (η)} dη + 2u  (−1) + u (−1)

+ 2u  (−1) + u (−1)

= u  (1) + u  (−1) −

2a1 3

+ u (−1).

From (48) we have, using integration by parts,

(51)

JID:APNUM AID:2992 /FLA

[m3G; v1.172; Prn:15/02/2016; 13:24] P.12 (1-17)

P. Panaseti et al. / Applied Numerical Mathematics ••• (••••) •••–•••

12

a1 =

=

1

3 2

−1



3⎨  ξ u (ξ )dξ = u (1) + u  (−1) − 2⎩ 

1

⎫ ⎬



u (ξ )dξ

−1



3  I (1) + I p (−1) − u (1) + I p (−1) . 2 p

(52)

Substituting (52) into (51) we get

I p u (1) = u  (1) − u  (−1) − (I p u (1) +

+ I p u (−1) − u (1) + I p u (−1)) + 2u  (−1) + u (−1) = u (1), and this establishes (43). Equation (46) follows from (47) and Lemma 3.10 in [13], so let us show (45). Observe that

I p u (x) − u  (x) =

x ∞

ai P i (t )dt =

−1 i = p +1

x



ai ψi (x),

i = p +1

where ψi (x) = −1 P i (t )dt. Integrating the Legendre differential equation, we see that

ψi (x) = −

1 i ( i + 1)

(1 − x2 ) P i (x),

hence

1

1 1 − x2

−1

ψi (x)ψ j (x)dx =

δi j i (i + 1)(2 j + 1)

,

by (3.3.9) in [13], where δi j is the Kronecker delta. Therefore,

1

   I u (x) − u  (x)2 dx ≤ p

−1

1 ≤ −1

1

⎧ ∞ ⎨

1 − x2 ⎩

i = p +1

1

1 1−

−1

x2

   I u (x) − u  (x)2 dx p

⎫2 ∞ ⎬ 2 ai ψi (x) dx ≤ |ai |2 ⎭ i (i + 1)(2i + 1) i = p +1

1 ∞ ( p − s)! ( p − s)! 1 1 2 2 (i + p )! ≤ |u (s+1) (x)|2 dx, |ai | ≤ ( p + s)! i (i + 1) (2i + 1) (i − p )! ( p + s)! p 2 i = p +1

−1

where (3.3.8) from [13] was used in the last step. This establishes (45). Finally, for (44) we have in a similar fashion

I p u (x) − u (x) =

 x η ∞

ai P i (ξ )dξ dη =

−1 −1 i = p +2

∞ i = p +2



2 2i + 3

ai N i (x),

with N i (x) given by (10). After some calculations (see [20]), we find that



1 N i (x) N j (x)dx −1

≤ C i13 =0

2 2i +1

otherwise

So we have

1 |I p u (x) − u (x)|2 dx ≤ −1

if j = i , i + 2, or i − 2

⎧ 1 ⎨ ∞ 

−1



i = p +2

⎫2 ∞ ⎬ 1 2 dx ≤ |ai |2 ai N i (x) 4 ⎭ 2i + 3 i 2i + 1 2

i = p +2

1 ∞ ( p − s)! 1 2 ( p − s)! 1 2 (i + p )! ≤C |ai | |u (s+1) (x)|2 dx. ≤C ( p + s)! p 4 2i + 1 (i − p )! ( p + s)! p 4 i = p +2

This completes the proof.

2

−1

JID:APNUM AID:2992 /FLA

[m3G; v1.172; Prn:15/02/2016; 13:24] P.13 (1-17)

P. Panaseti et al. / Applied Numerical Mathematics ••• (••••) •••–•••

13

We now present our main result: Theorem 8. Let u be the solution to (3) and let u N ∈ S (κ , p ) be the solution of (8) with S (κ , p ) given by Definition 1. Then there exist constants κ0 , C , σ > 0 depending only on the input data α , β, f , such that for any 0 < κ ≤ κ0

|u − u N | E , I ≤ C e −σ κ p . Proof. The proof is separated in two cases: Case I: κ p ε ≥ 1/2 In this case, often called the “asymptotic range of p”, the mesh consists of one element, = {0, 1}. Note that the assumption κ p ε ≥ 1/2 gives (2κε )−1 ≤ p. The solution u satisfies

   (n)  u 

L∞ (I )

≤ C γ n max{nn , ε 1−n } ∀ n = 0, 1, 2, . . .

(53)

By Lemma 7, there exists a polynomial I p u ∈  p , p ≥ 3, such that (43)–(46) hold. This gives

|u

− I p u |2E , I

( p − s)! ≤C ( p + s)!

#

1 p4

+

1 p2

$ 2   + ε  u ( s +1 )  . 2

0, I

Choosing s = λ p ∈ N for some λ ∈ (0, 1) to be selected shortly, (53) gives

   2   ( s +1 )  2  u  = u (λ p +1)  ≤ C γ λ p +1 max{(λ p + 1)λ p +1 , ε (−λ p ) } 0, I

0, I

≤ Cγ

λ p +1

(λ p + 1)λ p +1 ,

κ ≤ λ/2. Thus, using the fact that (see, e.g. Lemma 3.1 in [19]) &p % (1 − λ)(1−λ) ( p − λ p )! p −2 λ p e 2 λ p +1 , ≤ ( p + λ p )! (1 + λ)(1+λ)

provided

we have

(54)

# $ ( p − λ p )! 1 1 2 + + ε C γ 2(λ p +1) (λ p + 1)2(λ p +1) ( p + λ p )! p 4 p2  p (1 − λ)(1−λ) ≤C p −λ p e λ p +1 γ 2(λ p +1) (λ p + 1)2(λ p +1) (1 + λ)(1+λ)  p (1 − λ)(1−λ) ≤C p 2 (γ e )λ p +1 (λ + 1/ p )λ p (1 + λ)(1+λ)  p (1−λ) 2 (1 − λ) λ ≤ Cp (γ e λ) . (1 + λ)(1+λ)

|u − I p u |2E , I ≤ C

Choosing λ ≤ (γ e )−1 ∈ (0, 1) gives us

|u − I p u | E , I ≤ C e −σ κ p , where

σ := | ln β|, β :=

(1−λ)(1−λ) . We note that (1+λ)(1+λ)

κ in the definition of the Spectral Boundary Layer Mesh must satisfy κ ≤

(γ e )−1 /2. The desired result then follows from (9). Case II: κ p ε < 1/2 In this case the mesh consists of three elements, := I 1 ∪ I 2 ∪ I 3 := {0, κ p ε , 1 − κ p ε , 1} and we decompose the solution u as in Theorem 6:

u = w M + u˜ M + uˆ M + r M . Since the remainder r M above is small (by Theorem 6), it will not be approximated. So we only need to consider w M , uˆ M and u˜ M . For w M , we have by Theorem 7 that there exists I p w M such that

   1 1 ( p − s)!   (s+1) 2 | w M − I p w M | E , I ≤ ε 2 + 2 + 4 w M  0, I ( p + s)! p p

JID:APNUM AID:2992 /FLA

[m3G; v1.172; Prn:15/02/2016; 13:24] P.14 (1-17)

P. Panaseti et al. / Applied Numerical Mathematics ••• (••••) •••–•••

14

˜ p, for some λ˜ ∈ (0, 1) to be selected shortly, such that By (54) and Theorem 6 we have, with s = λ

( p − λ˜ p )! λ˜ p +1 K (λ˜ p + 5)! ( p + λ˜ p )! % &p ˜ (1−λ˜ ) (1 − λ) ˜ ˜ ˜ ˜ p + 5)2λ˜ p +5 p −2λ p e 2λ p +1 K λ p +1 (2λ ≤C ˜) ( 1 + λ ˜ (1 + λ) % &p ˜ (1−λ˜ ) 5 ˜ ˜ 5 (1 − λ) ˜ ≤ C (2λ p + 5) (e K )2λ p +1 (2λ˜ + )2λ p ˜) ( 1 + λ ˜ p (1 + λ) % & p ˜ (1−λ˜ ) (1 − λ) ˜ ≤ Cp 5 (2e λ˜ K )2λ p +1 ˜ (1+λ˜ ) (1 + λ)

| w M − I p w M |2E , I ≤ C

˜ ≤ (2K e )−1 ∈ (0, 1) gives us Choosing λ

| w M − I p w M | E , I ≤ Cp 5 e −σ κ p , where

˜ β˜ := σ := | ln β|,

˜ ˜ (1−λλ) (1−λ) 5 ˜ . The term p can be absorbed into the exponential by adjusting the constants C and ˜ (1+λ) (1+λ)

σ.

Next, we consider u˜ M (only, since uˆ M is similar), keeping in mind that κ p ε < 1/2. We will construct two separate approximations for u˜ M , one on the interval I 1 = (0, κ p ε ) and one on I ∗ := I 2 ∪ I 3 . By Theorem 7, there exists I p u˜ M such that

   1 1 ( p − s)!   ( s +1 )  2 |u˜ M − I p u˜ M |2E , I 1 ≤ (κ p ε )2s ε 2 + 2 + 4 . u˜ M  0, I 1 ( p + s)! p p

ˆ p for some λˆ ∈ (0, 1) and using (54) gives Choosing s = λ ˆ

|u˜ M − I p u˜ M |2E , I 1 ≤ C (κ p ε )2λ p

κ p ε ( p − λˆ p )! 2λˆ p +2 2 ˆ K ε ε−2λ p−2 e−dist(x,∂ I )/ε dx ( p + λˆ p )! 0

ˆ

≤ C (κ p )2λ p

ˆ (1 − λˆ p )(1−λ p )

ˆ

ˆ

ˆ

p −2 λ p e 2 λ p +1 K 2 λ p +2 κ p ε

(1 + λˆ p )(1+λˆ p ) % &p ˆ ˆ (1−λ) ˆ p (1 − λ) ˆ ˆ 2λ ≤ C (κ ) e 2 λ p +1 K 2 λ p +2 ˆ ˆ (1+λ) (1 + λ) % &p ˆ (1 − λˆ )(1−λ) ˆ 2λ ≤C (κ e K ) ˆ (1 + λˆ )(1+λ) ≤ C e −σˆ κ p , ˆ ˆ βˆ := (1−λ) where σˆ := | ln β|, , provided ˆ (1+λ)

κ ≤ (e K )−1 . Hence we should choose κ ≤ e−1 min{ K −1 , γ −1 /2}.

Now, on I ∗ we approximate u˜ M by its linear interpolant I1 u˜ M and we have

|u˜ M − I1 u˜ M |2E , I ∗ ≤ |u˜ M |2E , I ∗ + |I1 u˜ M |2E , I ∗ . Note that the definition of the energy norm (7) and (39) give

1 |u˜ M |2E , I ∗ =

2 2 2 {ε 2 u˜ M + u˜ M + u˜ M }dx

κ pε

≤ C e −κ p . The term |I1 u˜ M |2E , I ∗ satisfies the same bound since

|I1 u˜ M |2E , I ∗ ≤ C u˜ M (κ p ε ) ≤ C e −κ p , so, combining the bounds on I 1 and I ∗ we have that

JID:APNUM AID:2992 /FLA

[m3G; v1.172; Prn:15/02/2016; 13:24] P.15 (1-17)

P. Panaseti et al. / Applied Numerical Mathematics ••• (••••) •••–•••

15

Fig. 2. Energy norm convergence for the hp version (constant coefficients BVP).

Fig. 3. Maximum norm convergence for the hp version (constant coefficients BVP). Left: Error in u. Right: Error in u  .

|u˜ M − I p u˜ M |2E , I ≤ C e −κ p . The desired result follows from (9).

2

In the next section we illustrate our main result through two numerical experiments. 5. Numerical results 5.1. Constant coefficients First we present the results of numerical computations for the BVP (1), (2) when the data is α = β = f = 1; an exact solution is available, hence our reported results are reliable. Moreover, the theory developed in the previous sections covers this case. Fig. 2 shows the percentage relative error in the energy norm,

Error := 100 ×

|u − u N | E , I , |u | E , I

(55)

versus the number of degrees of freedom, N, in a semi-log scale. As N is increased the error curves become straight, something that verifies the exponential convergence of the proposed method. Moreover, as ε → 0 the method not only does not deteriorate, but it actually performs better. This suggests that in

JID:APNUM AID:2992 /FLA

[m3G; v1.172; Prn:15/02/2016; 13:24] P.16 (1-17)

P. Panaseti et al. / Applied Numerical Mathematics ••• (••••) •••–•••

16

Fig. 4. Energy norm convergence for the hp version (variable coefficients BVP).

Fig. 5. Estimated maximum norm convergence for the hp version (variable coefficients BVP). Left: Error in u. Right: Error in u  .

the error estimate of Theorem 8, there is a positive power of ε present, something that, we believe, is due to the fact that the problem has constant coefficients and the right hand side belongs to the subspace (see [14] for more details on this for second order scalar problems). The use of a more balanced norm (see Remark 1) would eliminate this phenomenon. Fig. 3 shows the maximum norm convergence, both for u and its derivative. The robustness of the proposed method in this norm as well, is readily visible. Moreover, the positive power of ε does not seem to be present in the derivative error, something that is consistent with our observations in Remark 1. 5.2. Variable coefficients We now consider the BVP (1), (2) when the data is α (x) = e −x , β(x) = 0, f (x) = e −x +1 ; an exact solution is not available so for our computations we use a reference solution obtained with a high number of degrees of freedom (of the order 103 ). This case is not covered by the theory presented in the previous sections but, as already mentioned, extending it is possible and as we will see the method performs equally well. This is what Fig. 4 shows, i.e. the robustness and exponential convergence of the method. We note that the method performs better as ε gets smaller, in this case as well, even though the data does not belong to the finite dimensional subspace. Finally, in Fig. 5 we show the convergence of the solution and its derivative, in the maximum norm, for various values of ε . The robustness of the proposed method is, once again, readily visible, as is the absence of the positive power of ε in the derivative error. 2

JID:APNUM AID:2992 /FLA

[m3G; v1.172; Prn:15/02/2016; 13:24] P.17 (1-17)

P. Panaseti et al. / Applied Numerical Mathematics ••• (••••) •••–•••

17

6. Conclusions We considered fourth order singularly perturbed BVPs in one-dimension and the approximation of their solution by the hp version of the FEM. We constructed a C 1 -conforming method that, when used in conjunction with the Spectral Boundary Layer mesh, yields a robust approximation that converges exponentially, when the error is measured in the energy norm. Through numerical experiments we verified the theoretical results and illustrated how well the method performs even in the cases not covered by our theory. Acknowledgement The last author would like to thank Prof. Serge Nicaise (Univ. Valenciennes) for helpful discussions and suggestions regarding this work. References [1] N.S. Bakhvalov, Towards optimization of methods for solving boundary value problems in the presence of boundary layers, Zh. Vychisl. Mat. Mat. Fiz. 9 (1969) 841–859 (in Russian). [2] S. Franz, H.-G. Roos, Error estimation in a balanced norm for a convection–diffusion problem with two different boundary layers, Calcolo 51 (2014) 423–440. [3] N. Madden, M. Stynes, G.P. Thomas, On the application of robust numerical methods to a complete-flow wave-current model, in: International Conference on Boundary and Interior Layers, BAIL 2004, Toulouse, France, 2004. [4] 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. [5] J.M. Melenk, Finite Element Methods for Singular Perturbations, Springer-Verlag, 2002. [6] J.M. Melenk, C. Xenophontos, Robust exponential convergence of hp-FEM in balanced norms for singularly perturbed reaction–diffusion equations, Calcolo (2016), http://dx.doi.org/10.1007/s10092-015-0139-y, in press. [7] J.M. Melenk, C. Xenophontos, L. Oberbroeckling, Analytic regularity for a singularly perturbed system of reaction–diffusion equations with multiple scales: proofs, arXiv:1108.2002 [math.NA], 2011. [8] J.M. Melenk, C. Xenophontos, L. Oberbroeckling, Robust exponential convergence of hp-FEM for singularly perturbed systems of reaction–diffusion equations with multiple scales, IMA J. Numer. Anal. 33 (2013) 609–628. [9] Z.-X. Pan, The difference and asymptotic methods for a fourth order equation with a small parameter, in: BAIL IV Proceedings, Book Press, Dublin, 1986, pp. 392–397. [10] H.-G. Roos, A uniformly convergent discretization method for a singularly perturbed boundary value problem of the fourth order, Rev. Res., Fac. Sci., Math. Ser., Univ. Novi Sad 19 (1989) 51–64. [11] H.-G. Roos, M. Schopf, Convergence and stability in balanced norms for finite element methods on Shishkin meshes for reaction–diffusion problems, ZAMM 95 (2015) 551–565. [12] H.-G. Roos, M. Stynes, A uniformly convergent discretization method for a fourth order singular perturbation problem, Bonner Math. Schr. 228 (1991) 30–40. [13] C. Schwab, p- and hp-Finite Element Methods, Oxford Science Publications, 1998. [14] C. Schwab, M. Suri, The p and hp versions of the finite element method for problems with boundary layers, Math. Comput. 65 (1996) 1403–1429. [15] G.I. Shishkin, A difference scheme for an ordinary differential equation of the fourth order with a small parameter at the highest derivative, Differ. Equ. 21 (1985) 1734–1742 (in Russian). [16] G.I. Shishkin, Grid approximation of singularly perturbed boundary value problems with a regular boundary layer, Sov. J. Numer. Anal. Math. Model. 4 (1989) 397–417. [17] G. Sun, M. Stynes, Finite-element methods for singularly perturbed high order elliptic two point boundary value problems I: reaction–diffusion-type problems, IMA J. Numer. Anal. 15 (1995) 117–139. [18] P. Constantinou, C. Varnava, C. Xenophontos, An hp FEM for fourth order singularly perturbed problems, Numer. Algorithms (2016), http://dx.doi.org/ 10.1007/s11075-016-0108-9, in press. [19] C. Xenophontos, L. Oberbroeckling, On the finite element approximation of systems of reaction–diffusion equations by p /hp methods, J. Comput. Math. 28 (2010) 386–400. [20] A. Zouvani, The hp finite element method for singularly perturbed 4th order problems, Masters thesis, Department of Mathematics and Statistics, University of Cyprus, 2012.