- Email: [email protected]
Parameter uniform numerical methods for singularly perturbed elliptic problems with parabolic boundary layers
Applied Numerical Mathematics 58 (2008) 1761–1772 www.elsevier.com/locate/apnum
Parameter uniform numerical methods for singularly perturbed elliptic problems with parabolic boundary layers E. O’Riordan a,∗ , G.I. Shishkin b,1 a School of Mathematical Sciences, Dublin City University, Ireland b Institute of Mathematics and Mechanics, Russian Academy of Sciences, Ekaterinburg, Russia
Available online 22 November 2007
Abstract Parameter-uniform numerical methods for a singularly perturbed elliptic problem with parabolic boundary layers in the solution are analyzed. A priori parameter explicit bounds on the solution and its derivatives are obtained using a suitable decomposition of the solution into regular and layer components. Based on this decomposition a numerical method is constructed and analyzed. © 2007 IMACS. Published by Elsevier B.V. All rights reserved. MSC: 65N15 Keywords: Parabolic boundary layers; Elliptic equation; Piecewise-uniform mesh
1. Introduction Hemker [3] formulated a singular perturbation test problem, which can be viewed as a simple mathematical model of flow past a circular obstacle. This problem also acts as a mathematical challenge to numerical analysts interested in designing parameter-uniform numerical methods for singularly perturbed differential equations. In this paper, we examine a related simpler problem, where the circular object is absent and we confine the discussion to problems where no interior parabolic layers are present in the solution. In this paper, we will examine numerical methods for the problem Lu ≡ −εu + aux = f, u = g,
(x, y) ∈ Ω = (0, 1)2 ,
(x, y) ∈ ∂Ω,
a > α > 0,
∞
a ∈ C (D),
(1a) (1b)
f ∈C
5,γ
(D),
Ω ⊂ D.
(1c)
Asymptotic expansions of the constant coefficient version of this problem [10] indicate the complicated analytical nature of the solution, which will in general contain a regular exponential boundary layer, two parabolic (or character* Corresponding author.
E-mail addresses: [email protected] (E. O’Riordan), [email protected] (G.I. Shishkin). 1 The research of the second author was supported by the Boole Centre for Research in Informatics at the National University of Ireland in
Cork, by the Mathematics Applications Consortium for Science and Industry in Ireland (MACSI) under the Science Foundation Ireland (SFI) Mathematics Initiative and by the Russian Foundation for Basic Research under grant No. 07-01-00729. 0168-9274/$30.00 © 2007 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2007.11.003
1762
E. O’Riordan, G.I. Shishkin / Applied Numerical Mathematics 58 (2008) 1761–1772
istic) boundary layers and several types of corner layer functions in the vicinity of the inflow and outflow corners. To analyze the convergence of any set of numerical approximations, we require bounds on the solution and its derivatives. In the case of constant coefficients, Kellogg and Stynes [4] supply these parameter explicit bounds on the derivatives of the solution in the general case of g ∈ / C 0 (∂Ω). Again, these bounds on the derivatives display the non-trivial nature of the solution even in the case of constant coefficients. Roos [9] also examines a constant coefficient version of problem (1), where a decomposition of the solution into a regular component and several layer components is given. In [9] bounds on the derivatives of the layer components are established by first deriving derivative bounds on the boundaries and then extending these throughout the domain by simply differentiating the differential equation. It is not easy to extend this argument to variable coefficient problems. In this paper, we examine the case of variable coefficients. However, here we impose sufficient regularity and compatibility on the data so that the amplitude of the inflow corner layer functions are negligible. u = g = 0,
(x, y) ∈ ∂Ω,
(1d)
f (0, 0) = f (1, 0) = f (0, 1) = f (1, 1) = 0, ∂ i+j f ∂ i+j f (0, 0) = (0, 1) = 0, 0 i + j 5. ∂x i ∂y j ∂x i ∂y j
(1e) (1f)
The solution of (1) is decomposed [11] into a sum of components u = v + wR + wT + wB + wT R + wBR where v is the regular component, wR is a regular boundary layer function associated with the right edge x = 1, wT , wB are parabolic layer functions associated with the top y = 1 and bottom y = 0 sides of the square and wT R , wBR are corner layer functions near the outflow corners (0, 1), (1, 1). The decomposition isolates each layer component as a solution of the homogeneous differential equation with non-zero boundary conditions on only one side or near one single corner of the domain. The allure of computational approaches to solving partial differential equations is that they can be applied to problems for which there is no known closed form representation of the exact solution available. Of course, although a computational algorithm is widely applicable, the accuracy of the computed approximations is not necessarily guaranteed. As always, there is a wide gap between the complicated systems of non-linear partial differential equations which practical modellers wish to solve and the simpler mathematical problems for which theoretical error bounds exist. The problem Hemker [3] formulated lies somewhere in the middle of these positions. The tools of analysis presented here are applicable to a wide class of singularly perturbed partial differential equations [11]. In the establishment of the asymptotic error bounds, we obtain upper bounds on the components of the solution rather than explicit closed form representations of these components. We hope that the tools of analysis used here may prove useful in designing and analyzing a new appropriate numerical algorithm for the benchmark problem proposed by Hemker [3]. 2. Extensions and compatibility issues Note that throughout this paper C denotes a constant independent of ε. Let D ⊂ R2 be an open set containing the closed unit square Ω = [0, 1] × [0, 1]. Let D1 , D2 be open sets so that D ⊂ D1 ⊂ D2 . For any f ∈ C n (D) there exists an extension f ∗ ∈ C n (D2 ) of f such that f ∗ (x, y) = f (x, y), (x, y) ∈ Ω; f ∗ (x, y) = C, (x, y) ∈ D2 \ D1 . The space C γ (D) is the set of all functions that are Hölder continuous of degree γ with respect to the Euclidean norm · e . That is f ∈ C γ (D) if f 0,γ ,D =
|f (u) − f (v)| γ u − ve u =v, u,v∈D sup
is finite. The space C n,γ (D) is the set of all functions in C n (D) whose derivatives of order n are Hölder continuous of ∂ i+j z γ degree γ . That is, C n,γ (D) = {z: ∂x i ∂y j ∈ C (D), 0 i + j n}. Also · n,γ are the associated norms and · n,γ is the associated Hölder semi-norm defined by ∂ nv , v = |v|k + v n,γ , v n,γ = n,γ ∂x i ∂y j 0,γ i+j =n
0kn
E. O’Riordan, G.I. Shishkin / Applied Numerical Mathematics 58 (2008) 1761–1772
|v|k =
i+j =k
i+j ∂ v sup i j , ∂x ∂y
vk =
1763
|v|j .
0j k
In the decomposition of the solution of (1), some of the components will be solutions of parabolic differential equations. Thus we need to introduce the relevant function spaces. The space C 0+γ (D) is the set of all functions that are Hölder continuous of degree γ with respect to the metric d(u, v) , where for all u, v ∈ R2 d(u, v) = (u1 − v1 )2 + |u2 − v2 |. For f to be in C 0+γ (D) the following semi-norm needs to be finite |f (u) − f (v)| f 0+γ ,D = sup . γ u =v, u,v∈D (d(u, v)) The space C n+γ (D) is defined by ∂ i+j z 0+γ C n+γ (D) = z: ∈ C (D), 0 i + 2j n ∂x i ∂y j and ∂ nv , v = |v|k + v n+γ , v n+γ = n+γ ∂x i ∂y j 0+γ i+j =n
0kn
are the associated semi-norms and norms. From Volkov [12] and Ladyzhenskaya and Ural’tseva [5] (see also [2, Theorem 3.2]) if f ∈ C 1,γ (Ω), f (0, 0) = f (1, 0) = f (0, 1) = f (1, 1) = 0 then the solution of (1) u is in C 3,γ (Ω). Using the stretching transformations η = x/ε, ζ = y/ε, problem (1) transforms to −(u˜ ηη + u˜ ζ ζ ) + a˜ u˜ η = ε f˜, (η, ζ ) ∈ Ωε = (0, 1/ε) × (0, 1/ε) where z˜ (η, ζ ) = z(x, y). Note that the dimension of the stretched domain depend on ε. From Ladyzhenskaya and Ural’tseva [5, p. 110] applied to the problem posed on the stretched domain we have for k = 0, 1
u ˜ 2+k,γ C εf˜k,γ + |u|0 . Transforming back to the original variables, we deduce the following global a priori bounds on the solution of (1). For k = 0, 1 |u|0 C|f |0 ,
|u|1 + ε u 1,γ C |f |0 + ε γ f 0,γ + ε −1 |u|0 , k γ −1+i−k −1+γ −(k+2) |u|2+k + ε u 2+k,γ C ε |f |i + ε f k,γ + ε |u|0 . γ
(2a) (2b) (2c)
i=0
The derivative estimates in (2) can be localized by using the argument in Ladyzhenskaya and Ural’tseva [5, pp. 132– 134]. Note that this argument is applicable to points on the boundary. Hence, given any point (η0 , ζ0 ) ∈ Ω ε and any δ˜ = O(1) > 0, let
N˜ ˜ (η0 , ζ0 ) = B (η0 , ζ0 ); δ˜ ∩ Ω ε δ
where B(a; R) = {x ∈ R2 such that x − ae < R}. Then
u ˜ C εf˜ ˜ + |u| ˜ , k = 0, 1. ˜ 2+k,Nδ˜ ,γ
k,N2δ˜ ,γ
0,N2δ˜
Hence, given any point x0 ∈ Ω and any δ = O(ε) > 0, let Nδ (x0 ) = B(x0 ; δ) ∩ Ω. Then we have the following bounds on the solution of problem (1). For k = 0, 1 and u = 0, (x, y) ∈ ∂Ω, we have
(3a) |u|1,Nδ + ε γ u 1,Nδ ,γ C |f |0,N2δ + ε γ f 0,N2δ ,γ + ε −1 |u|0,N2δ , k |u|2+k,Nδ + ε γ u 2+k,Nδ ,γ Cε −1 ε i−k |f |i,N2δ + ε γ f k,N2δ ,γ + ε −(k+1) |u|0,N2δ . (3b) i=0
1764
E. O’Riordan, G.I. Shishkin / Applied Numerical Mathematics 58 (2008) 1761–1772
Note that if u = g = 0, (x, y) ∈ ∂Ω the following additional boundary terms (see [7, Theorem 3.2]) 2 −1 i 2+γ +ε ε |g|i,N2δ ∩∂Ω + ε g 2,γ ,N2δ ∩∂Ω , i=0
+ε
−(2+k)
2+k
(4a)
ε |g|i,N2δ ∩∂Ω + ε i
k+γ
g 2+k,γ ,N2δ ∩∂Ω
(4b)
i=0
are added to the right-hand side of (3a), (3b), respectively. 3. Bounds on the regular component Define the extended domain Ω ∗ = (0, 1 + d) × (−d, 1 + d), d > 0. The reduced solution v0∗ is the solution of the first order problem a∗
∂v0∗ = f ∗, ∂x
x > 0, (x, y) ∈ Ω ∗ ;
v0∗ (0, y) = 0.
The extension f ∗ is constructed to be zero in the neighbourhood of the four corners of Ω ∗ , which implies that v0∗ is constant in a neighbourhood of the four corners of Ω ∗ . The first order correction to the reduced solution is the solution of ∂v1∗ v1∗ (0, y) = 0 = − v0∗ , x > 0, (x, y) ∈ Ω ∗ ; ∂x and the second order correction is the solution of the elliptic problem a∗
−ε v2∗ + a ∗ (v2∗ )x = − v1∗ ,
(x, y) ∈ Ω ∗ ;
v2∗ = 0,
(x, y) ∈ ∂Ω ∗ .
The extensions are such that v2∗ ∈ C 3,γ (Ω ∗ ). Let v ∗ = v0∗ + εv1∗ + ε 2 v2∗ and define the regular component v to be the restriction of v ∗ to the closed domain Ω. Hence it is the solution of −ε v + avx = f,
(x, y) ∈ Ω,
v = v∗,
(x, y) ∈ ∂Ω.
(5)
Note that v0 and v1 do not depend on ε and using the bounds (2) we deduce the following result. Lemma 1. The regular component v defined in (5) satisfies the following bounds on its derivatives |v|0 C
and |v|k + ε γ v k,γ C(1 + ε 2−k ),
k = 1, 2, 3.
(6)
Note further that if f ∈ C k+2,γ (Ω) and ∂ i+j f (0, 1) = 0, ∂x i ∂y j
∀i + j k + 2
then for the first two components of v, for m = 0, 1, ∂ i+j vm (0, 1) = 0, i + j k ∂x i ∂y j
and
∂ k+1 vm (0, 1) = 0. ∂x k+1
(7)
4. Bounds on the regular layer component For the regular layer component wR , define the extended domain Ω ∗∗ = (0, 1) × (−d, 1 + d), d > 0. The extended regular layer component is defined to be ∗ = 0, L∗ w R
(x, y) ∈ Ω ∗∗
∗ wR (1, y) = (u − v)∗ (1, y),
∗ ∗ ∗ wR (0, y) = wR (x, −d) = wR (x, 1 + d) = 0
E. O’Riordan, G.I. Shishkin / Applied Numerical Mathematics 58 (2008) 1761–1772
1765
∗ ∈ C 3,γ (Ω ∗∗ ) and (u − v)∗ (1, y) is defined in such a way that (u − v)∗ (1, y) = 0, y < −d/2 or y > 1 + d/2. Thus wR and using a maximum principle we deduce that ∗ α(1−x) w (x, y) Ce− ε , (x, y) ∈ Ω ∗∗ . (8) R ∗ Consider the following expansion of wR ∗ wR (x, y) = (u − v)∗ (1, y)φ(x, y) + εzR (x, y)
(9)
where for all y ∈ [−d, 1] the function φ is the solution of the boundary value problem −εφxx + a ∗ (1, y)φx = 0, t n e−t
x ∈ (0, 1),
φ(0, y) = 0,
φ(1, y) = 1.
Ce−t/2 , n 1, t
0, we have Note that, by using ∂φ a ∗ (1,y) a ∗ (1,y)(1−x) a ∗ (1,y)(1−x) C 1 e− ε + 1 − x e− ε 2ε Ce− , ∂y ε ε 2 ∗ ∗ ∂ φ ∂φ C a∗ (1,y)(1−x) − a (1,y)(1−x) − a (1,y)(1−x) ε 2ε ε |φ| e , 2 Ce , e− . ∂x ε ∂y Note that zR = 0 on ∂Ω ∗∗ and ∗ α(1−x) α(1−x) L zR (x, y) C 1 + 1 − x e− 2ε C e− 4ε . ε ε ε α(1−x)
One can deduce that |zR (x, y)| Ce− 4ε and then using the bounds (2) and the fact that ε γ e−α(1−x)/ε 0,γ ,[0,1] C, we have |zR |k + ε γ zR k,γ Cε −k , k = 1, 2, 3. Using the local bounds given in (4), we have that for i = 1, 2, 3 and j = 2, 3 i ∗ j ∗ ∂ wR ∂ wR −i −α(1−x)/4ε 1−j −α(1−x)/4ε , . (10) ∂x i (x, y) Cε e ∂y j (x, y) Cε e ∗ to the original domain and is the solution of The regular layer component wR is defined to be the restriction of wR
LwR = 0,
(x, y) ∈ Ω,
(11a)
wR (1, y) = (u − v)(1, y), ∗ (x, 0), wR (x, 0) = wR
Given the above bounds on
∗ wR
wR (0, y) = 0,
(11b)
∗ wR (x, 1) = wR (x, 1).
(11c)
and its derivatives, we have the following result.
Lemma 2. The regular layer component wR , which is the solution of (11), satisfies the following bounds at all points (x, y) ∈ Ω i ∂ wR α(1−x) Cε −i e−α(1−x)/4ε , i = 1, 2, 3, wR (x, y) Ce− ε , (x, y) ∂x i j ∂ wR 1−j −α(1−x)/4ε , j = 2, 3. ∂y j (x, y) Cε e 5. Bounds on the parabolic layer component Define the extended domain Ω ∗∗∗ = (0, 1 + d) × (0, 1), d > 0. The extended top layer component is defined to be L∗ wT∗ = 0,
(x, y) ∈ Ω ∗∗∗ ,
wT∗ (x, 1) = (u − v)∗ (x, 1), wT∗ (1 + d, y)
wT∗ (0, y) = wT∗ (x, 0) = 0,
is specified below.
The extension (u − v)∗ (x, 1) is defined in such a way that (u − v)∗ (x, 1) = 0, x > 1 + d/2. Note that L(u − v) = 0 and thus u − v is sufficiently compatible with 0 at (0, 0) and (0, 1) so that the extension can be arranged for wT∗ ∈
1766
E. O’Riordan, G.I. Shishkin / Applied Numerical Mathematics 58 (2008) 1761–1772
C 3,γ (Ω¯ ∗∗∗ ). The top parabolic layer component wT is defined to be the restriction of wT∗ to the original domain and is the solution of LwT = 0 and wT (x, 1) = (u − v)(x, 1), wT (x, 0) = 0, wT (0, y) = 0, wT (1, y) = wT∗ (1, y). Consider the following parabolic problems (where x is a time-like variable) L∗p w0∗ ≡ −ε(w0∗ )yy + a ∗ (x, y)(w0∗ )x = 0, w0∗ (0, y) = w0∗ (x, 0) = 0,
(x, y) ∈ (0, 1 + d) × (0, 1),
w0∗ (x, 1) = −(v0∗
(12a)
+ εv1∗ )(x, 1),
(12b)
and L∗p w1∗ = (w0∗ )xx ,
w1∗ (0, y) = 0 = w1∗ (x, 0) = w1∗ (x, 1).
(12c)
Let w2∗ be the solution of the elliptic problem L∗ w2∗ = (w1∗ )xx ,
w2∗ (0, y) = w2∗ (x, 0) = w2∗ (1 + d, 0) = 0, w2∗ (x, 1) = −(v2∗ )(x, 1).
(12d)
If the boundary value for wT∗ (1 + d, y) is chosen to be wT∗ (1 + d, y) = w0∗ (1 + d, y) + εw1∗ (1 + d, y) then wT = w0 + εw1 + ε 2 w2 ,
(x, y) ∈ Ω.
(12e)
We need to address the issue of sufficient compatibility and regularity for patibility conditions for a parabolic problem of the form (εzxx − cz − f1 ) = bzt From [6],
¯ z ∈ C 2+γ (G)
on G = (0, 1) × (0, T ],
if c, b, f1
¯ ∈ C 0+γ (G)
w2∗
∈
C 3,γ (Ω ∗∗∗ ).
z(x, 0) = z(0, t) = z(1, t) = 0.
First we state com(13)
and
f1 (0, 0) = f1 (1, 0) = 0.
(14a)
Differentiate the differential equation (13) with respect to time to get an equation for zt . From this we can deduce that ¯ if c, b, f1 ∈ C 2+γ (G) ¯ and in addition to (14a) we have z ∈ C 4+γ (G) −ε(f1 /b)xx (0, 0) = (f1 )t (0, 0),
−ε(f1 /b)xx (1, 0) = (f1 )t (1, 0).
(14b)
Differentiate the differential equation (13) twice with respect to time yields an equation for ztt . From this we deduce ¯ if c, b, f ∈ C 4+γ (G) ¯ and in addition to (14a) and (14b) we have that z ∈ C 6+γ (G)
ε ztt (x, 0) xx (0) = (f1 )tt (0, 0), ε ztt (x, 0) xx (1) = (f1 )tt (1, 0) (14c) where
f1 ε f1 c ztt (x, 0) = − − + f1 (x, 0). b t b b xx b2
Consider the particular parabolic problem εyxx − cy − byt = 0, y(x, 0) = 0,
(x, t) ∈ G = (0, 1) × (0, T ], b > 0,
y(0, t) = g0 (t),
(15a)
y(1, t) = 0.
(15b) bg0 (t)(1
Let y = z + (1 − x)g0 (t) where εzxx − cz − bzx = c(1 − x)g0 (t) + − x). In this case, the compatibility ¯ are satisfied if g0 (0) = g (0) = g (0) = 0. From above, we see that if in addition we conditions for y ∈ C 4+γ (G) 0 0 ¯ and if g (i) (0) = 0, i k then y ∈ C 2k+γ (G). ¯ have g0 (0) = 0 then y ∈ C 6+γ (G)) 0 Lemma 3. Assuming the following compatibility ∂ i+j f (0, 1) = 0, ∂x i ∂y j
0 i + j 5,
then w2∗ ∈ C 3,γ (Ω ∗∗∗ ), γ ∈ (0, 0.5), where w2∗ is the solution of (12e). Proof. Let w2∗ = z2∗ − yv2∗ (x, 1). Then z2∗ = 0 on the boundary ∂Ω ∗∗∗ and
L∗ z2∗ = L∗ w2∗ + yL∗ v2∗ (x, 1) = (w1∗ )xx + yL∗ v2∗ (x, 1) , (x, y) ∈ Ω ∗∗∗ .
(16)
E. O’Riordan, G.I. Shishkin / Applied Numerical Mathematics 58 (2008) 1761–1772
1767
Note that (w1∗ )xx (0, 1) = 0 and w2∗ is the solution of an elliptic problem. So to have sufficient compatibility at the corner (0, 1) we also require that L(v2∗ (x, 1))(0, 1) = 0. This follows from the fact that v2∗ ∈ C 3,γ (Ω) and v2∗ (0, y) = 0, 0 y 1. For sufficient regularity of the data we hence require that (w1∗ )xx ∈ C 1,γ (Ω ∗∗∗ ), γ ∈ (0, 0.5). Since w1∗ is a solution of a parabolic equation (where x is the time-like variable) we require that ∂ i+j w1∗ ∈ C γ (Ω ∗∗∗ ), ∂x i ∂y j
0 2i + j 6.
This places constraints on w0∗ . From (7), if we assume (16) then ∂ i (v0∗ + εv1∗ ) (0, 1) = 0, ∂x i
i 4.
On ΓT∗ , −w0∗ = v0∗ + εv1∗ so from the discussion on problem (15), the assumption (16) suffices for w0∗ ∈ C 8+2γ (Ω ∗∗∗ ). Recall that L∗p w1∗ = (w0∗ )xx ≡ f3 . By virtue of (16) ∂ i+2 (v0∗ + εv1∗ ) ∂ i f3 (0, 1) = (0, 1) = 0, i 2. ∂x i ∂x i+2 Note also that ε(f3 )yy = (a ∗ (w0∗ )x )xx which implies that
(17a)
∂ 2 f3 (0, 1) = 0. ∂y 2
(17b)
Note further that since w0 (0, y) = 0 ∂ ε ∂ 2 w0 ∂ 2 w0 (0, 1) = 0 (0, 1) = ∂y∂x ∂y a ∗ ∂y 2 and since w0∗ ∈ C 8+2γ (Ω ∗∗∗ ) we also have that ∂ 3 w0 ∂ 3 w0 ∂ 4 w0 (0, 1) = (0, 1) = (0, 1) = 0. ∂y 2 ∂x ∂y∂x 2 ∂y∂x 3 Combining these with ε 2 (f3 )yyyy = ε(a ∗ (w0 )x )xxyy implies that ∂f3 (0, 1) = 0, ∂y Finally ε
∂ 4 f3 ∂ 2 f3 (0, 1) = 0. (0, 1) = 4 ∂x∂y ∂y
∂ 5 w0 ∂3 ∗ ∂w0 , a = ∂x ∂y 5 ∂y 3
ε
∂ 7 w0 ∂5 ∗ ∂w0 , a = ∂y 7 ∂x ∂y 5
(17c)
ε
∂ 5 w0 ∂3 ∗ ∂w0 a = ∂x ∂x 3 ∂y 2 ∂x 3
will result in the following be satisfied ∂ 3 f3 ∂ 3 f3 (0, 1) = (0, 1) = 0. 2 ∂x∂y ∂y 3
(17d)
Together the identities in (17) yield sufficient compatibility (see (14)) for ∂ i+j w1∗ ∈ C 0+2γ (Ω ∗∗∗ ), ∂x i ∂y j
0 2i + j 6.
Note that for f sufficiently regular, f 0,γ Cf 0 + f 0+2γ , and thus ∂ i+j w1∗ ∈ C γ (Ω ∗∗∗ ), ∂x i ∂y j
0 2i + j 6.
2
Using the above expansion for wT we have the following bounds on the parabolic layer function associated with the edge y = 1.
1768
E. O’Riordan, G.I. Shishkin / Applied Numerical Mathematics 58 (2008) 1761–1772
Lemma 4. For the solution wT of (12) and i, j = 2, 3 √ wT (x, y) Ce−(1−y)/ ε , j i √ √ ∂ wT ∂ wT Cε −j/2 e−(1−y)/ ε , Cε 2−i e−(1−y)/ ε . ∂y j ∂x i
(18a) (18b)
Proof. By choosing the appropriate barrier function √ ∗ w (x, y) Ce αx e−(1−y)/ ε . 0
Note also that w0∗ = −(v0∗ + εv1∗ ) on ΓT∗ . Assuming sufficient regularity on a and f then we have that i ∗ ∂ (v0 + εv1∗ ) C, 0 i 4. ∂x i ∗ √ and noting that L∗ Introducing the stretched variable ζ = 1−y p w0 = 0 we deduce from the corresponding bounds to ε (3) for the parabolic problem (using the argument in [8, Theorem 4]) i+j ∂ w˜ 0 (x, ζ ) Ce−ζ , 0 2i + j 8. ∂x i ∂ζ j 0+2γ
This implies that i+j 2 ∂ ∂ w˜ 0 Ce−ζ , ∂x i ∂ζ j ∂x 2 0+2γ Noting that L∗p w1∗ = (w0∗ )xx , (x, y) ∈ G, x α
√ e−(1−y)/ ε .
0 2i + j 4. w1 = 0, (x, y) ∈ S and by choosing appropriate barrier function we have
that |w1 (x, y)| Ce Again using the argument from [8] we have for 0 2i + j 6 i+j ∂ w˜ 1 (x, ζ ) Ce−ζ , 0 2i + j 6, ∂x i ∂ζ j 0+2γ which implies that i+j 2 ∂ ∂ w˜ 1 Ce−ζ , ∂x i ∂ζ j ∂x 2 0+2γ
0 2i + j 2.
The above inequalities imply that i+j ∗ i+j ∗ √ ∂ w0 ∂ w1 −j/2 −(1−y)/ ε e , ∂x i ∂y j + ∂x i ∂y j Cε
(19)
0 i + j 3. 2x
√
Finally for ε sufficiently small (4ε < α 2 ) |w2∗ (x, y)| Ce α e−(1−y)/ ε . Note that from (19), we have that for γ ∈ (0, 0.5) and k = 0, 1, k = i + j i+j i+j (1−y) ∂ ∂ γ −j/2 − √ε ˜ ˜ + ε Cε e , 0 2i + j 2, f f ∂x i ∂y j ∂x i ∂y j k k+2γ where f˜ = (w1∗ )xx . From the bounds (2) and the derivative bounds on v2 , we deduce that
|w2∗ |2 C ε −1 |f˜|0 + ε −1+γ f˜ 0,γ + ε −2 |w2∗ |0 + Cε −2 Cε −2 ,
|w2∗ |3 C ε −2 |f˜|0 + ε −1 |f˜|1 + ε −1+γ f˜ 1,γ + ε −3 |w2∗ |0 + Cε −3 Cε −3 since for any f ∈ C 0+2γ , f 0,γ Cf 0 + C f 0+2γ . This completes the proof. The bottom parabolic layer component wB is defined analogously and
2
E. O’Riordan, G.I. Shishkin / Applied Numerical Mathematics 58 (2008) 1761–1772 y wB (x, y) Ce− √ε ,
1769
3 y ∂ wB Cε −3/2 e− √ε , (x, y) ∂y 3
i y ∂ wB Cε 2−i e− √ε , (x, y) ∂x i
i = 2, 3.
6. Corner layer functions Finally we define the outflow corner layer function wT R (and wBR is defined analogously) which is defined to be the solution of LwT R = 0,
(x, y) ∈ Ω = (0, 1) × (0, 1);
wT R (1, y) = −wT (1, y),
wT R (x, 0) = 0,
wT R (0, y) = 0,
wT R (x, 1) = −wR (x, 1).
Note that wT is compatible with u − v at the corner (1, 1) and, in turn, u − v is compatible with wR at the corner (1, 1). Hence wT R ∈ C 3,γ (Ω). From the comparison principle and the bounds on wT and wR , we have that for ε sufficiently small (1−y) α(1−x) wT R (x, y) Ce− 2ε e− √ε .
(20)
We can write this corner layer function as the following sum
wT R (x, y) = wT R (1, y)φ1 (x) + wT R (x, 1) − wT R (1, 1)φ1 (x) φ2 (y) + εzT R (x, y) where −εφ1 (x) + a(1, 1)φ1 (x) = 0,
φ1 (0) = 0,
φ1 (1) = 1,
φ2 (y) =
1 − e−y/
√
ε √ . 1 − e−1/ ε
Then zT R (x, y) = 0 on the boundary of Ω and noting that −ε
∂2 ∂ ∂2 wR (x, 1) + a(x, y) wR (x, 1) = ε 2 wR (x, 1) 2 ∂x ∂x ∂y
and using the bounds on wR corresponding to (10) we get that ε LzT R (x, y) Ce−α(1−x)/(4ε) . As before, we deduce the derivative bounds |zT R |k Cε −k , k = 1, 2, 3. One can then establish the bounds i i ∂ wT R Cε −i , ∂ wT R Cε −i , i = 2, 3. ∂x i ∂y i Also with the above decomposition of wT R we have the additional bound 3 ∂ wT R −2 ∂y 3 Cε .
(21)
(22)
Lemma 5. The corner layer component wT R satisfies the bounds on its derivatives given in (20), (21) and (22). 7. Discrete problem and error analysis We present error bounds for the approximations generated by using a standard upwind difference operator LN U = −ε(δx2 + δy2 )U + aDx− U = f,
(xi , yj ) ∈ Ω N
(23a)
on a mesh Ω N = ωx × ωy , which is a tensor product of two piecewise-uniform one-dimensional meshes ωx , ωy . The finite difference operators Dx− and δx2 are the standard first order backward difference and the second order centered difference on a non-uniform mesh [1]. Here the mesh ωx places N/2 mesh intervals into both [0, σx ] and [σx , 1] and
1770
E. O’Riordan, G.I. Shishkin / Applied Numerical Mathematics 58 (2008) 1761–1772
the mesh ωy places N/4 mesh intervals into the intervals [0, σy ] and [1 − σy , 1] and N/2 mesh intervals into the region [σy , 1 − σy ]. The transition parameters are taken to be √ ε σx = min 0.5, ln N (23b) and σy = min{0.25, ε ln N }. α The discrete solution is decomposed in an analogous fashion to the continuous solution. That is U = V + WR + WT + WBR + WT R where LN V = f, LN W = 0 (xi , yj ) ∈ Ω N and V = v, W = w (xi , yj ) ∈ ∂Ω N . On an arbitrary mesh using the bounds (6) one has the truncation error estimate 3 3 2 N L (v − V ) CN −1 ε ∂ v + ∂ v + CN −1 ε ∂ v CN −1 ∂y 3 ∂x 3 ∂x 2 and hence v − V Ω N CN −1 . If σx = 0.5 then the same argument coupled with the bounds (2) yields u − U Ω N CN −1 (ln N)2 . (If σy = 0.25 and σx < 0.5, then a minor modification to the argument below yields the same error bound.) In the remainder of this section we assume that √ ε σx = ln N and σy = ε ln N. α Note that from (8) and by the choice of transition point σx wR (x, y) CN −1 , x 1 − σx . Consider the discrete one-dimensional barrier function Ψ (xi ) which is the solution of −εδx2 Ψ (xi ) + αD − Ψ (xi ) = 0,
xi ∈ ω x ,
Ψ (0) = 0,
Ψ (1) = 1.
Note also that for the top y = 1 and bottom y = 0 edges, by (8), α(1−x ) WR (xi , yj ) = wR (xi , yj ) Ce− ε i CΨ (xi ) + CN −1 ,
yj = 0, 1.
Using CΨ (xi ) + CN −1 as a barrier function we deduce that |WR (xi , yj )| CΨ (xi ) + CN −1 and WR (xi , yj ) CN −1 , xi 1 − σx . Thus |(WR − wR )(xi , yj )| CN −1 , xi 1 − σx . In the fine mesh region (1 − σx , 1) × (0, 1), using the bounds (10) we have the truncation error bound 2 3 3 −1 N L (wR − WR ) CN −1 σx ε ∂ wR + ∂ wR + CN −1 ε ∂ wR C N ln N . ∂x 3 ∂x 2 ∂y 3 ε Use the barrier function
N −1 ln N C xi − (1 − σx ) + CN −1 ε in the region (1 − σx , 1) × (0, 1) to derive the error bound WR − wR Ω N CN −1 (ln N )2 . Next we examine the top parabolic layer function. Consider the discrete barrier function Υ (xi )Φ(yj ) where Φ(0) = 0, −εδy2 Φ(yj ) + Φ(yj ) = 0, yj ∈ ωy , i 2hj 1+ , hj = xj − xj −1 . Υ (xi ) = α
Φ(1) = 1,
j =1
Assuming the strict inequality 8ε < α 2 and for N sufficiently large (independently of ε) we have that −εδx2 Υ (xi ) + αDx− Υ (xi ) − Υ (xi ) 0,
xi ∈ ω x .
E. O’Riordan, G.I. Shishkin / Applied Numerical Mathematics 58 (2008) 1761–1772
1771
Note also that √ WT (1, yj ) Ce−(1−yj )/ ε CΦ(yj ) + CN −1 . Then |WT (xi , yj )| CΥ (xi )Φ(yj ) + CN −1 and hence |WT (xi , yj )| CN −1 , yj 1 − σy . In the fine mesh region (0, 1) × (1 − σy , 1) use the bounds (18) to obtain 2 3 3 N L (wT − WT ) CN −1 ε ∂ wT + ∂ wT + Cσy N −1 ε ∂ wT CN −1 ln N. ∂y 3 ∂x 3 ∂x 2 Hence wT − WT Ω N CN −1 ln N. As for the discrete regular and parabolic layer components we have that WT R (xi , yj ) CΨ (xi ) + CN −1 CN −1 , xi 1 − σx , WT R (xi , yj ) CΥ (xi )Φ(yj ) + CN −1 CN −1 , yj 1 − σy . In the fine mesh corner region (1 − σx , 1) × (1 − σy , 1), use the bounds given in (21) and (22) to obtain −1 N L (wT R − WT R ) C N ln N . ε
Use the barrier function C(εN )−1 ln N (xi − (1 − σx )) + CN −1 in the corner region (1 − σx , 1) × (1 − σy , 1) to derive the error bound WT R − wT R Ω N CN −1 (ln N )2 . Theorem 6. If U is the solution of (23) and u is the solution of the elliptic problem (1) then for 8ε < α 2 and N sufficiently large (independently of ε) U − uΩ N CN −1 (ln N )2 . Numerical results in [1] agree with this theoretical bound on the nodal error. Furthermore, numerical experiments in [1] demonstrate that √ the numerical method (23) also produces parameter-uniform approximations to the scaled derivatives εux and εuy . 8. Conclusions A proof of parameter-uniform pointwise convergence of a set of numerical approximations to a class of singularly perturbed elliptic equations is given. The tools of analysis used involve a special decomposition of the solution into a sum of a regular component and several types of layer components. In particular, the solutions contain parabolic boundary layers. Parameter-explicit pointwise bounds on the derivatives of all of these components are established. The establishment of these bounds on the derivatives is the main component in the proof of first order (up to logarithmic factors) parameter-uniform convergence of the numerical approximations. References [1] P.A. Farrell, A.F. Hegarty, J.J.H. Miller, E. O’Riordan, G.I. Shishkin, Robust Computational Techniques for Boundary Layers, Chapman and Hall/CRC Press, Boca Raton, FL, 2000. [2] H. Han, R.B. Kellogg, Differentiability properties of solutions of the equation −ε 2 u + ru = f (x, y) in a square, SIAM J. Math. Anal. 21 (1990) 394–408. [3] P.W. Hemker, A singularly perturbed model problem for numerical computation, J. Comput. Appl. Math. 76 (1996) 277–285. [4] R.B. Kellogg, M. Stynes, Corner singularities and boundary layers in a simple convection diffusion problem, J. Differential Equations 213 (2005) 81–120. [5] O.A. Ladyzhenskaya, N.N. Ural’tseva, Linear and Quasilinear Elliptic Equations, Academic Press, New York and London, 1968. [6] O.A. Ladyzhenskaya, V.A. Solonnikov, N.N. Ural’tseva, Linear and Quasilinear Equations of Parabolic Type, Transactions of Mathematical Monographs, vol. 23, American Mathematical Society, Providence, RI, 1968. [7] T. Linß, M. Stynes, Asymptotic analysis and Shishkin-type decomposition for an elliptic convection–diffusion problem, J. Math. Anal. Appl. 261 (2001) 604–632. [8] J.J.H. Miller, E. O’Riordan, G.I. Shishkin, L.P. Shishkina, Fitted mesh methods for problems with parabolic boundary layers, Math. Proc. Roy. Irish Acad. A 98 (2) (1998) 173–190.
1772
E. O’Riordan, G.I. Shishkin / Applied Numerical Mathematics 58 (2008) 1761–1772
[9] H.-G. Roos, Optimal convergence of basic schemes for elliptic boundary value problems with strong parabolic layers, J. Math. Anal. Appl. 267 (2002) 194–208. [10] S. Shih, R.B. Kellogg, Asymptotic analysis of a singular perturbation problem, SIAM J. Math. Anal. 18 (1987) 1467–1511. [11] G.I. Shishkin, Discrete Approximation of Singularly Perturbed Elliptic and Parabolic Equations, Russian Academy of Sciences, Ural Section, Ekaterinburg, 1992 (in Russian). [12] E.A. Volkov, Differentiability properties of solutions of boundary value problems for the Laplace and Poisson equations, Proc. Steklov Inst. Math. 77 (1965) 101–126.
Parameter uniform numerical methods for singularly perturbed elliptic problems with parabolic boundary layers E. O’Riordan a,∗ , G.I. Shishkin b,1 a School of Mathematical Sciences, Dublin City University, Ireland b Institute of Mathematics and Mechanics, Russian Academy of Sciences, Ekaterinburg, Russia
Available online 22 November 2007
Abstract Parameter-uniform numerical methods for a singularly perturbed elliptic problem with parabolic boundary layers in the solution are analyzed. A priori parameter explicit bounds on the solution and its derivatives are obtained using a suitable decomposition of the solution into regular and layer components. Based on this decomposition a numerical method is constructed and analyzed. © 2007 IMACS. Published by Elsevier B.V. All rights reserved. MSC: 65N15 Keywords: Parabolic boundary layers; Elliptic equation; Piecewise-uniform mesh
1. Introduction Hemker [3] formulated a singular perturbation test problem, which can be viewed as a simple mathematical model of flow past a circular obstacle. This problem also acts as a mathematical challenge to numerical analysts interested in designing parameter-uniform numerical methods for singularly perturbed differential equations. In this paper, we examine a related simpler problem, where the circular object is absent and we confine the discussion to problems where no interior parabolic layers are present in the solution. In this paper, we will examine numerical methods for the problem Lu ≡ −εu + aux = f, u = g,
(x, y) ∈ Ω = (0, 1)2 ,
(x, y) ∈ ∂Ω,
a > α > 0,
∞
a ∈ C (D),
(1a) (1b)
f ∈C
5,γ
(D),
Ω ⊂ D.
(1c)
Asymptotic expansions of the constant coefficient version of this problem [10] indicate the complicated analytical nature of the solution, which will in general contain a regular exponential boundary layer, two parabolic (or character* Corresponding author.
E-mail addresses: [email protected] (E. O’Riordan), [email protected] (G.I. Shishkin). 1 The research of the second author was supported by the Boole Centre for Research in Informatics at the National University of Ireland in
Cork, by the Mathematics Applications Consortium for Science and Industry in Ireland (MACSI) under the Science Foundation Ireland (SFI) Mathematics Initiative and by the Russian Foundation for Basic Research under grant No. 07-01-00729. 0168-9274/$30.00 © 2007 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2007.11.003
1762
E. O’Riordan, G.I. Shishkin / Applied Numerical Mathematics 58 (2008) 1761–1772
istic) boundary layers and several types of corner layer functions in the vicinity of the inflow and outflow corners. To analyze the convergence of any set of numerical approximations, we require bounds on the solution and its derivatives. In the case of constant coefficients, Kellogg and Stynes [4] supply these parameter explicit bounds on the derivatives of the solution in the general case of g ∈ / C 0 (∂Ω). Again, these bounds on the derivatives display the non-trivial nature of the solution even in the case of constant coefficients. Roos [9] also examines a constant coefficient version of problem (1), where a decomposition of the solution into a regular component and several layer components is given. In [9] bounds on the derivatives of the layer components are established by first deriving derivative bounds on the boundaries and then extending these throughout the domain by simply differentiating the differential equation. It is not easy to extend this argument to variable coefficient problems. In this paper, we examine the case of variable coefficients. However, here we impose sufficient regularity and compatibility on the data so that the amplitude of the inflow corner layer functions are negligible. u = g = 0,
(x, y) ∈ ∂Ω,
(1d)
f (0, 0) = f (1, 0) = f (0, 1) = f (1, 1) = 0, ∂ i+j f ∂ i+j f (0, 0) = (0, 1) = 0, 0 i + j 5. ∂x i ∂y j ∂x i ∂y j
(1e) (1f)
The solution of (1) is decomposed [11] into a sum of components u = v + wR + wT + wB + wT R + wBR where v is the regular component, wR is a regular boundary layer function associated with the right edge x = 1, wT , wB are parabolic layer functions associated with the top y = 1 and bottom y = 0 sides of the square and wT R , wBR are corner layer functions near the outflow corners (0, 1), (1, 1). The decomposition isolates each layer component as a solution of the homogeneous differential equation with non-zero boundary conditions on only one side or near one single corner of the domain. The allure of computational approaches to solving partial differential equations is that they can be applied to problems for which there is no known closed form representation of the exact solution available. Of course, although a computational algorithm is widely applicable, the accuracy of the computed approximations is not necessarily guaranteed. As always, there is a wide gap between the complicated systems of non-linear partial differential equations which practical modellers wish to solve and the simpler mathematical problems for which theoretical error bounds exist. The problem Hemker [3] formulated lies somewhere in the middle of these positions. The tools of analysis presented here are applicable to a wide class of singularly perturbed partial differential equations [11]. In the establishment of the asymptotic error bounds, we obtain upper bounds on the components of the solution rather than explicit closed form representations of these components. We hope that the tools of analysis used here may prove useful in designing and analyzing a new appropriate numerical algorithm for the benchmark problem proposed by Hemker [3]. 2. Extensions and compatibility issues Note that throughout this paper C denotes a constant independent of ε. Let D ⊂ R2 be an open set containing the closed unit square Ω = [0, 1] × [0, 1]. Let D1 , D2 be open sets so that D ⊂ D1 ⊂ D2 . For any f ∈ C n (D) there exists an extension f ∗ ∈ C n (D2 ) of f such that f ∗ (x, y) = f (x, y), (x, y) ∈ Ω; f ∗ (x, y) = C, (x, y) ∈ D2 \ D1 . The space C γ (D) is the set of all functions that are Hölder continuous of degree γ with respect to the Euclidean norm · e . That is f ∈ C γ (D) if f 0,γ ,D =
|f (u) − f (v)| γ u − ve u =v, u,v∈D sup
is finite. The space C n,γ (D) is the set of all functions in C n (D) whose derivatives of order n are Hölder continuous of ∂ i+j z γ degree γ . That is, C n,γ (D) = {z: ∂x i ∂y j ∈ C (D), 0 i + j n}. Also · n,γ are the associated norms and · n,γ is the associated Hölder semi-norm defined by ∂ nv , v = |v|k + v n,γ , v n,γ = n,γ ∂x i ∂y j 0,γ i+j =n
0kn
E. O’Riordan, G.I. Shishkin / Applied Numerical Mathematics 58 (2008) 1761–1772
|v|k =
i+j =k
i+j ∂ v sup i j , ∂x ∂y
vk =
1763
|v|j .
0j k
In the decomposition of the solution of (1), some of the components will be solutions of parabolic differential equations. Thus we need to introduce the relevant function spaces. The space C 0+γ (D) is the set of all functions that are Hölder continuous of degree γ with respect to the metric d(u, v) , where for all u, v ∈ R2 d(u, v) = (u1 − v1 )2 + |u2 − v2 |. For f to be in C 0+γ (D) the following semi-norm needs to be finite |f (u) − f (v)| f 0+γ ,D = sup . γ u =v, u,v∈D (d(u, v)) The space C n+γ (D) is defined by ∂ i+j z 0+γ C n+γ (D) = z: ∈ C (D), 0 i + 2j n ∂x i ∂y j and ∂ nv , v = |v|k + v n+γ , v n+γ = n+γ ∂x i ∂y j 0+γ i+j =n
0kn
are the associated semi-norms and norms. From Volkov [12] and Ladyzhenskaya and Ural’tseva [5] (see also [2, Theorem 3.2]) if f ∈ C 1,γ (Ω), f (0, 0) = f (1, 0) = f (0, 1) = f (1, 1) = 0 then the solution of (1) u is in C 3,γ (Ω). Using the stretching transformations η = x/ε, ζ = y/ε, problem (1) transforms to −(u˜ ηη + u˜ ζ ζ ) + a˜ u˜ η = ε f˜, (η, ζ ) ∈ Ωε = (0, 1/ε) × (0, 1/ε) where z˜ (η, ζ ) = z(x, y). Note that the dimension of the stretched domain depend on ε. From Ladyzhenskaya and Ural’tseva [5, p. 110] applied to the problem posed on the stretched domain we have for k = 0, 1
u ˜ 2+k,γ C εf˜k,γ + |u|0 . Transforming back to the original variables, we deduce the following global a priori bounds on the solution of (1). For k = 0, 1 |u|0 C|f |0 ,
|u|1 + ε u 1,γ C |f |0 + ε γ f 0,γ + ε −1 |u|0 , k γ −1+i−k −1+γ −(k+2) |u|2+k + ε u 2+k,γ C ε |f |i + ε f k,γ + ε |u|0 . γ
(2a) (2b) (2c)
i=0
The derivative estimates in (2) can be localized by using the argument in Ladyzhenskaya and Ural’tseva [5, pp. 132– 134]. Note that this argument is applicable to points on the boundary. Hence, given any point (η0 , ζ0 ) ∈ Ω ε and any δ˜ = O(1) > 0, let
N˜ ˜ (η0 , ζ0 ) = B (η0 , ζ0 ); δ˜ ∩ Ω ε δ
where B(a; R) = {x ∈ R2 such that x − ae < R}. Then
u ˜ C εf˜ ˜ + |u| ˜ , k = 0, 1. ˜ 2+k,Nδ˜ ,γ
k,N2δ˜ ,γ
0,N2δ˜
Hence, given any point x0 ∈ Ω and any δ = O(ε) > 0, let Nδ (x0 ) = B(x0 ; δ) ∩ Ω. Then we have the following bounds on the solution of problem (1). For k = 0, 1 and u = 0, (x, y) ∈ ∂Ω, we have
(3a) |u|1,Nδ + ε γ u 1,Nδ ,γ C |f |0,N2δ + ε γ f 0,N2δ ,γ + ε −1 |u|0,N2δ , k |u|2+k,Nδ + ε γ u 2+k,Nδ ,γ Cε −1 ε i−k |f |i,N2δ + ε γ f k,N2δ ,γ + ε −(k+1) |u|0,N2δ . (3b) i=0
1764
E. O’Riordan, G.I. Shishkin / Applied Numerical Mathematics 58 (2008) 1761–1772
Note that if u = g = 0, (x, y) ∈ ∂Ω the following additional boundary terms (see [7, Theorem 3.2]) 2 −1 i 2+γ +ε ε |g|i,N2δ ∩∂Ω + ε g 2,γ ,N2δ ∩∂Ω , i=0
+ε
−(2+k)
2+k
(4a)
ε |g|i,N2δ ∩∂Ω + ε i
k+γ
g 2+k,γ ,N2δ ∩∂Ω
(4b)
i=0
are added to the right-hand side of (3a), (3b), respectively. 3. Bounds on the regular component Define the extended domain Ω ∗ = (0, 1 + d) × (−d, 1 + d), d > 0. The reduced solution v0∗ is the solution of the first order problem a∗
∂v0∗ = f ∗, ∂x
x > 0, (x, y) ∈ Ω ∗ ;
v0∗ (0, y) = 0.
The extension f ∗ is constructed to be zero in the neighbourhood of the four corners of Ω ∗ , which implies that v0∗ is constant in a neighbourhood of the four corners of Ω ∗ . The first order correction to the reduced solution is the solution of ∂v1∗ v1∗ (0, y) = 0 = − v0∗ , x > 0, (x, y) ∈ Ω ∗ ; ∂x and the second order correction is the solution of the elliptic problem a∗
−ε v2∗ + a ∗ (v2∗ )x = − v1∗ ,
(x, y) ∈ Ω ∗ ;
v2∗ = 0,
(x, y) ∈ ∂Ω ∗ .
The extensions are such that v2∗ ∈ C 3,γ (Ω ∗ ). Let v ∗ = v0∗ + εv1∗ + ε 2 v2∗ and define the regular component v to be the restriction of v ∗ to the closed domain Ω. Hence it is the solution of −ε v + avx = f,
(x, y) ∈ Ω,
v = v∗,
(x, y) ∈ ∂Ω.
(5)
Note that v0 and v1 do not depend on ε and using the bounds (2) we deduce the following result. Lemma 1. The regular component v defined in (5) satisfies the following bounds on its derivatives |v|0 C
and |v|k + ε γ v k,γ C(1 + ε 2−k ),
k = 1, 2, 3.
(6)
Note further that if f ∈ C k+2,γ (Ω) and ∂ i+j f (0, 1) = 0, ∂x i ∂y j
∀i + j k + 2
then for the first two components of v, for m = 0, 1, ∂ i+j vm (0, 1) = 0, i + j k ∂x i ∂y j
and
∂ k+1 vm (0, 1) = 0. ∂x k+1
(7)
4. Bounds on the regular layer component For the regular layer component wR , define the extended domain Ω ∗∗ = (0, 1) × (−d, 1 + d), d > 0. The extended regular layer component is defined to be ∗ = 0, L∗ w R
(x, y) ∈ Ω ∗∗
∗ wR (1, y) = (u − v)∗ (1, y),
∗ ∗ ∗ wR (0, y) = wR (x, −d) = wR (x, 1 + d) = 0
E. O’Riordan, G.I. Shishkin / Applied Numerical Mathematics 58 (2008) 1761–1772
1765
∗ ∈ C 3,γ (Ω ∗∗ ) and (u − v)∗ (1, y) is defined in such a way that (u − v)∗ (1, y) = 0, y < −d/2 or y > 1 + d/2. Thus wR and using a maximum principle we deduce that ∗ α(1−x) w (x, y) Ce− ε , (x, y) ∈ Ω ∗∗ . (8) R ∗ Consider the following expansion of wR ∗ wR (x, y) = (u − v)∗ (1, y)φ(x, y) + εzR (x, y)
(9)
where for all y ∈ [−d, 1] the function φ is the solution of the boundary value problem −εφxx + a ∗ (1, y)φx = 0, t n e−t
x ∈ (0, 1),
φ(0, y) = 0,
φ(1, y) = 1.
Ce−t/2 , n 1, t
0, we have Note that, by using ∂φ a ∗ (1,y) a ∗ (1,y)(1−x) a ∗ (1,y)(1−x) C 1 e− ε + 1 − x e− ε 2ε Ce− , ∂y ε ε 2 ∗ ∗ ∂ φ ∂φ C a∗ (1,y)(1−x) − a (1,y)(1−x) − a (1,y)(1−x) ε 2ε ε |φ| e , 2 Ce , e− . ∂x ε ∂y Note that zR = 0 on ∂Ω ∗∗ and ∗ α(1−x) α(1−x) L zR (x, y) C 1 + 1 − x e− 2ε C e− 4ε . ε ε ε α(1−x)
One can deduce that |zR (x, y)| Ce− 4ε and then using the bounds (2) and the fact that ε γ e−α(1−x)/ε 0,γ ,[0,1] C, we have |zR |k + ε γ zR k,γ Cε −k , k = 1, 2, 3. Using the local bounds given in (4), we have that for i = 1, 2, 3 and j = 2, 3 i ∗ j ∗ ∂ wR ∂ wR −i −α(1−x)/4ε 1−j −α(1−x)/4ε , . (10) ∂x i (x, y) Cε e ∂y j (x, y) Cε e ∗ to the original domain and is the solution of The regular layer component wR is defined to be the restriction of wR
LwR = 0,
(x, y) ∈ Ω,
(11a)
wR (1, y) = (u − v)(1, y), ∗ (x, 0), wR (x, 0) = wR
Given the above bounds on
∗ wR
wR (0, y) = 0,
(11b)
∗ wR (x, 1) = wR (x, 1).
(11c)
and its derivatives, we have the following result.
Lemma 2. The regular layer component wR , which is the solution of (11), satisfies the following bounds at all points (x, y) ∈ Ω i ∂ wR α(1−x) Cε −i e−α(1−x)/4ε , i = 1, 2, 3, wR (x, y) Ce− ε , (x, y) ∂x i j ∂ wR 1−j −α(1−x)/4ε , j = 2, 3. ∂y j (x, y) Cε e 5. Bounds on the parabolic layer component Define the extended domain Ω ∗∗∗ = (0, 1 + d) × (0, 1), d > 0. The extended top layer component is defined to be L∗ wT∗ = 0,
(x, y) ∈ Ω ∗∗∗ ,
wT∗ (x, 1) = (u − v)∗ (x, 1), wT∗ (1 + d, y)
wT∗ (0, y) = wT∗ (x, 0) = 0,
is specified below.
The extension (u − v)∗ (x, 1) is defined in such a way that (u − v)∗ (x, 1) = 0, x > 1 + d/2. Note that L(u − v) = 0 and thus u − v is sufficiently compatible with 0 at (0, 0) and (0, 1) so that the extension can be arranged for wT∗ ∈
1766
E. O’Riordan, G.I. Shishkin / Applied Numerical Mathematics 58 (2008) 1761–1772
C 3,γ (Ω¯ ∗∗∗ ). The top parabolic layer component wT is defined to be the restriction of wT∗ to the original domain and is the solution of LwT = 0 and wT (x, 1) = (u − v)(x, 1), wT (x, 0) = 0, wT (0, y) = 0, wT (1, y) = wT∗ (1, y). Consider the following parabolic problems (where x is a time-like variable) L∗p w0∗ ≡ −ε(w0∗ )yy + a ∗ (x, y)(w0∗ )x = 0, w0∗ (0, y) = w0∗ (x, 0) = 0,
(x, y) ∈ (0, 1 + d) × (0, 1),
w0∗ (x, 1) = −(v0∗
(12a)
+ εv1∗ )(x, 1),
(12b)
and L∗p w1∗ = (w0∗ )xx ,
w1∗ (0, y) = 0 = w1∗ (x, 0) = w1∗ (x, 1).
(12c)
Let w2∗ be the solution of the elliptic problem L∗ w2∗ = (w1∗ )xx ,
w2∗ (0, y) = w2∗ (x, 0) = w2∗ (1 + d, 0) = 0, w2∗ (x, 1) = −(v2∗ )(x, 1).
(12d)
If the boundary value for wT∗ (1 + d, y) is chosen to be wT∗ (1 + d, y) = w0∗ (1 + d, y) + εw1∗ (1 + d, y) then wT = w0 + εw1 + ε 2 w2 ,
(x, y) ∈ Ω.
(12e)
We need to address the issue of sufficient compatibility and regularity for patibility conditions for a parabolic problem of the form (εzxx − cz − f1 ) = bzt From [6],
¯ z ∈ C 2+γ (G)
on G = (0, 1) × (0, T ],
if c, b, f1
¯ ∈ C 0+γ (G)
w2∗
∈
C 3,γ (Ω ∗∗∗ ).
z(x, 0) = z(0, t) = z(1, t) = 0.
First we state com(13)
and
f1 (0, 0) = f1 (1, 0) = 0.
(14a)
Differentiate the differential equation (13) with respect to time to get an equation for zt . From this we can deduce that ¯ if c, b, f1 ∈ C 2+γ (G) ¯ and in addition to (14a) we have z ∈ C 4+γ (G) −ε(f1 /b)xx (0, 0) = (f1 )t (0, 0),
−ε(f1 /b)xx (1, 0) = (f1 )t (1, 0).
(14b)
Differentiate the differential equation (13) twice with respect to time yields an equation for ztt . From this we deduce ¯ if c, b, f ∈ C 4+γ (G) ¯ and in addition to (14a) and (14b) we have that z ∈ C 6+γ (G)
ε ztt (x, 0) xx (0) = (f1 )tt (0, 0), ε ztt (x, 0) xx (1) = (f1 )tt (1, 0) (14c) where
f1 ε f1 c ztt (x, 0) = − − + f1 (x, 0). b t b b xx b2
Consider the particular parabolic problem εyxx − cy − byt = 0, y(x, 0) = 0,
(x, t) ∈ G = (0, 1) × (0, T ], b > 0,
y(0, t) = g0 (t),
(15a)
y(1, t) = 0.
(15b) bg0 (t)(1
Let y = z + (1 − x)g0 (t) where εzxx − cz − bzx = c(1 − x)g0 (t) + − x). In this case, the compatibility ¯ are satisfied if g0 (0) = g (0) = g (0) = 0. From above, we see that if in addition we conditions for y ∈ C 4+γ (G) 0 0 ¯ and if g (i) (0) = 0, i k then y ∈ C 2k+γ (G). ¯ have g0 (0) = 0 then y ∈ C 6+γ (G)) 0 Lemma 3. Assuming the following compatibility ∂ i+j f (0, 1) = 0, ∂x i ∂y j
0 i + j 5,
then w2∗ ∈ C 3,γ (Ω ∗∗∗ ), γ ∈ (0, 0.5), where w2∗ is the solution of (12e). Proof. Let w2∗ = z2∗ − yv2∗ (x, 1). Then z2∗ = 0 on the boundary ∂Ω ∗∗∗ and
L∗ z2∗ = L∗ w2∗ + yL∗ v2∗ (x, 1) = (w1∗ )xx + yL∗ v2∗ (x, 1) , (x, y) ∈ Ω ∗∗∗ .
(16)
E. O’Riordan, G.I. Shishkin / Applied Numerical Mathematics 58 (2008) 1761–1772
1767
Note that (w1∗ )xx (0, 1) = 0 and w2∗ is the solution of an elliptic problem. So to have sufficient compatibility at the corner (0, 1) we also require that L(v2∗ (x, 1))(0, 1) = 0. This follows from the fact that v2∗ ∈ C 3,γ (Ω) and v2∗ (0, y) = 0, 0 y 1. For sufficient regularity of the data we hence require that (w1∗ )xx ∈ C 1,γ (Ω ∗∗∗ ), γ ∈ (0, 0.5). Since w1∗ is a solution of a parabolic equation (where x is the time-like variable) we require that ∂ i+j w1∗ ∈ C γ (Ω ∗∗∗ ), ∂x i ∂y j
0 2i + j 6.
This places constraints on w0∗ . From (7), if we assume (16) then ∂ i (v0∗ + εv1∗ ) (0, 1) = 0, ∂x i
i 4.
On ΓT∗ , −w0∗ = v0∗ + εv1∗ so from the discussion on problem (15), the assumption (16) suffices for w0∗ ∈ C 8+2γ (Ω ∗∗∗ ). Recall that L∗p w1∗ = (w0∗ )xx ≡ f3 . By virtue of (16) ∂ i+2 (v0∗ + εv1∗ ) ∂ i f3 (0, 1) = (0, 1) = 0, i 2. ∂x i ∂x i+2 Note also that ε(f3 )yy = (a ∗ (w0∗ )x )xx which implies that
(17a)
∂ 2 f3 (0, 1) = 0. ∂y 2
(17b)
Note further that since w0 (0, y) = 0 ∂ ε ∂ 2 w0 ∂ 2 w0 (0, 1) = 0 (0, 1) = ∂y∂x ∂y a ∗ ∂y 2 and since w0∗ ∈ C 8+2γ (Ω ∗∗∗ ) we also have that ∂ 3 w0 ∂ 3 w0 ∂ 4 w0 (0, 1) = (0, 1) = (0, 1) = 0. ∂y 2 ∂x ∂y∂x 2 ∂y∂x 3 Combining these with ε 2 (f3 )yyyy = ε(a ∗ (w0 )x )xxyy implies that ∂f3 (0, 1) = 0, ∂y Finally ε
∂ 4 f3 ∂ 2 f3 (0, 1) = 0. (0, 1) = 4 ∂x∂y ∂y
∂ 5 w0 ∂3 ∗ ∂w0 , a = ∂x ∂y 5 ∂y 3
ε
∂ 7 w0 ∂5 ∗ ∂w0 , a = ∂y 7 ∂x ∂y 5
(17c)
ε
∂ 5 w0 ∂3 ∗ ∂w0 a = ∂x ∂x 3 ∂y 2 ∂x 3
will result in the following be satisfied ∂ 3 f3 ∂ 3 f3 (0, 1) = (0, 1) = 0. 2 ∂x∂y ∂y 3
(17d)
Together the identities in (17) yield sufficient compatibility (see (14)) for ∂ i+j w1∗ ∈ C 0+2γ (Ω ∗∗∗ ), ∂x i ∂y j
0 2i + j 6.
Note that for f sufficiently regular, f 0,γ Cf 0 + f 0+2γ , and thus ∂ i+j w1∗ ∈ C γ (Ω ∗∗∗ ), ∂x i ∂y j
0 2i + j 6.
2
Using the above expansion for wT we have the following bounds on the parabolic layer function associated with the edge y = 1.
1768
E. O’Riordan, G.I. Shishkin / Applied Numerical Mathematics 58 (2008) 1761–1772
Lemma 4. For the solution wT of (12) and i, j = 2, 3 √ wT (x, y) Ce−(1−y)/ ε , j i √ √ ∂ wT ∂ wT Cε −j/2 e−(1−y)/ ε , Cε 2−i e−(1−y)/ ε . ∂y j ∂x i
(18a) (18b)
Proof. By choosing the appropriate barrier function √ ∗ w (x, y) Ce αx e−(1−y)/ ε . 0
Note also that w0∗ = −(v0∗ + εv1∗ ) on ΓT∗ . Assuming sufficient regularity on a and f then we have that i ∗ ∂ (v0 + εv1∗ ) C, 0 i 4. ∂x i ∗ √ and noting that L∗ Introducing the stretched variable ζ = 1−y p w0 = 0 we deduce from the corresponding bounds to ε (3) for the parabolic problem (using the argument in [8, Theorem 4]) i+j ∂ w˜ 0 (x, ζ ) Ce−ζ , 0 2i + j 8. ∂x i ∂ζ j 0+2γ
This implies that i+j 2 ∂ ∂ w˜ 0 Ce−ζ , ∂x i ∂ζ j ∂x 2 0+2γ Noting that L∗p w1∗ = (w0∗ )xx , (x, y) ∈ G, x α
√ e−(1−y)/ ε .
0 2i + j 4. w1 = 0, (x, y) ∈ S and by choosing appropriate barrier function we have
that |w1 (x, y)| Ce Again using the argument from [8] we have for 0 2i + j 6 i+j ∂ w˜ 1 (x, ζ ) Ce−ζ , 0 2i + j 6, ∂x i ∂ζ j 0+2γ which implies that i+j 2 ∂ ∂ w˜ 1 Ce−ζ , ∂x i ∂ζ j ∂x 2 0+2γ
0 2i + j 2.
The above inequalities imply that i+j ∗ i+j ∗ √ ∂ w0 ∂ w1 −j/2 −(1−y)/ ε e , ∂x i ∂y j + ∂x i ∂y j Cε
(19)
0 i + j 3. 2x
√
Finally for ε sufficiently small (4ε < α 2 ) |w2∗ (x, y)| Ce α e−(1−y)/ ε . Note that from (19), we have that for γ ∈ (0, 0.5) and k = 0, 1, k = i + j i+j i+j (1−y) ∂ ∂ γ −j/2 − √ε ˜ ˜ + ε Cε e , 0 2i + j 2, f f ∂x i ∂y j ∂x i ∂y j k k+2γ where f˜ = (w1∗ )xx . From the bounds (2) and the derivative bounds on v2 , we deduce that
|w2∗ |2 C ε −1 |f˜|0 + ε −1+γ f˜ 0,γ + ε −2 |w2∗ |0 + Cε −2 Cε −2 ,
|w2∗ |3 C ε −2 |f˜|0 + ε −1 |f˜|1 + ε −1+γ f˜ 1,γ + ε −3 |w2∗ |0 + Cε −3 Cε −3 since for any f ∈ C 0+2γ , f 0,γ Cf 0 + C f 0+2γ . This completes the proof. The bottom parabolic layer component wB is defined analogously and
2
E. O’Riordan, G.I. Shishkin / Applied Numerical Mathematics 58 (2008) 1761–1772 y wB (x, y) Ce− √ε ,
1769
3 y ∂ wB Cε −3/2 e− √ε , (x, y) ∂y 3
i y ∂ wB Cε 2−i e− √ε , (x, y) ∂x i
i = 2, 3.
6. Corner layer functions Finally we define the outflow corner layer function wT R (and wBR is defined analogously) which is defined to be the solution of LwT R = 0,
(x, y) ∈ Ω = (0, 1) × (0, 1);
wT R (1, y) = −wT (1, y),
wT R (x, 0) = 0,
wT R (0, y) = 0,
wT R (x, 1) = −wR (x, 1).
Note that wT is compatible with u − v at the corner (1, 1) and, in turn, u − v is compatible with wR at the corner (1, 1). Hence wT R ∈ C 3,γ (Ω). From the comparison principle and the bounds on wT and wR , we have that for ε sufficiently small (1−y) α(1−x) wT R (x, y) Ce− 2ε e− √ε .
(20)
We can write this corner layer function as the following sum
wT R (x, y) = wT R (1, y)φ1 (x) + wT R (x, 1) − wT R (1, 1)φ1 (x) φ2 (y) + εzT R (x, y) where −εφ1 (x) + a(1, 1)φ1 (x) = 0,
φ1 (0) = 0,
φ1 (1) = 1,
φ2 (y) =
1 − e−y/
√
ε √ . 1 − e−1/ ε
Then zT R (x, y) = 0 on the boundary of Ω and noting that −ε
∂2 ∂ ∂2 wR (x, 1) + a(x, y) wR (x, 1) = ε 2 wR (x, 1) 2 ∂x ∂x ∂y
and using the bounds on wR corresponding to (10) we get that ε LzT R (x, y) Ce−α(1−x)/(4ε) . As before, we deduce the derivative bounds |zT R |k Cε −k , k = 1, 2, 3. One can then establish the bounds i i ∂ wT R Cε −i , ∂ wT R Cε −i , i = 2, 3. ∂x i ∂y i Also with the above decomposition of wT R we have the additional bound 3 ∂ wT R −2 ∂y 3 Cε .
(21)
(22)
Lemma 5. The corner layer component wT R satisfies the bounds on its derivatives given in (20), (21) and (22). 7. Discrete problem and error analysis We present error bounds for the approximations generated by using a standard upwind difference operator LN U = −ε(δx2 + δy2 )U + aDx− U = f,
(xi , yj ) ∈ Ω N
(23a)
on a mesh Ω N = ωx × ωy , which is a tensor product of two piecewise-uniform one-dimensional meshes ωx , ωy . The finite difference operators Dx− and δx2 are the standard first order backward difference and the second order centered difference on a non-uniform mesh [1]. Here the mesh ωx places N/2 mesh intervals into both [0, σx ] and [σx , 1] and
1770
E. O’Riordan, G.I. Shishkin / Applied Numerical Mathematics 58 (2008) 1761–1772
the mesh ωy places N/4 mesh intervals into the intervals [0, σy ] and [1 − σy , 1] and N/2 mesh intervals into the region [σy , 1 − σy ]. The transition parameters are taken to be √ ε σx = min 0.5, ln N (23b) and σy = min{0.25, ε ln N }. α The discrete solution is decomposed in an analogous fashion to the continuous solution. That is U = V + WR + WT + WBR + WT R where LN V = f, LN W = 0 (xi , yj ) ∈ Ω N and V = v, W = w (xi , yj ) ∈ ∂Ω N . On an arbitrary mesh using the bounds (6) one has the truncation error estimate 3 3 2 N L (v − V ) CN −1 ε ∂ v + ∂ v + CN −1 ε ∂ v CN −1 ∂y 3 ∂x 3 ∂x 2 and hence v − V Ω N CN −1 . If σx = 0.5 then the same argument coupled with the bounds (2) yields u − U Ω N CN −1 (ln N)2 . (If σy = 0.25 and σx < 0.5, then a minor modification to the argument below yields the same error bound.) In the remainder of this section we assume that √ ε σx = ln N and σy = ε ln N. α Note that from (8) and by the choice of transition point σx wR (x, y) CN −1 , x 1 − σx . Consider the discrete one-dimensional barrier function Ψ (xi ) which is the solution of −εδx2 Ψ (xi ) + αD − Ψ (xi ) = 0,
xi ∈ ω x ,
Ψ (0) = 0,
Ψ (1) = 1.
Note also that for the top y = 1 and bottom y = 0 edges, by (8), α(1−x ) WR (xi , yj ) = wR (xi , yj ) Ce− ε i CΨ (xi ) + CN −1 ,
yj = 0, 1.
Using CΨ (xi ) + CN −1 as a barrier function we deduce that |WR (xi , yj )| CΨ (xi ) + CN −1 and WR (xi , yj ) CN −1 , xi 1 − σx . Thus |(WR − wR )(xi , yj )| CN −1 , xi 1 − σx . In the fine mesh region (1 − σx , 1) × (0, 1), using the bounds (10) we have the truncation error bound 2 3 3 −1 N L (wR − WR ) CN −1 σx ε ∂ wR + ∂ wR + CN −1 ε ∂ wR C N ln N . ∂x 3 ∂x 2 ∂y 3 ε Use the barrier function
N −1 ln N C xi − (1 − σx ) + CN −1 ε in the region (1 − σx , 1) × (0, 1) to derive the error bound WR − wR Ω N CN −1 (ln N )2 . Next we examine the top parabolic layer function. Consider the discrete barrier function Υ (xi )Φ(yj ) where Φ(0) = 0, −εδy2 Φ(yj ) + Φ(yj ) = 0, yj ∈ ωy , i 2hj 1+ , hj = xj − xj −1 . Υ (xi ) = α
Φ(1) = 1,
j =1
Assuming the strict inequality 8ε < α 2 and for N sufficiently large (independently of ε) we have that −εδx2 Υ (xi ) + αDx− Υ (xi ) − Υ (xi ) 0,
xi ∈ ω x .
E. O’Riordan, G.I. Shishkin / Applied Numerical Mathematics 58 (2008) 1761–1772
1771
Note also that √ WT (1, yj ) Ce−(1−yj )/ ε CΦ(yj ) + CN −1 . Then |WT (xi , yj )| CΥ (xi )Φ(yj ) + CN −1 and hence |WT (xi , yj )| CN −1 , yj 1 − σy . In the fine mesh region (0, 1) × (1 − σy , 1) use the bounds (18) to obtain 2 3 3 N L (wT − WT ) CN −1 ε ∂ wT + ∂ wT + Cσy N −1 ε ∂ wT CN −1 ln N. ∂y 3 ∂x 3 ∂x 2 Hence wT − WT Ω N CN −1 ln N. As for the discrete regular and parabolic layer components we have that WT R (xi , yj ) CΨ (xi ) + CN −1 CN −1 , xi 1 − σx , WT R (xi , yj ) CΥ (xi )Φ(yj ) + CN −1 CN −1 , yj 1 − σy . In the fine mesh corner region (1 − σx , 1) × (1 − σy , 1), use the bounds given in (21) and (22) to obtain −1 N L (wT R − WT R ) C N ln N . ε
Use the barrier function C(εN )−1 ln N (xi − (1 − σx )) + CN −1 in the corner region (1 − σx , 1) × (1 − σy , 1) to derive the error bound WT R − wT R Ω N CN −1 (ln N )2 . Theorem 6. If U is the solution of (23) and u is the solution of the elliptic problem (1) then for 8ε < α 2 and N sufficiently large (independently of ε) U − uΩ N CN −1 (ln N )2 . Numerical results in [1] agree with this theoretical bound on the nodal error. Furthermore, numerical experiments in [1] demonstrate that √ the numerical method (23) also produces parameter-uniform approximations to the scaled derivatives εux and εuy . 8. Conclusions A proof of parameter-uniform pointwise convergence of a set of numerical approximations to a class of singularly perturbed elliptic equations is given. The tools of analysis used involve a special decomposition of the solution into a sum of a regular component and several types of layer components. In particular, the solutions contain parabolic boundary layers. Parameter-explicit pointwise bounds on the derivatives of all of these components are established. The establishment of these bounds on the derivatives is the main component in the proof of first order (up to logarithmic factors) parameter-uniform convergence of the numerical approximations. References [1] P.A. Farrell, A.F. Hegarty, J.J.H. Miller, E. O’Riordan, G.I. Shishkin, Robust Computational Techniques for Boundary Layers, Chapman and Hall/CRC Press, Boca Raton, FL, 2000. [2] H. Han, R.B. Kellogg, Differentiability properties of solutions of the equation −ε 2 u + ru = f (x, y) in a square, SIAM J. Math. Anal. 21 (1990) 394–408. [3] P.W. Hemker, A singularly perturbed model problem for numerical computation, J. Comput. Appl. Math. 76 (1996) 277–285. [4] R.B. Kellogg, M. Stynes, Corner singularities and boundary layers in a simple convection diffusion problem, J. Differential Equations 213 (2005) 81–120. [5] O.A. Ladyzhenskaya, N.N. Ural’tseva, Linear and Quasilinear Elliptic Equations, Academic Press, New York and London, 1968. [6] O.A. Ladyzhenskaya, V.A. Solonnikov, N.N. Ural’tseva, Linear and Quasilinear Equations of Parabolic Type, Transactions of Mathematical Monographs, vol. 23, American Mathematical Society, Providence, RI, 1968. [7] T. Linß, M. Stynes, Asymptotic analysis and Shishkin-type decomposition for an elliptic convection–diffusion problem, J. Math. Anal. Appl. 261 (2001) 604–632. [8] J.J.H. Miller, E. O’Riordan, G.I. Shishkin, L.P. Shishkina, Fitted mesh methods for problems with parabolic boundary layers, Math. Proc. Roy. Irish Acad. A 98 (2) (1998) 173–190.
1772
E. O’Riordan, G.I. Shishkin / Applied Numerical Mathematics 58 (2008) 1761–1772
[9] H.-G. Roos, Optimal convergence of basic schemes for elliptic boundary value problems with strong parabolic layers, J. Math. Anal. Appl. 267 (2002) 194–208. [10] S. Shih, R.B. Kellogg, Asymptotic analysis of a singular perturbation problem, SIAM J. Math. Anal. 18 (1987) 1467–1511. [11] G.I. Shishkin, Discrete Approximation of Singularly Perturbed Elliptic and Parabolic Equations, Russian Academy of Sciences, Ural Section, Ekaterinburg, 1992 (in Russian). [12] E.A. Volkov, Differentiability properties of solutions of boundary value problems for the Laplace and Poisson equations, Proc. Steklov Inst. Math. 77 (1965) 101–126.