WEIGHTED-NORM FIRST-ORDER SYSTEM LEAST SQUARES ...
WEIGHTED-NORM FIRST-ORDER SYSTEM LEAST SQUARES (FOSLS) FOR PROBLEMS WITH CORNER SINGULARITIES E. LEE∗ , T. A. MANTEUFFEL∗‡ , AND C. R. WESTPHAL† § Abstract. A weighted-norm least-squares method is considered for the numerical approximation of solutions that have singularities at the boundary. While many methods suffer from a global loss of accuracy due to boundary singularities, the least-squares method can be particularly sensitive to a loss of regularity. The method we describe here requires only a rough lower bound on the power of the singularity and can be applied to a wide range of elliptic equations. Optimal order discretization accuracy is achieved in weighted H 1 , and functional norms and L2 accuracy is retained for boundary value problems with a dominant div/curl operator. Our analysis, including interpolation bounds and several Poincar´ e type inequalities, are carried out in appropriately weighted Sobolev spaces. Numerical results confirm the error bounds predicted in the analysis.
1. Introduction. In this paper, we develop a method for treating div/curl systems with reduced regularity. While motivated by first-order systems that arise from second-order elliptic boundary value problems, div/curl systems appear in many contexts, for example, in Maxwell’s equations and in the vorticity from of Stokes equations. These problems have the fortunate property of a guaranteed smooth solution as long as the data and domain are smooth. However, many problems of interest are posed in non-smooth domains and, as a consequence, lose this property at a finite number of points on the boundary in two dimensions or along curves on the boundary in three dimensions. In the present paper, we study two-dimensional problems that have non-smooth solutions at irregular boundary points, that is, points that are corners of polygonal domains, locations of changing boundary condition type, or both. Similar behavior occurs in problems with discontinuous material coefficients and the methods presented here can easily be extended to that situation. Standard solution techniques that attempt to approximately solve a div/curl system with reduced regularity using H 1 -conforming finite elements will, in general, fail to converge. This phenomenon can be explained by noting that the Sobolev space (H 1 )2 is a closed subspace of H(div) ∩ H(curl) (see section 2), which implies either that (H 1 )2 = H(div) ∩ H(curl), the case of full regularity, or (H 1 )2 is a proper subspace. In this case, the co-dimension is finite and is spanned by so-called singular functions. In the presence of reduced regularity, the solution will, in general, not be in (H 1 )2 . A standard finite element method using H 1 -conforming elements will converge to the element of (H 1 )2 closest to the true solution in the H(div) ∩ H(curl) norm. Local mesh refinement will not alter this outcome. If a basis for the singular functions is known, it can be incorporated directly into the finite element space (c.f.[10, 16, 24]). In [3, 4], this approach is applied to div/curl systems and shown to restore optimal convergence throughout the domain at a minimal additional cost. For some two-dimensional problems, the singular basis functions are known and can be included in the finite element space. For the other ∗ Department of Applied Mathematics, University of Colorado at Boulder; Boulder,CO 80309– 0526, USA. Email: {tmanteuf,eunjung}@colorado.edu. † Department of Mathematics and Computer Science, Wabash College, PO Box 352, Crawfordsville, IN 47933, USA. Email: [email protected]. ‡ This work was sponsored by the National Science Foundation under grant number DMS– 0420873, and the Department of Energy under grant numbers DE–FC02-01ER25479 and DE–FG0203ER25574. § This work was sponsored by the National Science Foundation under grant DMS–9810751
1
2
LEE, MANTEUFFEL, AND WESTPHAL
problems, or in three dimensions, the exact character of the singular functions is less well understood which makes this approach more difficult to implement. As a different type of remedy for this so-called pollution effect, least-squares methods based on inverse norms can be effective for problems with irregular boundary points, discontinuous coefficients, and data in H −1 . For example, in [5, 8, 11–13, 19], the functional is posed in terms of H −1 -norms rather than L2 -norms, resulting in optimal L2 -approximations to the solution. A more recent approach, called FOSLL∗ , uses an inverse norm induced by the equations, and is shown in [14, 22] to be more efficient than the H −1 -norm methods. Graded mesh refinement in weighted Sobolev spaces has been shown to be effective in restoring optimal convergence, in L2 - and H 1 -norms, in the context of a Galerkin formulation of second-order elliptic problems with reduced regularity (c.f. [2], [1]). However, in that context, the solution is in H 1 . Convergence would occur, although more slowly, without weighting and mesh refinement. In [15], a weighted norm and a sequence of graded meshes is used in an H(div) least-squares functional, arising from a second-order elliptic problem, to restore optimal convergence for the primal variable, in both L2 - and H 1 -norms, if the flux variable is approximated in a finite element space satisfying the grid decomposition property, for example, Raviart-Thomas elements (c.f.[6]). The flux variable converges in a weighted H(div)-norm. In this paper we examine div/curl systems that lack regularity within the firstorder system least-squares (FOSLS) framework. The basic FOSLS approach is to recast the original system as an appropriate first-order system and apply an L 2 minimization principle over the residual of the equations. If possible, this reformulation is done by minimizing a functional whose quadratic part is equivalent to the product H 1 -norm, indicating that the process is similar to solving a weakly coupled system of Poisson-like equations. This equivalence also guarantees optimal H 1 -accuracy for standard discretizations. For div/curl systems with reduced regularity, as briefly mentioned above, the L2 based functional fails to be H 1 -equivalent and, as a consequence, standard discretizations suffer from the pollution effect. Here, a weighted-norm leastsquares method is developed that restores optimal convergence using H 1 -conforming finite elements without graded mesh refinement. It replaces the L2 -norms in the FOSLS functional with weighted L2 -norms, making the functional norm equivalent to a weighted H 1 -norm. With an appropriate weighting function, this method recovers optimal order accuracy in the weighted L2 - and H 1 -norms, and retains optimal L2 convergence even near the singularity. Our method requires only the power of the singularity (not the actual singular solution) to be known a priori and, in practice, can be used with only a rough estimate of the power of the singularity, which can be adaptively determined if unknown (c.f.[4]). The method developed here has some similarity to [15], but considers an H(div)∩ H(curl) functional and considers more general boundary conditions. Most importantly, our analysis admits a more aggressive weighting, resulting in optimal order accuracy in the weighted norms without mesh refinement. In addition, we prove several Poincar´e type inequalities in weighted Sobolev spaces under a variety of boundary conditions that, in addition to being necessary for our main result, may be of independent interest. We use Poisson’s equation on a domain with a reentrant corner as a model problem and as the formal setting for analysis. The resulting div/curl system is the focus. The analysis in this paper is restricted to two dimensions. However, the approach suggests
FOSLS FOR PROBLEMS WITH SINGULARITIES
3
a natural generalization to problems with reduced regularity due to discontinuous coefficients and to problems in three dimensions. This paper is organized as follows. Section 2 contains notation and preliminary discussion. In section 3, the weighted FOSLS functionals are described. Section 4 contains Poincar´e inequalities and regularity results in weighted Sobolev spaces. Error bounds for the weighted FOSLS functional are presented in section 5 and section 6 contains computational results. 2. Singular solutions and preliminaries. For vector function u = (u 1 , u2 )t , let the divergence and curl of u be defined in the standard way: ∇ · u = ∂x u1 + ∂y u2 and ∇ × u = ∂x u2 − ∂y u1 . Further, define the formal adjoint of the curl operator by ∂y q . ∇⊥ q = −∂x q We use standard notation for Sobolev spaces H k (Ω)d , corresponding inner product (·, ·)k,Ω , and norm k · kk,Ω , for k ≥ 0. We drop subscript Ω and superscript d when the domain and dimension are clear by context. Since H 0 (Ω) coincides with L2 (Ω) we often denote k · k0 by k · k. Define the subspaces of L2 (Ω) induced by the divergence and curl of u by H(div) = {u ∈ L2 (Ω) : k∇ · uk < ∞},
H(curl) = {u ∈ L2 (Ω) : k∇ × uk < ∞}. We also make use of the following general inequalities for nonnegative a and b: |a|2 + |b|2 ≤ |a + b|2 ≤ 2(|a|2 + |b|2 ).
(2.1)
Consider the function f (r, θ) = r a in two dimensional polar coordinates. Assume that the origin lies on the boundary of domain Ωw = {(r, θ) : 0 < r < R, 0 < θ < ω < 2π}, as pictured in figure 2.1. By a direct computation it is clear that f ∈ H k (Ω) only for k < a + 1.
ω R
Fig. 2.1. Simple wedge-shaped domain, Ωw .
Now, consider Poisson’s equation on a domain in R2 with a corner of interior angle ω. It is well known that, for the case of Dirichlet or Neumann boundary conditions, the solutions of this boundary value problem may include those with radial part of
4
LEE, MANTEUFFEL, AND WESTPHAL π
the form p ∼ r ω in a local polar coordinate system centered at the corner. Thus, for the case of reentrant corners, ω > π, the solution fails to be in H 2 (Ω) and we say that the problem has a singularity (or singular solution). For problems with Dirichlet and Neumann boundaryπ conditions meeting at the corner, solutions may have components of the form p ∼ r (2ω) . Thus, for mixed boundary conditions, singularities may occur at corners with ω > π/2. We now explore this issue in more detail. Define the power of the singularity to be α = π/ω for Dirichlet or Neumann boundary conditions and α = π/(2ω) for mixed boundary conditions. The solution to Poisson’s equation may be written as p(r, θ) = p0 (r, θ) + s(r, θ), where p0 (r, θ) ∈ H 2 (Ω) and s(r, θ) ∈ H 1+m (Ω) for m < α. The singular part of the solution has the form s(r, θ) = r α (κ1 sin(αθ) + κ2 cos(αθ)), where the values of κ1 and κ2 depend on boundary conditions (see [17, 18]). For the FOSLS formulation of this problem, we may similarly decompose unknown u = ∇p as u(r, θ) = u0 (r, θ) + ∇s(r, θ), where ∇s(r, θ) has the form ∇s(r, θ) = αr α−1
κ1 sin(α − 1)θ + κ2 cos(α − 1)θ . κ1 cos(α − 1)θ − κ2 sin(α − 1)θ
Thus, the unknown u(r, θ) is in H k (Ω) only for k < α. For example, consider Poisson’s equation posed on the simple domain in figure 2.1. Let the solution to this boundary value problem in polar coordinates be 2 p = χ(r)r 3 sin(2θ/3), where χ(r) is a smooth transition function that is 1 on a platform near the origin and vanishes at the boundaries not adjacent to the origin. Then, p = 0 on ∂Ω and ∆p =
1 1 2 ∂r (r∂r p) + 2 ∂θθ p r r 2
= r 3 χ00 (r) + 73 r
−1 3
χ0 (r) sin(2θ/3),
and, thus, it is clear that ∆p ∈ L2 (Ω), but p ∈ / H 2 (Ω). We say this problem fails to 2 provide full lifting of the data (from L (Ω) to H 2 (Ω) e.g.). The solution, u = ∇p, is, thus, not in H 1 (Ω). 3. Weighted-norm least squares. As before, let Ω be a domain with a corner of interior angle ω at the origin, and we may, without loss of generality, further assume diam(Ω) ≤ 1. For f ∈ L2 (Ω), let p satisfy in Ω, −∆p = f, p = 0, on ΓD , (3.1) n · ∇p = 0, on ΓN ,
FOSLS FOR PROBLEMS WITH SINGULARITIES
5
¯D ∪ Γ ¯ N . When this problem where n is the outward unit normal to Ω and ∂Ω = Γ 2 is H regular, the normal FOSLS methodology is to introduce the new unknown, u = ∇p, and rewrite system (3.1) as u − ∇p = 0, in Ω, −∇ · u = f, in Ω, in Ω, ∇ × u = 0, τ · u = 0, on ΓD , (3.2) n · u = 0, on Γ , N p = 0, on Γ D, n · ∇p = 0, on ΓN . Here, τ is the counter-clockwise unit tangental vector to Ω. Since this system can be posed completely in terms of u, we may decouple the equations in (3.2), solve for u first, and then recover p from u. To this end, define the two L2 -norm functionals, G(u; f ) = k∇ · u + f k2 + k∇ × uk2 ,
G2 (p; u) = ku − ∇pk2 , and the spaces,
V = {v ∈ H 1 (Ω) : τ · v = 0 on ΓD , n · v = 0 on ΓN },
W = {q ∈ H 1 (Ω) : q = 0 on ΓD }.
Thus, the two-stage solution process is to minimize G(v; f ) over V and then, given the minimizer, u, minimize G2 over W: (1)
G(u; f ) = inf G(v; f ),
(2)
G2 (p; u) = inf G2 (q; u).
v∈V
(3.3)
q∈W
The goal of the FOSLS methodology is, generally, to formulate functionals whose quadratic part is equivalent to the H 1 norm whenever possible. The second stage functional is H 1 equivalent and the solution we seek is always in H 1 . The first stage functional, however, is not always H 1 equivalent. For domains with reentrant corners, there is no H 1 sequence of functions that converges to the solution in the H(div) ∩ 2 H(curl) norm. To illustrate, consider the example above where p = χ(r)r 3 sin(2θ/3) and u = ∇p. A simple computation reveals that ∇·u, ∇×u ∈ L2 (Ω), but u ∈ / H 1 (Ω). Define the weighted functional by Gw (u; f ) = kw(∇ · u + f )k2 + kw∇ × uk2 ,
(3.4)
where the weight function has the form w = r β for some β > 0. Define the weighted Sobolev norm, k · kk,β , on Ω in terms of the standard L2 norm, k · k0 , by kqkk,β = (
X
|j|≤k
1/
krβ−k+|j| Dj qk20 ) 2 ,
(3.5)
6
LEE, MANTEUFFEL, AND WESTPHAL
where D j is the standard distributional derivative of order j. Similarly, define the weighted seminorm by |q|k,β = (
X
|j|=k
krβ−k+|j| Dj qk20 )
1/2
(3.6)
and the associated weighted Sobolev space by Hβk (Ω) = {q : kqkk,β < ∞}. Define the div/curl operator, L, and vector f by L = may now write the weighted functional from (3.4) as
(3.7)
∇· ∇×
and f =
f . We 0
Gw (u; f ) = kLu − f k20,β . The weighted-norm least-squares minimization problem for the first-stage solution is then: find u ∈ V such that Gw (u; f ) = inf Gw (v; f ). v∈V
The second-stage solution for p remains as described above. We seek values of β that make H 1 (Ω) dense in H(div) ∩ H(curl) in the weighted functional norm and result in the most accurate discretizations possible. For the discrete problem, we may choose any finite dimensional subset of H 1 over which to minimize the weighted functional. Let P h denote the space of C 0 piecewise polynomial (or tensor product) elements on triangles (or quadrilaterals) of meshsize h and V h the subspace of P h that satisfies the appropriate boundary conditions on Ω: V h = {vh ∈ P h : τ · vh = 0 on ΓD , n · vh = 0 on ΓN }. The discrete weighted-norm least-squares minimization problem is, then, to minimize the discrete functional: find uh ∈ V h such that Gw (uh ; f ) = min Gw (vh ; f ). vh ∈V h
(3.8)
By unweighting the equations near the singularity, the functional is freed from trying to approximate the solution (which is not in H 1 (Ω)) in the H 1 sense near the singularity. But, away from the singularity, the weighted functional retains the same character as the normal non-weighted functional. We now consider the choice of weight parameter β and its relation to weighted and nonweighted a priori error bounds on the approximated solution. 4. Poincar´ e Bounds and Regularity Estimates. In this section, we establish several theoretical results in weighted Sobolev spaces and error bounds for the weighted-norm method. Here, we establish several Poincar´e bounds in the domain Ωw . We first prove a result for the scalar pure Neumann and pure Dirichlet problems, and then for the scalar mixed boundary condition problem. These results lead to a Poincar´e inequality for the vector case.
7
FOSLS FOR PROBLEMS WITH SINGULARITIES
Lemma 4.1 Take Ω = Ωw and let β > 0, > 0 and γ = β − 3/2 − . Further, assume that γ 6= −2. For functions q ∈ Hβ1 (Ω), with r β− ∇q, r γ+1 ∇q ∈ L2 (Ω), that can be chosen to satisfy ZZ rγ q(r, θ) r dr dθ = 0, (4.1) Ω
we have the bound: kqk0,β−1 ≤ C()k∇qk0,β− ,
(4.2)
where C() depends on , β, and Ω and C() → ∞ as → 0. Proof. For any two points (r, θ) and (r0 , θ0 ) in Ω, write q(r, θ) as q(r, θ) = q(r0 , θ0 ) +
Z
r r0
∂q(ˆ r, θ) dˆ r+ ∂ rˆ
Z
θ θ0
ˆ ∂q(r0 , θ) ˆ dθ. ∂ θˆ
Multiplying both sides of the equation by r0γ+1 , integrating with respect to r0 and θ0 over Ω, and Fubini’s theorem yield Rγ+2 ω q(r, θ) γ+2
=
Z ωZ RZ 0
0
r0
∂q(ˆ r, θ) r0γ+1 ∂ rˆ
dˆ r dr0 dθ0 +
Z ωZ RZ 0
0
θ θ0
r0γ+1
ˆ ∂q(r0 , θ) dθˆ dr0 dθ0 ∂ θˆ
Z RZ R ∂q(ˆ r, θ) ∂q(ˆ r, θ) dr0 dˆ r−ω dr0 dˆ r r0γ+1 ∂ rˆ ∂ rˆ 0 r 0 rˆ Z RZ ωZ ω Z RZ θZ θˆ ˆ ˆ ∂q(r0 , θ) γ+1 ∂q(r0 , θ) ˆ dθ0 dθ dr0 − dθ0 dθˆ dr0 + r0 r0γ+1 ∂ θˆ ∂ θˆ 0 θ θˆ 0 0 0
= ω
=
Z
r
r
Z
rˆ
r0γ+1
Z R Z ωRγ+2 R ∂q(ˆ ω ∂q(ˆ r, θ) r, θ) dˆ r− dˆ r rˆγ+2 γ+2 0 ∂ rˆ γ+2 r ∂ rˆ Z RZ ω Z RZ ω ˆ ˆ ∂q(r0 , θ) γ+1 ∂q(r0 , θ) ˆ ˆ + dθ dr0 − ω dθˆ dr0 . θ r0 r0γ+1 ∂ θˆ ∂ θˆ 0 0 0 0
Note that above q is, by assumption, sufficiently smooth to apply Fubini’s theorem. By the triangle inequality we have that γ+2 R ω γ + 2 q(r, θ) Z R Z ω r, θ) ωRγ+2 R ∂q(ˆ r, θ) γ+2 ∂q(ˆ ≤ rˆ r+ dˆ r ∂ rˆ dˆ γ+2 0 γ + 2 r ∂ rˆ Z RZ ω ∂q(r , θ) ˆ 0 + 2ω r0γ+1 dθˆ dr0 . ∂ θˆ 0 0
8
LEE, MANTEUFFEL, AND WESTPHAL
Now, squaring each side and using (a + b + c)2 ≤ 3(a2 + b2 + c2 ) we get, !2 !2 Z R Z R ∂q(ˆ r, θ) r, θ) 3 2 γ+2 ∂q(ˆ r +3 r rˆ |q(r, θ)| ≤ 2(γ+2) ∂ rˆ dˆ ∂ rˆ dˆ R 0 r !2 Z RZ ω ˆ 12(γ + 2)2 γ+1 ∂q(r0 , θ) ˆ + r0 dθ dr0 . ∂ θˆ R2(γ+2) 0 0
(4.3)
Multiply each side of (4.3) by r 2β−1 and integrate with respect to r and θ over Ω. We consider each of the terms on the resulting right-hand side separately. First, !2 Z ωZ R Z R ∂q(ˆ r , θ) dˆ dr dθ r2β−1 rˆγ+2 r ∂ r ˆ 0 0 0 2 Z ωZ RZ R r, θ) 2β−1 2γ+4 ∂q(ˆ r rˆ ≤ R r dr dθ ∂ rˆ dˆ 0 0 0 2 Z Z r, θ) R2β+1 ω R 2γ+3 ∂q(ˆ rˆ = r dθ ∂ rˆ rˆ dˆ 2β 0
0
R2β+1 k∇qk20,γ+ 3 . ≤ 2 2β
We now consider the second term in (4.3). Since by the Schwarz inequality, !2 Z R ∂q(ˆ r, θ) r ∂ rˆ dˆ r !2 Z R r, θ) (−1)/2 (1−)/2 ∂q(ˆ = rˆ rˆ r ∂ rˆ dˆ r 2 Z R Z R ∂q(ˆ r, θ) dˆ r rˆ1− rˆ−1 dˆ r ≤ ∂ rˆ r r 2 Z r, θ) R R 1− ∂q(ˆ rˆ r, = ∂ rˆ dˆ r
we can apply Fubini’s theorem to bound the second term: !2 Z R Z ωZ R ∂q(ˆ r , θ) 2β−1 r r dr dθ ∂ rˆ dˆ r 0 0 2 Z Z Z r, θ) R ω R R 2β−1 1− ∂q(ˆ r dr dθ ≤ r rˆ ∂ rˆ dˆ 0 0 r 2 Z Z Z R ω R rˆ 2β−1 1− ∂q(ˆ r, θ) = r dθ r rˆ ∂ rˆ dr dˆ 0 0 0 2 Z ωZ R ∂q(ˆ r, θ) R = rˆ dˆ r dθ rˆ2β− 2β 0 0 ∂ rˆ R = k∇qk20,β− 2 2β R k∇qk20,β− , ≤ 2β
9
FOSLS FOR PROBLEMS WITH SINGULARITIES
The third term can be bounded similarly: !2 Z RZ ω Z ωZ R ∂q(r , θ) ˆ 0 γ+1 r0 r2β−1 dr dθ dθˆ dr0 ∂ θˆ 0 0 0 0 2 ! Z ωZ R Z RZ ω ∂q(r , θ) ˆ 0 2γ+2 ≤ r0 r2β−1 dr dθ Rω dr0 dθˆ ∂ θˆ 0 0 0 0 2 Z Z ˆ R2β+1 ω 2 ω R 2γ+3 1 ∂q(r0 , θ) r0 = r0 dr0 dθˆ r0 2β ∂ θˆ 0 0 =
R2β+1 ω 2 k∇qk20,γ+ 3 . 2 2β
Putting the three terms together and substituting γ = β − 3/2 − we may now write the bound Z ωZ R 2 r2β−2 |q(r, θ)| r dr dθ 0 0 −2 6(β + 21 − )2 R−2 ω 2 3R 3R ≤ + + k∇qk20,β− , 2β 2β β and the lemma follows by taking the square root of both sides. In what follows, let χ be a smooth function of r where χ = 1 for r < η and χ = 0 for r > 2η. We take η to be sufficiently small to ensure that supp(χ) ⊂ Ω. Lemma 4.2 Take Ω = Ωw and let q be a scalar function in Hβ1 (Ω), where β > 0. The following bound holds for χ as defined above: kχqk0,β−1 ≤
1 k∇(χq)k0,β . β
(4.4)
Proof. Hardy’s inequality for f (t) defined for t > 0 with lim f (t) = 0 gives t→0
(see [20]): Z
∞ 0
f2 dt ≤ 4 t2
Z
∞ 0
2
|f 0 | dt.
(4.5)
The lemma follows after a change of variables, t = r −2β , a substitution f (r) = χq(r, θ) for fixed θ, and an integration on both sides with respect to θ. Lemma 4.3 Take Ω = Ωw and let either q ∈ Hβ1 (Ω) with q = 0 on ∂Ω or q ∈ Hβ1 (Ω)/R with n · ∇q = 0 on ∂Ω. Then kqk0,β−1 ≤ Ck∇qk0,β for β > 0 where C depends only on Ω, and β. Proof. First, if q = 0 on ∂Ω, write q = q(r, θ) as q(r, θ) =
Z
θ 0
ˆ ∂q(r, θ) ˆ dθ. ∂ θˆ
(4.6)
10
LEE, MANTEUFFEL, AND WESTPHAL
Square both sides and multiply by r 2β−1 : r
2β−1
Z 2 θ ∂q(r, θ) ˆ dθˆ |q(r, θ)| = r ˆ 0 ∂θ 2 Z ω ˆ 1 ∂q(r, θ) 2β+1 ˆ ≤r ω dθ. ∂ θˆ 0 r 2
2β−1
Integrate both sides with respect to r and θ over Ω: Z
ω 0
Z
R 0
2 1 ∂q(r, θ) ˆ r2β+1 r2β−2 |q(r, θ)| r dr dθ ≤ ω dθˆ dr dθ ∂ θˆ 0 r 0 0 2 Z ω Z R0 1 ∂q(r, θ) ˆ ˆ r2β ≤ ω2 r dr dθ. r ∂ θˆ 0 0 Z
2
ω
Z
R0
Z
ω
The lemma follows since the right side is bounded by Ck∇qk20,β . Now, if q ∈ Hβ1 (Ω)/R then it may be chosen to satisfy ZZ rγ q r dr dθ = 0
(4.7)
Ω
for γ chosen as in lemma 4.1. By the triangle inequality and lemma 4.2 we get kqk0,β−1 ≤ kχ qk0,β−1 + k(1 − χ)qk0,β−1 ! 12 Z ω Z R 2 r 2β−2 2 2 ≤ Ck∇(χq)k0,β + r (1 − χ) q r dr dθ η 0 η ≤ C (k∇(q)k0,β + kqk0,β ) . Apply lemma 4.1 with = 1 to the kqk0,β term on the right side and the lemma follows. For the problem with mixed boundary conditions, consider Ωw partitioned into subdomains Ω0 = {(r, θ) : r ≤ 21 R0 , 0 ≤ θ ≤ ω} and Ω1 = Ω\Ω0 , as shown in figure 4.1.
Ω1 Ω0 1 R 2 0
R0 R
Fig. 4.1. Wedge-shaped domain, Ωw , partitioned into subdomains Ω0 and Ω1 .
FOSLS FOR PROBLEMS WITH SINGULARITIES
11
Lemma 4.4 Consider domain Ω = Ωw , as pictured in figure 4.1, and let q ∈ Hβ1 (Ω) vanish on the line segment of ∂Ω corresponding to θ = 0 and r < R0 . Then, there is a constant, C, dependent only on Ω, β, and R0 , such that kqk0,β−1 ≤ Ck∇qk0,β
(4.8)
for β > 0. Proof. For points (r, θ) in Ω0 we may derive the bound, kqk0,β−1,Ω0 ≤ Ck∇qk0,β,Ω0 ,
(4.9)
completely analogous to the proof of lemma 4.3. Now, consider points (r, θ) in Ω 1 . We may write q = q(r, θ) as q(r, θ) =
Z
r r˜
∂q(ˆ r, θ) dˆ r+ ∂ rˆ
Z
θ 0
ˆ ∂q(˜ r, θ) ˆ dθ, ∂ θˆ
where the point (˜ r, 0) is on the part of ∂Ω1 where q vanishes. By the Schwarz inequality, the triangle inequality and inequality (2.1) we have the bound 1 |q(r, θ)| ≤ 2(R − R0 ) 2 2
Z
r r˜
2 2 Z θ ˆ ∂q(ˆ r, θ) r, θ) ∂q(˜ ˆ r + 2ω dθ. ∂ rˆ dˆ ∂ θˆ 0
(4.10)
We now expand the limits in the integrals, multiply each side by r 2β−1 , integrate with respect to r over ( 21 R0 , R), integrate with respect to θ over (0, ω), and integrate with respect to r˜ over ( 12 R0 , R0 ), by Fubini’s theorem: 1 ( R0 ) 2
Z
ω 0
Z
R 1 2 R0
r2β−1 |q(r, θ)|2 dr dθ
2 Z R0Z ωZ R Z R 1 r, θ) 2β−1 ∂q(ˆ ≤ 2(R − R0 ) r r dr dθ d˜ r ∂ rˆ dˆ 1 1 1 2 2 R0 0 2 R0 2 R0 2 Z R0Z ωZ R Z θ ˆ r, θ) ∂q(˜ r. + 2ω dθˆ dr dθ d˜ 1 1 ∂ θˆ 0 2 R0 0 2 R0
(4.11)
We use the inequalities,
1 R0 ≤ r˜ ≤ R0 , 2
r˜ ≤ rˆ ≤ r ≤ R,
to derive the following simple bounds: 1≤
2ˆ r R0
,
r≤
2R R0
rˆ,
r≤
2R R0
r˜.
(4.12)
12
LEE, MANTEUFFEL, AND WESTPHAL
By applying the bounds in (4.12) and Fubini’s theorem, we may now write (4.11) as 1 ( R0 ) 2
Z
ω 0
Z
R 1 2 R0
r2β−2 |q(r, θ)|2 r dr dθ
2 2 2β−1 r, θ) 2ˆ r 2R 2β−1 ∂q(ˆ r dθ rˆ ∂ rˆ dˆ 1 R0 R0 0 2 R0 2 2β−1 Z R0 Z ω ∂q(˜ ˆ r , θ) 1 2R + 2ω 2 (R − R0 ) r r˜2β−1 dθˆ d˜ 1 ∂ θˆ 2 R 0 0 R 0 2 2 Z ωZ R r, θ) 1 2β+1 2β ∂q(ˆ 2 −2β 2β−1 rˆ dˆ r dθ ≤2 (R − R0 ) R0 R rˆ 1 2 ∂ rˆ 0 2 R0 2 2β−1 Z R Z ω 1 ∂q(˜ ˆ 1 r , θ) R r, + 22β ω 2 (R − R0 ) r˜2β r˜ dθˆ d˜ 1 r˜ ∂ θˆ 2 R0 0 R 0 2
1 ≤ (R − R0 )2 R0 2
Z
ω
Z
R
which directly implies
kqk0,β−1,Ω1 ≤ Ck∇qk0,β,Ω1 ,
(4.13)
where 1 ω2 1 C = 22β+1 (R − R0 )R2β−1 R0−2β−1 2(R − R0 ) + . 2 2 R0 Combining inequalities (4.9) and (4.13) completes the lemma. We now consider a similar Poincar´e inequality for the vector case. Again, consider Ω = Ωw , where ∂Ω is partitioned into Dirichlet and Neumann boundaries, ΓD and ΓN respectively. The following lemma is valid for the pure Dirichlet and Neumann cases and for the mixed boundary condition cases when ΓD includes a part of the boundary adjacent to the origin and ω 6= 3π 2 . Lemma 4.5 Take Ω = Ωw and let u ∈ Hβ1 (Ω)2 satisfy τ · u = 0 on ΓD and n · u = 0 on ΓN . Assume that for the mixed boundary condition case that ω 6= 3π/2. Then there is a constant, C, dependent only on Ω, β and the length of the segments of Γ D and ΓN adjacent to the origin, such that kuk0,β−1 ≤ Ck∇uk0,β
(4.14)
for β > 0. Proof. First, consider the case when τ · u = 0 on ∂Ω. Denote the part of ∂Ω aligned with θ = 0 as Γ1 and the part of ∂Ω aligned with θ = ω as Γ2 . Thus, u1 = 0 on Γ1 and τx u1 + τy u2 = 0 on Γ2 . Since u1 and τx u1 + τy u2 satisfy the conditions in lemma 4.4, we may use ku1 k0,β−1 ≤ Ck∇u1 k0,β and kτx u1 + τy u2 k0,β−1 ≤ Ck∇(τx u1 + τy u2 )k0,β .
FOSLS FOR PROBLEMS WITH SINGULARITIES
13
Further, take τy 6= 0, since τy = 0 corresponds to either ω = π, for which the result holds trivially since the boundary is smooth, or ω = 2π, which we do not consider. Now, kuk20,β−1 = ku1 k20,β−1 + ku2 k20,β−1 1 = ku1 k20,β−1 + 2 kτx u1 − τx u1 + τy u2 k20,β−1 τy ≤ (1 + 2
2 τx2 )ku1 k20,β−1 + 2 kτx u1 + τy u2 k20,β−1 τy2 τy
≤ C(k∇u1 k20,β + k∇(τx u1 + τy u2 )k20,β ) ≤ Ck∇uk20,β .
The case when n·u = 0 on ∂Ω is analogous since u2 = 0 on Γ1 and nx u1 +ny u2 = 0 on Γ2 . Also, when ω 6= π2 , 3π 2 , the case for mixed boundary conditions follows similarly using the result of lemma 4.4. The case for mixed boundary conditions when ω = π/2, follows from appealing to symmetry in the the pure Dirichlet problem for ω = π. Remark 4.6 Lemma 4.5 can be directly extended to more generally shaped domains. The proof of the scalar Poincar´e bounds in lemmas 4.1, 4.2, 4.3 and 4.4 are simplified when the domain has the shape of Ωw with only one irregular boundary point. Since we are primarily interested in a local result, proving lemma 4.5 in the simple domain is sufficient for our purposes. Consider the following scalar Poisson problem ∆p = f, in p = 0, on n · ∇p = 0, on
in Ωw : Ω, ΓD ,
(4.15)
ΓN .
We refer to system (4.15) as the pure Dirichlet problem when ∂Ω = ΓD ; the pure Neumann problem when ∂Ω = ΓN ; and the mixed boundary condition problem when ΓD includes the part of ∂Ω coinciding with one of either θ = 0 or θ = ω, and ΓN = ∂Ω\ΓD with ΓN 6= ∅. The following regularity results can be found in [23] and [21]. Lemma 4.7 Assume |1 − β| < π/ω for the pure Dirichlet problem, 0 < |1 − β| < π/ω for the pure Neumann problem and |1 − β| < π/2ω for the mixed boundary condition problem. Then, for every f ∈ Hβ0 (Ω), there exists a unique solution to (4.15), p ∈ Hβ2 (Ω) for the pure Dirichlet and mixed boundary condition cases and p ∈ Hβ2 (Ω)/R for the pure Neumann problem. Moreover, there exists a constant, C, independent of p, such that kpk2,β ≤ Ckf k0,β .
(4.16)
Proof. See Chapter 1 of [23] for the Dirichlet and Neumann problems and Chapter 2 of [21] for the mixed boundary problem. Define the subspace of functions in Hβ1 (Ω) satisfying the appropriate boundary conditions by Vβ = {v ∈ Hβ1 (Ω) : τ · v = 0 on ΓD , n · v = 0 on ΓN }.
14
LEE, MANTEUFFEL, AND WESTPHAL
We now prove a regularity result for functions in Vβ . Recall that the power of the singularity is defined as α = π/ω for Dirichlet or Neumann boundary conditions and α = π/(2ω) for mixed boundary conditions. Lemma 4.8 Consider domain Ω = Ωw . Then there is a positive constant, C, independent of u, such that, for |1 − β| < α, the following bound holds for all u ∈ V β : kuk1,β ≤ CkLuk0,β .
(4.17)
Proof. From lemma 4.7 we know that any u ∈ Vβ has the decomposition u = ∇φ + ∇⊥ ψ,
(4.18)
where φ, ψ ∈ Hβ2 (Ω) satisfy (
∆φ = ∇ · u, φ = 0,
in Ω, on ∂Ω,
(4.19)
and (
∆ψ = ∇ × u, n · ∇ψ = 0,
in Ω, on ∂Ω.
(4.20)
Then, by applying lemma 4.7 to problems (4.19) and (4.20) we have kuk1,β = k∇φ + ∇⊥ ψk1,β ≤ k∇φk1,β + k∇⊥ ψk1,β
≤ kφk2,β + kψk2,β ≤ C(k∇ · uk0,β + k∇ × uk0,β ) ≤ CkLuk0,β ,
which completes the proof. 5. Error Bounds. Let T h = ∪N i=1 τi be a quasi-uniform triangulation of polygonal domain Ω. Let I h represent standard interpolation onto a piecewise polynomial finite element space of degree k. From finite element theory, we have the following interpolation bounds. Lemma 5.1 Let Ω be a polygonal domain. There exists a constant, C, independent of v, such that, for all v ∈ H m (Ω), 1/2 X kv − I h vk2s,τ ≤ Chm−s |v|m (5.1) τ ∈T h
for 0 ≤ s ≤ m and 1 < m. Here, I h denotes interpolation by a piecewise polynomial of degree k = m − 1. (Note: here the norm k · ks,τ is the standard H s (τ ) norm.) Proof. See [9] or [7]. We now consider a weighted interpolation bound for functions on domains with a polygonal corner at the origin. Define the modified interpolation operator, I 0h , by h Pn I u = i=0 u(ai )φi , if τ does not intersect the origin, P I0h u|τ = n if τ intersects the origin, i=1 u(ai )φi ,
where I h is a standard polynomial interpolation operator, φi are basis functions corresponding to the n + 1 nodal points, ai , and a0 is the origin, (0, 0). Thus, the modified interpolation has a value of zero at the origin and resembles I h away from the origin.
15
FOSLS FOR PROBLEMS WITH SINGULARITIES
Lemma 5.2 Let Ω be a polygonal domain. There exists a constant, C, independent of u, such that, for all u ∈ Hβm (Ω) satisfying equation (4.8),
X
τ ∈T
h
1/2
ku − I0h uk21,β,τ
≤ Chm−1 kukm,β ,
(5.2)
for 1 < m and β > 0, where I0h is the modified interpolation operator onto piecewise polynomials of degree k = m − 1 defined above. Proof. Define K0 = {τ | τ¯ ∩ (0, 0) 6= ∅} as p the set of elements adjacent √ to the origin. On T h \K0 , we have h ≤ rmin ≤ r = x2 + y 2 ≤ rmax ≤ rmin + 2h with rmin = inf{r|(x, y) ∈ τ } and rmax = sup{r|(x, y) ∈ τ } in τ , and ku − I0h uk21,β,τ = ku − I h uk21,β,τ Z = r2β |∇(u − I h u)|2 + r2(β−1) |u − I h u|2 dτ τ Z Z −2 2β 2β rmin |u − I h u|2 dτ ≤ rmax |∇(u − I h u)|2 dτ + rmax ≤
=
τ τ −2 2β 2(m−1) 2 2β 2m Crmax h |u|m,0,τ + Crmax rmin h |u|2m,0,τ −2 2β 2β h2(m−1) |u|2m,0,τ h2(m−1) (1 + rmin Crmax h2 )|u|2m,0,τ ≤ Crmax
−2β 2β rmin ≤ Ch2(m−1) rmax
τ
√
r2β |Dm u|2 dτ
rmin + 2h rmin
≤ Ch2(m−1) ≤ Ch2(m−1)
Z
Z
τ
!2β Z
τ
r2β |Dm u|2 dτ
r2β |Dm u|2 dτ.
We now consider the case for which τ ∈ K0 . Let δ ∈ C ∞ be a cut-off function defined by
δ(r) =
1, 0,
if if
r ≤ h/3, r > 2h/3,
with |δ (m) | ≤ ch−m , where δ (m) is the mth derivative of δ. By the triangle inequality, ku − I0h uk1,β,τ ≤ kδu − I0h (δu)k1,β,τ + k(1 − δ)u − I0h ((1 − δ)u)k1,β,τ .
(5.3)
By the definition of δ we have I0h ((1 − δ)u) = I h ((1 − δ)u) and I0h (δu) = 0. For the second term in (5.3), we apply lemmas 4.4 and 5.1, the properties in δ, and Fubini’s
16
LEE, MANTEUFFEL, AND WESTPHAL
theorem to obtain k(1 − δ)u − I0h ((1 − δ)u)k21,β,τ = k(1 − δ)u − I h ((1 − δ)u)k21,β,τ Z ≤ Ck∇((1 − δ)u − I h ((1 − δ)u))k20,β,τ ≤ Ch2β |∇((1 − δ)u − I h ((1 − δ)u))|2 dτ τ Z X Z m 2 2(β+m−1) 2(β+m−1) |D ((1 − δ)u)| dτ ≤ Ch ≤ Ch |Dm−j (1 − δ)D j u|2 dτ τ |j|≤m
τ
≤ Ch2(β+m−1)
≤ Ch2(β+m−1)
ZZ
h 3
X
|j|≤m−1
X
|j|≤m−1
X Z
= Ch2(m−1)
2h 3
|j|≤m
τ
ZZ
|h|j|−m Dj u|2 dτ +
ZZ
h 3
2h 3 h 3
r(θ)
|(1 − δ)D m u|2 dτ
ZZ −2β 2(β+|j|−m) j 2 h r |D u| dτ +
r(θ)
h 3
(r/h)2β |Dm u|2 dτ
r2(β+|j|−m) |Dj u|2 dτ = Ch2(m−1) kuk2m,β,τ .
Using the properties of δ and Fubini’s theorem result in a similar bound for the first term in (5.3): I0h (δu)k21,β,τ
Z
kδuk21,β,τ
= = kδu − r2β |∇(δu)|2 + r2(β−1) |δu|2 dτ τ Z ≤ C r2β (|∇δ · u|2 + |δ∇u|2 ) + r 2(β−1) |δu|2 dτ τ
≤C
ZZ
≤C
ZZ
≤C
Z
2h 3
h 3 2h 3 h 3
r
r2β h−2 |u|2 dτ + C
ZZ
r2(β−1) |u|2 dτ + C
ZZ
2(m−1)
(r
τ
Thus we have X
τ ∈Th
2(β−m+1)
2h 3
0 2h 3
0
r2β |∇u|2 + r2(β−1) |u|2 dτ r2β |∇u|2 + r2(β−1) |u|2 dτ
2
|∇u| + r2(β−m) |u|2 )dτ ≤ Ch2(m−1) kuk2m,β,τ .
ku − I0h uk21,β,τ ≤ Ch2(m−1)
X
τ ∈Th
kuk2m,β,τ ≤ Ch2(m−1) kuk2m,β ,
and the lemma follows. Lemma 5.3 Assume equation (4.14) holds in Ω. Then, for all uh ∈ V h , kuh k0,β ≤ Ch−η kuh k0,β+η for β > −1 and η > 0. Proof. Using lemma 4.5 and an inverse inequality, we may write kuh k0,β ≤ Ck∇uh k0,β+1 ≤ Ch−1 kuh k0,β+1 ,
(5.4)
FOSLS FOR PROBLEMS WITH SINGULARITIES
17
which establishes (5.4) for η = 1. Repeated application of this inequality thus validates (5.4) for any positive integer. Now consider 1/
1/
kuh k20,β = hr β uh , rβ uh i = hr β− 2 uh , rβ+ 2 uh i
≤ kuh k0,β−1/2 kuh k0,β+1/2 ≤ Ch−1 kuh k20,β+1/2 .
Taking the square root establishes (5.4) for η = 1/2. Repeating these steps leads to (5.4) for all η = ηn = kn /2`n for any nonnegative integers kn and `n . For any η > 0, choose a monotonically decreasing sequence, {ηn }, such that ηn > η and limn→∞ |ηn − η| = 0. Now, gn = (r β+ηn uh )2 is a monotonically increasing function that converges to g = (r β+η uh )2 pointwise everywhere. Thus, by the Lebesgue monotone convergence theorem, we have Z Z h 2 ku k0,β+η = gdx = lim gn dx = lim kuh k20,β+ηn n
n
and, therefore, kuh k0,β = lim kuh k0,β n
≤ lim Ch−ηn kuh k0,β+ηn n
= (lim Ch−ηn )(lim kuh k0,β+ηn ) n
= Ch
n
−η
h
ku k0,β+η ,
which completes the proof. Define an irregular boundary point of polygonal domain Ω to be a point on ∂Ω where interior angle ω satisfies ω > π when Dirichlet or Neumann boundary conditions are applied on both sides of the point or ω > π/2 when one Dirichlet boundary and a Neumann boundary meet at the corner. We now present error bounds for the numerical solution in weighted and unweighted norms. Theorem 5.4 Let Ω be a polygonal domain with one irregular boundary point of interior angle ω and let f ∈ L2 (Ω). Suppose u ∈ V satisfy Lu = f . If uh ∈ V h is chosen to minimize the weighted functional, Gw (uh ; f ) = kLuh − f k20,β = inf kLvh − f k20,β , vh ∈V h
for |1 − β| < α, then the approximation error, u − uh , satisfies the following bounds: ku − uh k1,β ≤ Chα+β−1 kukα+β,β , Gw (u − uh ; 0)
1/2
≤ Chα+β−1 kukα+β,β ,
(5.5) (5.6)
ku − uh k0,β ≤ Chs+β kukα+β,β ,
(5.7)
ku − uh k0 ≤ Chs kukα+β,β ,
(5.8)
where s < α, for α + β ≤ k + 1 with k the degree of the piecewise polynomial elements in V h .
18
LEE, MANTEUFFEL, AND WESTPHAL
Proof. By lemmas 4.8 and 5.2, we have ku − uh k1,β ≤ CkL(u − uh )k0,β ≤ CkL(u − I0h u)k0,β ≤ Cku − I0h uk1,β ≤ Chα+β−1 kukα+β,β ,
which establishes both (5.5) and (5.6) since we may write 1/
kL(u − uh )k0,β = Gw (u − uh ; 0) 2 . Note that lemmas 4.8 and 5.2 are satisfied for |1 − β| < α and α + β ≤ 2. For the weighted L2 -norm, we write X ku − uh k20,β = ku − uh k20,β,τ . τ ∈T h
The Cauchy inequality yields, for any ∈ (0, 1), 1− Z Z Z 1dτ r2β |u − uh |2 dτ ≤ ku − uh k20,β,τ = τ
τ
≤ Ch2−2
Z
τ
r2β |u − uh |2
dτ
2 ·2 β h = Ch2−2 krβ (u − uh )k2 2 r |u − u | dτ 2
τ
1
L (τ )
.
Since H 1 (Ω) is continuously imbedded into Lq (Ω) for all q ∈ [1, ∞), X X ku − uh k20,β = ku − uh k20,β,τ ≤ Ch2−2 krβ (u − uh )k2 2 τ ∈T h 2−2
≤ Ch
krβ (u − uh )k2 2
L (Ω)
τ ∈T h 2−2
≤ Ch
L (τ )
krβ (u − uh )k2H 1 (Ω)
≤ Ch2−2 ku − uh k21,β . Thus, by (5.5), we have ku − uh k0,β ≤ Ch1− ku − uh k1,β ≤ Chs+β kukα+β,β , where any s < α. We now consider the bound on ku − uh k0 . Let K0 = {τ | τ¯ ∩ (0, 0) 6= ∅} and K1 = T h \K0 . First, we consider the case β < 1. If τ ∈ K0 , then r ≤ Ch and r1−β ≤ Ch1−β . Thus, we have ku − uh k20,τ ≤ Ch2(1−β) ku − uh k20,β−1,τ ≤ Ch2(1−β) ku − uh k21,β,τ .
(5.9)
If τ ∈ K1 , we use the technique above to get ku − uh k20,τ ≤ Ch2(1−) ku − uh k2 2
L (τ )
.
Again, since H 1 is continuously imbedded into Lq for all q ∈ [1, ∞), we have X X ku − uh k2 2 ku − uh k20,τ ≤ Ch2−2 τ ∈K1
τ ∈K1
L (τ )
≤ Ch2−2 ku − uh k2 2 ≤ Ch2−2 ku − uh k2H 1 (K1 ) L (K1 ) Z 2−2 r−2β r2β (|u − uh |2 + |∇(u − uh )|2 ) dΩ = Ch ≤ Ch
K1 2(1−β−)
ku − uh k21,β,K1 .
(5.10)
19
FOSLS FOR PROBLEMS WITH SINGULARITIES
Hence by (5.9), (5.10), and (5.5) we have ku − uh k20 =
X
τ ∈K0
ku − uh k20,τ +
≤ Ch2s kuk2α+β,β ,
X
τ ∈K1
ku − uh k20,τ ≤ h2(1−β−) ku − uh k21,β
where s < α. The proof for β ≥ 1 follows analogously. Remark 5.5 For the optimal finite element convergence of O(h) with respect to the weighted functional and H 1 norms, we select β = 2 − α. But theorem 5.4 also requires that β < 1 + α. Thus, when α ∈ [1/2, 1), we may use a weighting with β = 2 − α and expect optimal rates, but when α ∈ (0, 1/2), our theory only guarantees at best O(2α) convergence using β = 1 + α. Numerical results, however, indicate that values of β larger than the theory allows can be used to recover optimal rates. We explore this in the next section. 6. Computational results. In this section, we present some numerical examples of the weighted-norm procedure to validate the error bounds in the previous section. As a test problem, we minimize the weighted functional on the following L-shaped domain: Ω = (−0.5, 0.5)2 \ [0, 0.5) × (−0.5, 0], which yields α = π/ω = 2/3. Function 2 f is chosen so that the solution of this test problem is u = ∇(χ(r)r 3 sin(2θ/3)), where χ(r) = 1 for r < 1/8, χ(r) = 0 for r > 3/8, and χ(r) is C 2 smooth. Again, note that f ∈ L2 (Ω) but u ∈ / H 1 (Ω). Define the following measures of the accuracy of the computed solution, u h : nonweighted functional norm nonweighted L2 norm of the error nonweighted H 1 seminorm of the error weighted functional norm weighted L2 norm of the error weighted H 1 seminorm of the error
G
1/2
1/
= (k∇ · uh − f k20 + k∇ × uh k20 ) 2 ,
0 = ku − uh k0 ,
1 = |u − uh |1 ,
1/
1/
Gw2 = Gw (uh ; f ) 2 , 0w = ku − uh k0,β ,
1w = |u − uh |1,β .
Since α = 2/3, we choose the optimal weight parameter, β = 2 − α = 4/3, for our computations. Table 6.1 summarizes discretization error and convergence rates for β = 4/3. 1 Asymptotic convergence rates in Ω are found to be approximately O(h) for G w/2 2 and 1w , O(h2 ) for 0w and O(h 3 ) for 0 . The approximation does not converge in 1 / H 1 (Ω). either the 1 or G /2 measures since u ∈ To distinguish between behavior near to and away from the singularity, we consider the error of the solution above on a partitioning of Ω. Define Ω0 = Ω ∩ ( 83 , 58 )2 and Ω1 = Ω\Ω0 ; see figure 6.1. Table 6.2 summarizes the asymptotic discretization accuracy obtained at the finest mesh size in subdomains Ω0 and Ω1 . Away from the singularity we observe optimal accuracy in all measures. As expected, near the singularity, the solution fails to converge in the nonweighted functional and H 1 norms. The nonweighted L2 error 2 achieves accuracy of approximately O(h 3 ) near the singularity.
20
LEE, MANTEUFFEL, AND WESTPHAL 1
h 8−1 16−1 32−1 64−1 128−1 256−1 512−1
Gw/2 5.52 4.34 2.34 1.19 5.98e-01 3.00e-01 1.50e-01
Ratio
Rate
1.27 1.85 1.97 1.99 1.99 2.00
0.35 0.89 0.98 0.99 0.99 1.00
h 8−1 16−1 32−1 64−1 128−1 256−1 512−1
0w 3.08e-01 1.35e-01 4.07e-02 1.11e-02 2.98e-03 7.84e-04 2.06e-04
Ratio
Rate
2.28 3.32 3.67 3.72 3.80 3.81
1.19 1.73 1.88 1.90 1.93 1.93
1w 3.81 1.47 6.66e-01 2.97e-01 1.41e-01 6.74e-02 3.31e-02
Ratio
Rate
2.59 2.21 2.24 2.11 2.09 2.04
1.37 1.14 1.16 1.08 1.06 1.03
0 3.72e-01 1.93e-01 8.93e-02 4.93e-02 3.00e-02 1.87e-02 1.18e-02
Ratio
Rate
1.93 2.16 1.81 1.64 1.60 1.58
0.95 1.11 0.86 0.71 0.68 0.66
Table 6.1 Convergence of discretization error for weighted-norm FOSLS.
Ω1
Ω0
Fig. 6.1. L-shaped domain Ω and subdomains Ω0 and Ω1 .
Gw/2
1
G /2
1
1w
1
0w
0
Ω1
O(h)
O(h)
O(h)
O(h)
O(h2 )
O(h2 )
Ω0
O(h)
O(1)
O(h)
O(1)
O(h2 )
O(h 3 )
Ω
O(h)
O(1)
O(h)
O(1)
O(h2 )
O(h 3 )
2 2
Table 6.2 Accuracy in Ω0 , Ω1 , and Ω with β = 2 − α.
Figure 6.2 shows the first component of the exact solution, u1 , and the standard FOSLS approximation uh1 . Figure 6.3 shows the error of the first component of the approximated solution for the standard FOSLS and the weighted-norm FOSLS methods. We see that the error in the approximation in standard FOSLS is highest near the
21
FOSLS FOR PROBLEMS WITH SINGULARITIES
2
2
1
1 1
0
0.4
0.6
−2 0
0.4
0.2
0.8
−1
0.6
−2 0
1
0
0.8
−1
0.4
0.2
0.6 0.8 1
0.4
0.2 0.2
0.6 0.8
0
1
0
(b) Component uh 1 computed by standard FOSLS.
(a) Exact component u1 .
Fig. 6.2. Exact solution component u1 and solution component uh 1 approximated by standard FOSLS on the h = 32−1 mesh.
L2 error in u1
1
L2 error in u1
1
0.9
0.9 2.5
0.8
2.5 0.8
0.7
2
0.6
0.7
2
0.6
0.5
1.5
0.4
0.5
1.5
0.4 1
0.3 0.2
1
0.3 0.2
0.5 0.1
PSfrag replacements
0
0.5 0.1
0
0.2
0.4
PSfrag replacements 0.6 0.8 1
(a) Standard FOSLS.
0
0
0.2
0.4
0.6
0.8
1
(b) Weighted-norm FOSLS, β = 2 − α.
Fig. 6.3. Reduction of the pollution effect by the weighted-norm procedure. Each plot is the −1 mesh. error of solution component uh 1 on the h = 32
singularity, but remains large even away from the corner point. In the weighted-norm FOSLS implementation, the error remains large near the singularity, as we expect, but is now concentrated only near the corner point. The pollution effect is removed by the weighting procedure. There are many boundary value problems not directly covered by the theory presented here that are of interest. For example, Poisson’s equation with mixed boundary conditions on the domain used above has a value of α = 1/3. To recover optimal convergence for this problem, the weighted-norm method requires a value of β larger than theorem 5.4 allows. In other elliptic equations (e.g., Stokes or the linear elasticity equations), the value of α is generally smaller than for Poisson’s equation for the same domain and boundary condition type. In each of these cases, a larger β
22
LEE, MANTEUFFEL, AND WESTPHAL
value is necessary for optimal convergence. This leads us to consider using larger β than the theory allows. Consider the same example problem as above on uniform mesh sizes of h = 1/8, 1/16, ..., 1/512, and values of β ranging from 1/3 to 23/6. Figure 6.4 plots the convergence rate at the finest level for: the weighted functional 1 norm, Gw/2 ; the weighted L2 norm, 0w ; and the L2 norm, 0 . While the functional Convergence Rates vs. Β 2 weighted L2 error
weighted functional norm 1 0.67 L2 error
0 0
0.33
1
1.66
Β
3
4
Fig. 6.4. Convergence rates versus β. The shaded region indicates values of β for which the assumptions of theorem 5.4 are satisfied.
norm retains optimal accuracy for large values of β, the solution fails to converge in the weighted and nonweighted L2 measures for β & 3. This indicates that, although the weighted-norm approach seems to be more robust than the theory allows, large values of β should still be used with caution. The method presented here is applicable to a wide range of problems and provides an efficient alternative to more specialized techniques for treating singularities in boundary value problems. Further numerical results for other problems including systems and in three dimensions can be seen in a companion paper. REFERENCES [1] T. Apel and B. Heinrich, Mesh refinement and windowing near edges for some elliptic problem, SIAM J. Numer. Anal., 31 (1994), pp. 695–708. [2] I. Babuska, R. Kellogg, and J. Pitkaranta, Direct and inverse error estimates for finite elements with mesh refinements, Numer. Math., 33 (1979), pp. 447–471. [3] M. Berndt, T. Manteuffel, S. McCormick, and G. Starke, Analysis of first-order system least squares (fosls) for elliptic problems with discontinuous coefficients: part i, SIAM J. Numer. Anal., 43 (2005), pp. 386–408. [4] , Analysis of first-order system least squares (fosls) for elliptic problems with discontinuous coefficients: part ii, SIAM J. Numer. Anal., 43 (2005), pp. 409–436. [5] P. Bochev and M. Gunzburger, Analysis of least-squares finite element methods for the stokes equations, Math. Comp., 63 (1994), pp. 479–506. [6] , On least-squares finite element methods for the poisson equation and their connection to the dirichlet and kelvin principles, SIAM J. Numer. Anal., 43 (2005), pp. 340–362. [7] D. Braess, Finite Elements: Theory, Fast Solvers and Applications in Solid Mechanics, Cambridge, 2001. [8] J. Bramble, R. Lazarov, and J. Pasciak, A least-squares approach based on a discrete minus one inner product for first order systems, Math. Comp., 66 (1997), pp. 935–955.
FOSLS FOR PROBLEMS WITH SINGULARITIES
23
[9] S. Brenner and L. Scott, The Mathematical Theory of Finite Element Methods, SpringerVerlag, 1994. [10] Z. Cai and S. Kim, A finite element method using singular functions for the poisson equation: Corner singularities, SIAM J. Numer. Anal., 39 (2001), pp. 286–299. [11] Z. Cai, T. Manteuffel, and S. McCormick, First-order system least squares for the stokes equations, with application to linear elasticity, SIAM J. Numer. Anal., 34 (1997), pp. 1727– 1741. , First-order system least squares for velocity-vorticity-pressure form of the stokes [12] equations, with application to linear elasticity, Elect. Trans. Numer. Anal., 3 (1997), pp. 150–159. [13] Z. Cai, T. Manteuffel, S. McCormick, and S. Parter, First-order system least squares (FOSLS) for planar linear elasticity: Pure traction problem, SIAM J. Numer. Anal., 35 (1998), pp. 320–335. [14] Z. Cai, T. Manteuffel, S. McCormick, and J. Ruge, First-order system LL ∗ (FOSLL*): Scalar elliptic partial differential equations, SIAM J. Numer. Anal., 39 (2002), pp. 1418– 1445. [15] C. Cox and G. Fix, On the accuracy of least squares methods in the presence of corner singularities, Comp. and Maths. with Appls., 10 (1984), pp. 463–475. [16] G. Fix, S. Gulati, and G. Wakoff, On the use of singular functions with finite element approximations, J. Comput. Phys., 13 (1973), pp. 209–228. [17] P. Grisvard, Elliptic Problems in Nonsmooth Domains, Pitman, 1985. , Singularities in Boundary Value Problems, Springer-Verlag, 1992. [18] [19] S. Kim, T. Manteuffel, and S. McCormick, First-order system least squares (fosls) for spatial linear elasticity: Pure traction, SIAM J. Numer. Anal., 38 (2001), pp. 1454–1482. [20] V. Kondratiev, Boundary problems for elliptic equations in domains with conical or angular points, Trans. Moscow Math. Soc., 16 (1967), pp. 227–313. [21] V. Kozlov, V. Maz’ya, and J. Rossmann, Spectral Problems Associated with Corner Singularities of Solutions to Elliptic Equations, vol. 85, American Mathematical Society, 2001. [22] T. Manteuffel, S. McCormick, J. Ruge, and J. Schmidt, First-order system ll* (fosll*) for general scalar elliptic problems in the plane, SIAM J. Numer. Anal., 43 (2005), pp. 2098– 2120. [23] V. Maz’ya, S. Nazarov, and B. Plamenevskij, Asymptotic Theory of Elliptic Boundary Vaule Problems in Singularly Perturbed Domains, Volume I, vol. 111, Birkh¨ auser, 2000. Operator Theory Advances and Applications. [24] G. Strang and G. Fix, An analysis of the finite element method, Prentice-Hall, Englewood Cliffs, NJ, 1973.
1. Introduction. In this paper, we develop a method for treating div/curl systems with reduced regularity. While motivated by first-order systems that arise from second-order elliptic boundary value problems, div/curl systems appear in many contexts, for example, in Maxwell’s equations and in the vorticity from of Stokes equations. These problems have the fortunate property of a guaranteed smooth solution as long as the data and domain are smooth. However, many problems of interest are posed in non-smooth domains and, as a consequence, lose this property at a finite number of points on the boundary in two dimensions or along curves on the boundary in three dimensions. In the present paper, we study two-dimensional problems that have non-smooth solutions at irregular boundary points, that is, points that are corners of polygonal domains, locations of changing boundary condition type, or both. Similar behavior occurs in problems with discontinuous material coefficients and the methods presented here can easily be extended to that situation. Standard solution techniques that attempt to approximately solve a div/curl system with reduced regularity using H 1 -conforming finite elements will, in general, fail to converge. This phenomenon can be explained by noting that the Sobolev space (H 1 )2 is a closed subspace of H(div) ∩ H(curl) (see section 2), which implies either that (H 1 )2 = H(div) ∩ H(curl), the case of full regularity, or (H 1 )2 is a proper subspace. In this case, the co-dimension is finite and is spanned by so-called singular functions. In the presence of reduced regularity, the solution will, in general, not be in (H 1 )2 . A standard finite element method using H 1 -conforming elements will converge to the element of (H 1 )2 closest to the true solution in the H(div) ∩ H(curl) norm. Local mesh refinement will not alter this outcome. If a basis for the singular functions is known, it can be incorporated directly into the finite element space (c.f.[10, 16, 24]). In [3, 4], this approach is applied to div/curl systems and shown to restore optimal convergence throughout the domain at a minimal additional cost. For some two-dimensional problems, the singular basis functions are known and can be included in the finite element space. For the other ∗ Department of Applied Mathematics, University of Colorado at Boulder; Boulder,CO 80309– 0526, USA. Email: {tmanteuf,eunjung}@colorado.edu. † Department of Mathematics and Computer Science, Wabash College, PO Box 352, Crawfordsville, IN 47933, USA. Email: [email protected]. ‡ This work was sponsored by the National Science Foundation under grant number DMS– 0420873, and the Department of Energy under grant numbers DE–FC02-01ER25479 and DE–FG0203ER25574. § This work was sponsored by the National Science Foundation under grant DMS–9810751
1
2
LEE, MANTEUFFEL, AND WESTPHAL
problems, or in three dimensions, the exact character of the singular functions is less well understood which makes this approach more difficult to implement. As a different type of remedy for this so-called pollution effect, least-squares methods based on inverse norms can be effective for problems with irregular boundary points, discontinuous coefficients, and data in H −1 . For example, in [5, 8, 11–13, 19], the functional is posed in terms of H −1 -norms rather than L2 -norms, resulting in optimal L2 -approximations to the solution. A more recent approach, called FOSLL∗ , uses an inverse norm induced by the equations, and is shown in [14, 22] to be more efficient than the H −1 -norm methods. Graded mesh refinement in weighted Sobolev spaces has been shown to be effective in restoring optimal convergence, in L2 - and H 1 -norms, in the context of a Galerkin formulation of second-order elliptic problems with reduced regularity (c.f. [2], [1]). However, in that context, the solution is in H 1 . Convergence would occur, although more slowly, without weighting and mesh refinement. In [15], a weighted norm and a sequence of graded meshes is used in an H(div) least-squares functional, arising from a second-order elliptic problem, to restore optimal convergence for the primal variable, in both L2 - and H 1 -norms, if the flux variable is approximated in a finite element space satisfying the grid decomposition property, for example, Raviart-Thomas elements (c.f.[6]). The flux variable converges in a weighted H(div)-norm. In this paper we examine div/curl systems that lack regularity within the firstorder system least-squares (FOSLS) framework. The basic FOSLS approach is to recast the original system as an appropriate first-order system and apply an L 2 minimization principle over the residual of the equations. If possible, this reformulation is done by minimizing a functional whose quadratic part is equivalent to the product H 1 -norm, indicating that the process is similar to solving a weakly coupled system of Poisson-like equations. This equivalence also guarantees optimal H 1 -accuracy for standard discretizations. For div/curl systems with reduced regularity, as briefly mentioned above, the L2 based functional fails to be H 1 -equivalent and, as a consequence, standard discretizations suffer from the pollution effect. Here, a weighted-norm leastsquares method is developed that restores optimal convergence using H 1 -conforming finite elements without graded mesh refinement. It replaces the L2 -norms in the FOSLS functional with weighted L2 -norms, making the functional norm equivalent to a weighted H 1 -norm. With an appropriate weighting function, this method recovers optimal order accuracy in the weighted L2 - and H 1 -norms, and retains optimal L2 convergence even near the singularity. Our method requires only the power of the singularity (not the actual singular solution) to be known a priori and, in practice, can be used with only a rough estimate of the power of the singularity, which can be adaptively determined if unknown (c.f.[4]). The method developed here has some similarity to [15], but considers an H(div)∩ H(curl) functional and considers more general boundary conditions. Most importantly, our analysis admits a more aggressive weighting, resulting in optimal order accuracy in the weighted norms without mesh refinement. In addition, we prove several Poincar´e type inequalities in weighted Sobolev spaces under a variety of boundary conditions that, in addition to being necessary for our main result, may be of independent interest. We use Poisson’s equation on a domain with a reentrant corner as a model problem and as the formal setting for analysis. The resulting div/curl system is the focus. The analysis in this paper is restricted to two dimensions. However, the approach suggests
FOSLS FOR PROBLEMS WITH SINGULARITIES
3
a natural generalization to problems with reduced regularity due to discontinuous coefficients and to problems in three dimensions. This paper is organized as follows. Section 2 contains notation and preliminary discussion. In section 3, the weighted FOSLS functionals are described. Section 4 contains Poincar´e inequalities and regularity results in weighted Sobolev spaces. Error bounds for the weighted FOSLS functional are presented in section 5 and section 6 contains computational results. 2. Singular solutions and preliminaries. For vector function u = (u 1 , u2 )t , let the divergence and curl of u be defined in the standard way: ∇ · u = ∂x u1 + ∂y u2 and ∇ × u = ∂x u2 − ∂y u1 . Further, define the formal adjoint of the curl operator by ∂y q . ∇⊥ q = −∂x q We use standard notation for Sobolev spaces H k (Ω)d , corresponding inner product (·, ·)k,Ω , and norm k · kk,Ω , for k ≥ 0. We drop subscript Ω and superscript d when the domain and dimension are clear by context. Since H 0 (Ω) coincides with L2 (Ω) we often denote k · k0 by k · k. Define the subspaces of L2 (Ω) induced by the divergence and curl of u by H(div) = {u ∈ L2 (Ω) : k∇ · uk < ∞},
H(curl) = {u ∈ L2 (Ω) : k∇ × uk < ∞}. We also make use of the following general inequalities for nonnegative a and b: |a|2 + |b|2 ≤ |a + b|2 ≤ 2(|a|2 + |b|2 ).
(2.1)
Consider the function f (r, θ) = r a in two dimensional polar coordinates. Assume that the origin lies on the boundary of domain Ωw = {(r, θ) : 0 < r < R, 0 < θ < ω < 2π}, as pictured in figure 2.1. By a direct computation it is clear that f ∈ H k (Ω) only for k < a + 1.
ω R
Fig. 2.1. Simple wedge-shaped domain, Ωw .
Now, consider Poisson’s equation on a domain in R2 with a corner of interior angle ω. It is well known that, for the case of Dirichlet or Neumann boundary conditions, the solutions of this boundary value problem may include those with radial part of
4
LEE, MANTEUFFEL, AND WESTPHAL π
the form p ∼ r ω in a local polar coordinate system centered at the corner. Thus, for the case of reentrant corners, ω > π, the solution fails to be in H 2 (Ω) and we say that the problem has a singularity (or singular solution). For problems with Dirichlet and Neumann boundaryπ conditions meeting at the corner, solutions may have components of the form p ∼ r (2ω) . Thus, for mixed boundary conditions, singularities may occur at corners with ω > π/2. We now explore this issue in more detail. Define the power of the singularity to be α = π/ω for Dirichlet or Neumann boundary conditions and α = π/(2ω) for mixed boundary conditions. The solution to Poisson’s equation may be written as p(r, θ) = p0 (r, θ) + s(r, θ), where p0 (r, θ) ∈ H 2 (Ω) and s(r, θ) ∈ H 1+m (Ω) for m < α. The singular part of the solution has the form s(r, θ) = r α (κ1 sin(αθ) + κ2 cos(αθ)), where the values of κ1 and κ2 depend on boundary conditions (see [17, 18]). For the FOSLS formulation of this problem, we may similarly decompose unknown u = ∇p as u(r, θ) = u0 (r, θ) + ∇s(r, θ), where ∇s(r, θ) has the form ∇s(r, θ) = αr α−1
κ1 sin(α − 1)θ + κ2 cos(α − 1)θ . κ1 cos(α − 1)θ − κ2 sin(α − 1)θ
Thus, the unknown u(r, θ) is in H k (Ω) only for k < α. For example, consider Poisson’s equation posed on the simple domain in figure 2.1. Let the solution to this boundary value problem in polar coordinates be 2 p = χ(r)r 3 sin(2θ/3), where χ(r) is a smooth transition function that is 1 on a platform near the origin and vanishes at the boundaries not adjacent to the origin. Then, p = 0 on ∂Ω and ∆p =
1 1 2 ∂r (r∂r p) + 2 ∂θθ p r r 2
= r 3 χ00 (r) + 73 r
−1 3
χ0 (r) sin(2θ/3),
and, thus, it is clear that ∆p ∈ L2 (Ω), but p ∈ / H 2 (Ω). We say this problem fails to 2 provide full lifting of the data (from L (Ω) to H 2 (Ω) e.g.). The solution, u = ∇p, is, thus, not in H 1 (Ω). 3. Weighted-norm least squares. As before, let Ω be a domain with a corner of interior angle ω at the origin, and we may, without loss of generality, further assume diam(Ω) ≤ 1. For f ∈ L2 (Ω), let p satisfy in Ω, −∆p = f, p = 0, on ΓD , (3.1) n · ∇p = 0, on ΓN ,
FOSLS FOR PROBLEMS WITH SINGULARITIES
5
¯D ∪ Γ ¯ N . When this problem where n is the outward unit normal to Ω and ∂Ω = Γ 2 is H regular, the normal FOSLS methodology is to introduce the new unknown, u = ∇p, and rewrite system (3.1) as u − ∇p = 0, in Ω, −∇ · u = f, in Ω, in Ω, ∇ × u = 0, τ · u = 0, on ΓD , (3.2) n · u = 0, on Γ , N p = 0, on Γ D, n · ∇p = 0, on ΓN . Here, τ is the counter-clockwise unit tangental vector to Ω. Since this system can be posed completely in terms of u, we may decouple the equations in (3.2), solve for u first, and then recover p from u. To this end, define the two L2 -norm functionals, G(u; f ) = k∇ · u + f k2 + k∇ × uk2 ,
G2 (p; u) = ku − ∇pk2 , and the spaces,
V = {v ∈ H 1 (Ω) : τ · v = 0 on ΓD , n · v = 0 on ΓN },
W = {q ∈ H 1 (Ω) : q = 0 on ΓD }.
Thus, the two-stage solution process is to minimize G(v; f ) over V and then, given the minimizer, u, minimize G2 over W: (1)
G(u; f ) = inf G(v; f ),
(2)
G2 (p; u) = inf G2 (q; u).
v∈V
(3.3)
q∈W
The goal of the FOSLS methodology is, generally, to formulate functionals whose quadratic part is equivalent to the H 1 norm whenever possible. The second stage functional is H 1 equivalent and the solution we seek is always in H 1 . The first stage functional, however, is not always H 1 equivalent. For domains with reentrant corners, there is no H 1 sequence of functions that converges to the solution in the H(div) ∩ 2 H(curl) norm. To illustrate, consider the example above where p = χ(r)r 3 sin(2θ/3) and u = ∇p. A simple computation reveals that ∇·u, ∇×u ∈ L2 (Ω), but u ∈ / H 1 (Ω). Define the weighted functional by Gw (u; f ) = kw(∇ · u + f )k2 + kw∇ × uk2 ,
(3.4)
where the weight function has the form w = r β for some β > 0. Define the weighted Sobolev norm, k · kk,β , on Ω in terms of the standard L2 norm, k · k0 , by kqkk,β = (
X
|j|≤k
1/
krβ−k+|j| Dj qk20 ) 2 ,
(3.5)
6
LEE, MANTEUFFEL, AND WESTPHAL
where D j is the standard distributional derivative of order j. Similarly, define the weighted seminorm by |q|k,β = (
X
|j|=k
krβ−k+|j| Dj qk20 )
1/2
(3.6)
and the associated weighted Sobolev space by Hβk (Ω) = {q : kqkk,β < ∞}. Define the div/curl operator, L, and vector f by L = may now write the weighted functional from (3.4) as
(3.7)
∇· ∇×
and f =
f . We 0
Gw (u; f ) = kLu − f k20,β . The weighted-norm least-squares minimization problem for the first-stage solution is then: find u ∈ V such that Gw (u; f ) = inf Gw (v; f ). v∈V
The second-stage solution for p remains as described above. We seek values of β that make H 1 (Ω) dense in H(div) ∩ H(curl) in the weighted functional norm and result in the most accurate discretizations possible. For the discrete problem, we may choose any finite dimensional subset of H 1 over which to minimize the weighted functional. Let P h denote the space of C 0 piecewise polynomial (or tensor product) elements on triangles (or quadrilaterals) of meshsize h and V h the subspace of P h that satisfies the appropriate boundary conditions on Ω: V h = {vh ∈ P h : τ · vh = 0 on ΓD , n · vh = 0 on ΓN }. The discrete weighted-norm least-squares minimization problem is, then, to minimize the discrete functional: find uh ∈ V h such that Gw (uh ; f ) = min Gw (vh ; f ). vh ∈V h
(3.8)
By unweighting the equations near the singularity, the functional is freed from trying to approximate the solution (which is not in H 1 (Ω)) in the H 1 sense near the singularity. But, away from the singularity, the weighted functional retains the same character as the normal non-weighted functional. We now consider the choice of weight parameter β and its relation to weighted and nonweighted a priori error bounds on the approximated solution. 4. Poincar´ e Bounds and Regularity Estimates. In this section, we establish several theoretical results in weighted Sobolev spaces and error bounds for the weighted-norm method. Here, we establish several Poincar´e bounds in the domain Ωw . We first prove a result for the scalar pure Neumann and pure Dirichlet problems, and then for the scalar mixed boundary condition problem. These results lead to a Poincar´e inequality for the vector case.
7
FOSLS FOR PROBLEMS WITH SINGULARITIES
Lemma 4.1 Take Ω = Ωw and let β > 0, > 0 and γ = β − 3/2 − . Further, assume that γ 6= −2. For functions q ∈ Hβ1 (Ω), with r β− ∇q, r γ+1 ∇q ∈ L2 (Ω), that can be chosen to satisfy ZZ rγ q(r, θ) r dr dθ = 0, (4.1) Ω
we have the bound: kqk0,β−1 ≤ C()k∇qk0,β− ,
(4.2)
where C() depends on , β, and Ω and C() → ∞ as → 0. Proof. For any two points (r, θ) and (r0 , θ0 ) in Ω, write q(r, θ) as q(r, θ) = q(r0 , θ0 ) +
Z
r r0
∂q(ˆ r, θ) dˆ r+ ∂ rˆ
Z
θ θ0
ˆ ∂q(r0 , θ) ˆ dθ. ∂ θˆ
Multiplying both sides of the equation by r0γ+1 , integrating with respect to r0 and θ0 over Ω, and Fubini’s theorem yield Rγ+2 ω q(r, θ) γ+2
=
Z ωZ RZ 0
0
r0
∂q(ˆ r, θ) r0γ+1 ∂ rˆ
dˆ r dr0 dθ0 +
Z ωZ RZ 0
0
θ θ0
r0γ+1
ˆ ∂q(r0 , θ) dθˆ dr0 dθ0 ∂ θˆ
Z RZ R ∂q(ˆ r, θ) ∂q(ˆ r, θ) dr0 dˆ r−ω dr0 dˆ r r0γ+1 ∂ rˆ ∂ rˆ 0 r 0 rˆ Z RZ ωZ ω Z RZ θZ θˆ ˆ ˆ ∂q(r0 , θ) γ+1 ∂q(r0 , θ) ˆ dθ0 dθ dr0 − dθ0 dθˆ dr0 + r0 r0γ+1 ∂ θˆ ∂ θˆ 0 θ θˆ 0 0 0
= ω
=
Z
r
r
Z
rˆ
r0γ+1
Z R Z ωRγ+2 R ∂q(ˆ ω ∂q(ˆ r, θ) r, θ) dˆ r− dˆ r rˆγ+2 γ+2 0 ∂ rˆ γ+2 r ∂ rˆ Z RZ ω Z RZ ω ˆ ˆ ∂q(r0 , θ) γ+1 ∂q(r0 , θ) ˆ ˆ + dθ dr0 − ω dθˆ dr0 . θ r0 r0γ+1 ∂ θˆ ∂ θˆ 0 0 0 0
Note that above q is, by assumption, sufficiently smooth to apply Fubini’s theorem. By the triangle inequality we have that γ+2 R ω γ + 2 q(r, θ) Z R Z ω r, θ) ωRγ+2 R ∂q(ˆ r, θ) γ+2 ∂q(ˆ ≤ rˆ r+ dˆ r ∂ rˆ dˆ γ+2 0 γ + 2 r ∂ rˆ Z RZ ω ∂q(r , θ) ˆ 0 + 2ω r0γ+1 dθˆ dr0 . ∂ θˆ 0 0
8
LEE, MANTEUFFEL, AND WESTPHAL
Now, squaring each side and using (a + b + c)2 ≤ 3(a2 + b2 + c2 ) we get, !2 !2 Z R Z R ∂q(ˆ r, θ) r, θ) 3 2 γ+2 ∂q(ˆ r +3 r rˆ |q(r, θ)| ≤ 2(γ+2) ∂ rˆ dˆ ∂ rˆ dˆ R 0 r !2 Z RZ ω ˆ 12(γ + 2)2 γ+1 ∂q(r0 , θ) ˆ + r0 dθ dr0 . ∂ θˆ R2(γ+2) 0 0
(4.3)
Multiply each side of (4.3) by r 2β−1 and integrate with respect to r and θ over Ω. We consider each of the terms on the resulting right-hand side separately. First, !2 Z ωZ R Z R ∂q(ˆ r , θ) dˆ dr dθ r2β−1 rˆγ+2 r ∂ r ˆ 0 0 0 2 Z ωZ RZ R r, θ) 2β−1 2γ+4 ∂q(ˆ r rˆ ≤ R r dr dθ ∂ rˆ dˆ 0 0 0 2 Z Z r, θ) R2β+1 ω R 2γ+3 ∂q(ˆ rˆ = r dθ ∂ rˆ rˆ dˆ 2β 0
0
R2β+1 k∇qk20,γ+ 3 . ≤ 2 2β
We now consider the second term in (4.3). Since by the Schwarz inequality, !2 Z R ∂q(ˆ r, θ) r ∂ rˆ dˆ r !2 Z R r, θ) (−1)/2 (1−)/2 ∂q(ˆ = rˆ rˆ r ∂ rˆ dˆ r 2 Z R Z R ∂q(ˆ r, θ) dˆ r rˆ1− rˆ−1 dˆ r ≤ ∂ rˆ r r 2 Z r, θ) R R 1− ∂q(ˆ rˆ r, = ∂ rˆ dˆ r
we can apply Fubini’s theorem to bound the second term: !2 Z R Z ωZ R ∂q(ˆ r , θ) 2β−1 r r dr dθ ∂ rˆ dˆ r 0 0 2 Z Z Z r, θ) R ω R R 2β−1 1− ∂q(ˆ r dr dθ ≤ r rˆ ∂ rˆ dˆ 0 0 r 2 Z Z Z R ω R rˆ 2β−1 1− ∂q(ˆ r, θ) = r dθ r rˆ ∂ rˆ dr dˆ 0 0 0 2 Z ωZ R ∂q(ˆ r, θ) R = rˆ dˆ r dθ rˆ2β− 2β 0 0 ∂ rˆ R = k∇qk20,β− 2 2β R k∇qk20,β− , ≤ 2β
9
FOSLS FOR PROBLEMS WITH SINGULARITIES
The third term can be bounded similarly: !2 Z RZ ω Z ωZ R ∂q(r , θ) ˆ 0 γ+1 r0 r2β−1 dr dθ dθˆ dr0 ∂ θˆ 0 0 0 0 2 ! Z ωZ R Z RZ ω ∂q(r , θ) ˆ 0 2γ+2 ≤ r0 r2β−1 dr dθ Rω dr0 dθˆ ∂ θˆ 0 0 0 0 2 Z Z ˆ R2β+1 ω 2 ω R 2γ+3 1 ∂q(r0 , θ) r0 = r0 dr0 dθˆ r0 2β ∂ θˆ 0 0 =
R2β+1 ω 2 k∇qk20,γ+ 3 . 2 2β
Putting the three terms together and substituting γ = β − 3/2 − we may now write the bound Z ωZ R 2 r2β−2 |q(r, θ)| r dr dθ 0 0 −2 6(β + 21 − )2 R−2 ω 2 3R 3R ≤ + + k∇qk20,β− , 2β 2β β and the lemma follows by taking the square root of both sides. In what follows, let χ be a smooth function of r where χ = 1 for r < η and χ = 0 for r > 2η. We take η to be sufficiently small to ensure that supp(χ) ⊂ Ω. Lemma 4.2 Take Ω = Ωw and let q be a scalar function in Hβ1 (Ω), where β > 0. The following bound holds for χ as defined above: kχqk0,β−1 ≤
1 k∇(χq)k0,β . β
(4.4)
Proof. Hardy’s inequality for f (t) defined for t > 0 with lim f (t) = 0 gives t→0
(see [20]): Z
∞ 0
f2 dt ≤ 4 t2
Z
∞ 0
2
|f 0 | dt.
(4.5)
The lemma follows after a change of variables, t = r −2β , a substitution f (r) = χq(r, θ) for fixed θ, and an integration on both sides with respect to θ. Lemma 4.3 Take Ω = Ωw and let either q ∈ Hβ1 (Ω) with q = 0 on ∂Ω or q ∈ Hβ1 (Ω)/R with n · ∇q = 0 on ∂Ω. Then kqk0,β−1 ≤ Ck∇qk0,β for β > 0 where C depends only on Ω, and β. Proof. First, if q = 0 on ∂Ω, write q = q(r, θ) as q(r, θ) =
Z
θ 0
ˆ ∂q(r, θ) ˆ dθ. ∂ θˆ
(4.6)
10
LEE, MANTEUFFEL, AND WESTPHAL
Square both sides and multiply by r 2β−1 : r
2β−1
Z 2 θ ∂q(r, θ) ˆ dθˆ |q(r, θ)| = r ˆ 0 ∂θ 2 Z ω ˆ 1 ∂q(r, θ) 2β+1 ˆ ≤r ω dθ. ∂ θˆ 0 r 2
2β−1
Integrate both sides with respect to r and θ over Ω: Z
ω 0
Z
R 0
2 1 ∂q(r, θ) ˆ r2β+1 r2β−2 |q(r, θ)| r dr dθ ≤ ω dθˆ dr dθ ∂ θˆ 0 r 0 0 2 Z ω Z R0 1 ∂q(r, θ) ˆ ˆ r2β ≤ ω2 r dr dθ. r ∂ θˆ 0 0 Z
2
ω
Z
R0
Z
ω
The lemma follows since the right side is bounded by Ck∇qk20,β . Now, if q ∈ Hβ1 (Ω)/R then it may be chosen to satisfy ZZ rγ q r dr dθ = 0
(4.7)
Ω
for γ chosen as in lemma 4.1. By the triangle inequality and lemma 4.2 we get kqk0,β−1 ≤ kχ qk0,β−1 + k(1 − χ)qk0,β−1 ! 12 Z ω Z R 2 r 2β−2 2 2 ≤ Ck∇(χq)k0,β + r (1 − χ) q r dr dθ η 0 η ≤ C (k∇(q)k0,β + kqk0,β ) . Apply lemma 4.1 with = 1 to the kqk0,β term on the right side and the lemma follows. For the problem with mixed boundary conditions, consider Ωw partitioned into subdomains Ω0 = {(r, θ) : r ≤ 21 R0 , 0 ≤ θ ≤ ω} and Ω1 = Ω\Ω0 , as shown in figure 4.1.
Ω1 Ω0 1 R 2 0
R0 R
Fig. 4.1. Wedge-shaped domain, Ωw , partitioned into subdomains Ω0 and Ω1 .
FOSLS FOR PROBLEMS WITH SINGULARITIES
11
Lemma 4.4 Consider domain Ω = Ωw , as pictured in figure 4.1, and let q ∈ Hβ1 (Ω) vanish on the line segment of ∂Ω corresponding to θ = 0 and r < R0 . Then, there is a constant, C, dependent only on Ω, β, and R0 , such that kqk0,β−1 ≤ Ck∇qk0,β
(4.8)
for β > 0. Proof. For points (r, θ) in Ω0 we may derive the bound, kqk0,β−1,Ω0 ≤ Ck∇qk0,β,Ω0 ,
(4.9)
completely analogous to the proof of lemma 4.3. Now, consider points (r, θ) in Ω 1 . We may write q = q(r, θ) as q(r, θ) =
Z
r r˜
∂q(ˆ r, θ) dˆ r+ ∂ rˆ
Z
θ 0
ˆ ∂q(˜ r, θ) ˆ dθ, ∂ θˆ
where the point (˜ r, 0) is on the part of ∂Ω1 where q vanishes. By the Schwarz inequality, the triangle inequality and inequality (2.1) we have the bound 1 |q(r, θ)| ≤ 2(R − R0 ) 2 2
Z
r r˜
2 2 Z θ ˆ ∂q(ˆ r, θ) r, θ) ∂q(˜ ˆ r + 2ω dθ. ∂ rˆ dˆ ∂ θˆ 0
(4.10)
We now expand the limits in the integrals, multiply each side by r 2β−1 , integrate with respect to r over ( 21 R0 , R), integrate with respect to θ over (0, ω), and integrate with respect to r˜ over ( 12 R0 , R0 ), by Fubini’s theorem: 1 ( R0 ) 2
Z
ω 0
Z
R 1 2 R0
r2β−1 |q(r, θ)|2 dr dθ
2 Z R0Z ωZ R Z R 1 r, θ) 2β−1 ∂q(ˆ ≤ 2(R − R0 ) r r dr dθ d˜ r ∂ rˆ dˆ 1 1 1 2 2 R0 0 2 R0 2 R0 2 Z R0Z ωZ R Z θ ˆ r, θ) ∂q(˜ r. + 2ω dθˆ dr dθ d˜ 1 1 ∂ θˆ 0 2 R0 0 2 R0
(4.11)
We use the inequalities,
1 R0 ≤ r˜ ≤ R0 , 2
r˜ ≤ rˆ ≤ r ≤ R,
to derive the following simple bounds: 1≤
2ˆ r R0
,
r≤
2R R0
rˆ,
r≤
2R R0
r˜.
(4.12)
12
LEE, MANTEUFFEL, AND WESTPHAL
By applying the bounds in (4.12) and Fubini’s theorem, we may now write (4.11) as 1 ( R0 ) 2
Z
ω 0
Z
R 1 2 R0
r2β−2 |q(r, θ)|2 r dr dθ
2 2 2β−1 r, θ) 2ˆ r 2R 2β−1 ∂q(ˆ r dθ rˆ ∂ rˆ dˆ 1 R0 R0 0 2 R0 2 2β−1 Z R0 Z ω ∂q(˜ ˆ r , θ) 1 2R + 2ω 2 (R − R0 ) r r˜2β−1 dθˆ d˜ 1 ∂ θˆ 2 R 0 0 R 0 2 2 Z ωZ R r, θ) 1 2β+1 2β ∂q(ˆ 2 −2β 2β−1 rˆ dˆ r dθ ≤2 (R − R0 ) R0 R rˆ 1 2 ∂ rˆ 0 2 R0 2 2β−1 Z R Z ω 1 ∂q(˜ ˆ 1 r , θ) R r, + 22β ω 2 (R − R0 ) r˜2β r˜ dθˆ d˜ 1 r˜ ∂ θˆ 2 R0 0 R 0 2
1 ≤ (R − R0 )2 R0 2
Z
ω
Z
R
which directly implies
kqk0,β−1,Ω1 ≤ Ck∇qk0,β,Ω1 ,
(4.13)
where 1 ω2 1 C = 22β+1 (R − R0 )R2β−1 R0−2β−1 2(R − R0 ) + . 2 2 R0 Combining inequalities (4.9) and (4.13) completes the lemma. We now consider a similar Poincar´e inequality for the vector case. Again, consider Ω = Ωw , where ∂Ω is partitioned into Dirichlet and Neumann boundaries, ΓD and ΓN respectively. The following lemma is valid for the pure Dirichlet and Neumann cases and for the mixed boundary condition cases when ΓD includes a part of the boundary adjacent to the origin and ω 6= 3π 2 . Lemma 4.5 Take Ω = Ωw and let u ∈ Hβ1 (Ω)2 satisfy τ · u = 0 on ΓD and n · u = 0 on ΓN . Assume that for the mixed boundary condition case that ω 6= 3π/2. Then there is a constant, C, dependent only on Ω, β and the length of the segments of Γ D and ΓN adjacent to the origin, such that kuk0,β−1 ≤ Ck∇uk0,β
(4.14)
for β > 0. Proof. First, consider the case when τ · u = 0 on ∂Ω. Denote the part of ∂Ω aligned with θ = 0 as Γ1 and the part of ∂Ω aligned with θ = ω as Γ2 . Thus, u1 = 0 on Γ1 and τx u1 + τy u2 = 0 on Γ2 . Since u1 and τx u1 + τy u2 satisfy the conditions in lemma 4.4, we may use ku1 k0,β−1 ≤ Ck∇u1 k0,β and kτx u1 + τy u2 k0,β−1 ≤ Ck∇(τx u1 + τy u2 )k0,β .
FOSLS FOR PROBLEMS WITH SINGULARITIES
13
Further, take τy 6= 0, since τy = 0 corresponds to either ω = π, for which the result holds trivially since the boundary is smooth, or ω = 2π, which we do not consider. Now, kuk20,β−1 = ku1 k20,β−1 + ku2 k20,β−1 1 = ku1 k20,β−1 + 2 kτx u1 − τx u1 + τy u2 k20,β−1 τy ≤ (1 + 2
2 τx2 )ku1 k20,β−1 + 2 kτx u1 + τy u2 k20,β−1 τy2 τy
≤ C(k∇u1 k20,β + k∇(τx u1 + τy u2 )k20,β ) ≤ Ck∇uk20,β .
The case when n·u = 0 on ∂Ω is analogous since u2 = 0 on Γ1 and nx u1 +ny u2 = 0 on Γ2 . Also, when ω 6= π2 , 3π 2 , the case for mixed boundary conditions follows similarly using the result of lemma 4.4. The case for mixed boundary conditions when ω = π/2, follows from appealing to symmetry in the the pure Dirichlet problem for ω = π. Remark 4.6 Lemma 4.5 can be directly extended to more generally shaped domains. The proof of the scalar Poincar´e bounds in lemmas 4.1, 4.2, 4.3 and 4.4 are simplified when the domain has the shape of Ωw with only one irregular boundary point. Since we are primarily interested in a local result, proving lemma 4.5 in the simple domain is sufficient for our purposes. Consider the following scalar Poisson problem ∆p = f, in p = 0, on n · ∇p = 0, on
in Ωw : Ω, ΓD ,
(4.15)
ΓN .
We refer to system (4.15) as the pure Dirichlet problem when ∂Ω = ΓD ; the pure Neumann problem when ∂Ω = ΓN ; and the mixed boundary condition problem when ΓD includes the part of ∂Ω coinciding with one of either θ = 0 or θ = ω, and ΓN = ∂Ω\ΓD with ΓN 6= ∅. The following regularity results can be found in [23] and [21]. Lemma 4.7 Assume |1 − β| < π/ω for the pure Dirichlet problem, 0 < |1 − β| < π/ω for the pure Neumann problem and |1 − β| < π/2ω for the mixed boundary condition problem. Then, for every f ∈ Hβ0 (Ω), there exists a unique solution to (4.15), p ∈ Hβ2 (Ω) for the pure Dirichlet and mixed boundary condition cases and p ∈ Hβ2 (Ω)/R for the pure Neumann problem. Moreover, there exists a constant, C, independent of p, such that kpk2,β ≤ Ckf k0,β .
(4.16)
Proof. See Chapter 1 of [23] for the Dirichlet and Neumann problems and Chapter 2 of [21] for the mixed boundary problem. Define the subspace of functions in Hβ1 (Ω) satisfying the appropriate boundary conditions by Vβ = {v ∈ Hβ1 (Ω) : τ · v = 0 on ΓD , n · v = 0 on ΓN }.
14
LEE, MANTEUFFEL, AND WESTPHAL
We now prove a regularity result for functions in Vβ . Recall that the power of the singularity is defined as α = π/ω for Dirichlet or Neumann boundary conditions and α = π/(2ω) for mixed boundary conditions. Lemma 4.8 Consider domain Ω = Ωw . Then there is a positive constant, C, independent of u, such that, for |1 − β| < α, the following bound holds for all u ∈ V β : kuk1,β ≤ CkLuk0,β .
(4.17)
Proof. From lemma 4.7 we know that any u ∈ Vβ has the decomposition u = ∇φ + ∇⊥ ψ,
(4.18)
where φ, ψ ∈ Hβ2 (Ω) satisfy (
∆φ = ∇ · u, φ = 0,
in Ω, on ∂Ω,
(4.19)
and (
∆ψ = ∇ × u, n · ∇ψ = 0,
in Ω, on ∂Ω.
(4.20)
Then, by applying lemma 4.7 to problems (4.19) and (4.20) we have kuk1,β = k∇φ + ∇⊥ ψk1,β ≤ k∇φk1,β + k∇⊥ ψk1,β
≤ kφk2,β + kψk2,β ≤ C(k∇ · uk0,β + k∇ × uk0,β ) ≤ CkLuk0,β ,
which completes the proof. 5. Error Bounds. Let T h = ∪N i=1 τi be a quasi-uniform triangulation of polygonal domain Ω. Let I h represent standard interpolation onto a piecewise polynomial finite element space of degree k. From finite element theory, we have the following interpolation bounds. Lemma 5.1 Let Ω be a polygonal domain. There exists a constant, C, independent of v, such that, for all v ∈ H m (Ω), 1/2 X kv − I h vk2s,τ ≤ Chm−s |v|m (5.1) τ ∈T h
for 0 ≤ s ≤ m and 1 < m. Here, I h denotes interpolation by a piecewise polynomial of degree k = m − 1. (Note: here the norm k · ks,τ is the standard H s (τ ) norm.) Proof. See [9] or [7]. We now consider a weighted interpolation bound for functions on domains with a polygonal corner at the origin. Define the modified interpolation operator, I 0h , by h Pn I u = i=0 u(ai )φi , if τ does not intersect the origin, P I0h u|τ = n if τ intersects the origin, i=1 u(ai )φi ,
where I h is a standard polynomial interpolation operator, φi are basis functions corresponding to the n + 1 nodal points, ai , and a0 is the origin, (0, 0). Thus, the modified interpolation has a value of zero at the origin and resembles I h away from the origin.
15
FOSLS FOR PROBLEMS WITH SINGULARITIES
Lemma 5.2 Let Ω be a polygonal domain. There exists a constant, C, independent of u, such that, for all u ∈ Hβm (Ω) satisfying equation (4.8),
X
τ ∈T
h
1/2
ku − I0h uk21,β,τ
≤ Chm−1 kukm,β ,
(5.2)
for 1 < m and β > 0, where I0h is the modified interpolation operator onto piecewise polynomials of degree k = m − 1 defined above. Proof. Define K0 = {τ | τ¯ ∩ (0, 0) 6= ∅} as p the set of elements adjacent √ to the origin. On T h \K0 , we have h ≤ rmin ≤ r = x2 + y 2 ≤ rmax ≤ rmin + 2h with rmin = inf{r|(x, y) ∈ τ } and rmax = sup{r|(x, y) ∈ τ } in τ , and ku − I0h uk21,β,τ = ku − I h uk21,β,τ Z = r2β |∇(u − I h u)|2 + r2(β−1) |u − I h u|2 dτ τ Z Z −2 2β 2β rmin |u − I h u|2 dτ ≤ rmax |∇(u − I h u)|2 dτ + rmax ≤
=
τ τ −2 2β 2(m−1) 2 2β 2m Crmax h |u|m,0,τ + Crmax rmin h |u|2m,0,τ −2 2β 2β h2(m−1) |u|2m,0,τ h2(m−1) (1 + rmin Crmax h2 )|u|2m,0,τ ≤ Crmax
−2β 2β rmin ≤ Ch2(m−1) rmax
τ
√
r2β |Dm u|2 dτ
rmin + 2h rmin
≤ Ch2(m−1) ≤ Ch2(m−1)
Z
Z
τ
!2β Z
τ
r2β |Dm u|2 dτ
r2β |Dm u|2 dτ.
We now consider the case for which τ ∈ K0 . Let δ ∈ C ∞ be a cut-off function defined by
δ(r) =
1, 0,
if if
r ≤ h/3, r > 2h/3,
with |δ (m) | ≤ ch−m , where δ (m) is the mth derivative of δ. By the triangle inequality, ku − I0h uk1,β,τ ≤ kδu − I0h (δu)k1,β,τ + k(1 − δ)u − I0h ((1 − δ)u)k1,β,τ .
(5.3)
By the definition of δ we have I0h ((1 − δ)u) = I h ((1 − δ)u) and I0h (δu) = 0. For the second term in (5.3), we apply lemmas 4.4 and 5.1, the properties in δ, and Fubini’s
16
LEE, MANTEUFFEL, AND WESTPHAL
theorem to obtain k(1 − δ)u − I0h ((1 − δ)u)k21,β,τ = k(1 − δ)u − I h ((1 − δ)u)k21,β,τ Z ≤ Ck∇((1 − δ)u − I h ((1 − δ)u))k20,β,τ ≤ Ch2β |∇((1 − δ)u − I h ((1 − δ)u))|2 dτ τ Z X Z m 2 2(β+m−1) 2(β+m−1) |D ((1 − δ)u)| dτ ≤ Ch ≤ Ch |Dm−j (1 − δ)D j u|2 dτ τ |j|≤m
τ
≤ Ch2(β+m−1)
≤ Ch2(β+m−1)
ZZ
h 3
X
|j|≤m−1
X
|j|≤m−1
X Z
= Ch2(m−1)
2h 3
|j|≤m
τ
ZZ
|h|j|−m Dj u|2 dτ +
ZZ
h 3
2h 3 h 3
r(θ)
|(1 − δ)D m u|2 dτ
ZZ −2β 2(β+|j|−m) j 2 h r |D u| dτ +
r(θ)
h 3
(r/h)2β |Dm u|2 dτ
r2(β+|j|−m) |Dj u|2 dτ = Ch2(m−1) kuk2m,β,τ .
Using the properties of δ and Fubini’s theorem result in a similar bound for the first term in (5.3): I0h (δu)k21,β,τ
Z
kδuk21,β,τ
= = kδu − r2β |∇(δu)|2 + r2(β−1) |δu|2 dτ τ Z ≤ C r2β (|∇δ · u|2 + |δ∇u|2 ) + r 2(β−1) |δu|2 dτ τ
≤C
ZZ
≤C
ZZ
≤C
Z
2h 3
h 3 2h 3 h 3
r
r2β h−2 |u|2 dτ + C
ZZ
r2(β−1) |u|2 dτ + C
ZZ
2(m−1)
(r
τ
Thus we have X
τ ∈Th
2(β−m+1)
2h 3
0 2h 3
0
r2β |∇u|2 + r2(β−1) |u|2 dτ r2β |∇u|2 + r2(β−1) |u|2 dτ
2
|∇u| + r2(β−m) |u|2 )dτ ≤ Ch2(m−1) kuk2m,β,τ .
ku − I0h uk21,β,τ ≤ Ch2(m−1)
X
τ ∈Th
kuk2m,β,τ ≤ Ch2(m−1) kuk2m,β ,
and the lemma follows. Lemma 5.3 Assume equation (4.14) holds in Ω. Then, for all uh ∈ V h , kuh k0,β ≤ Ch−η kuh k0,β+η for β > −1 and η > 0. Proof. Using lemma 4.5 and an inverse inequality, we may write kuh k0,β ≤ Ck∇uh k0,β+1 ≤ Ch−1 kuh k0,β+1 ,
(5.4)
FOSLS FOR PROBLEMS WITH SINGULARITIES
17
which establishes (5.4) for η = 1. Repeated application of this inequality thus validates (5.4) for any positive integer. Now consider 1/
1/
kuh k20,β = hr β uh , rβ uh i = hr β− 2 uh , rβ+ 2 uh i
≤ kuh k0,β−1/2 kuh k0,β+1/2 ≤ Ch−1 kuh k20,β+1/2 .
Taking the square root establishes (5.4) for η = 1/2. Repeating these steps leads to (5.4) for all η = ηn = kn /2`n for any nonnegative integers kn and `n . For any η > 0, choose a monotonically decreasing sequence, {ηn }, such that ηn > η and limn→∞ |ηn − η| = 0. Now, gn = (r β+ηn uh )2 is a monotonically increasing function that converges to g = (r β+η uh )2 pointwise everywhere. Thus, by the Lebesgue monotone convergence theorem, we have Z Z h 2 ku k0,β+η = gdx = lim gn dx = lim kuh k20,β+ηn n
n
and, therefore, kuh k0,β = lim kuh k0,β n
≤ lim Ch−ηn kuh k0,β+ηn n
= (lim Ch−ηn )(lim kuh k0,β+ηn ) n
= Ch
n
−η
h
ku k0,β+η ,
which completes the proof. Define an irregular boundary point of polygonal domain Ω to be a point on ∂Ω where interior angle ω satisfies ω > π when Dirichlet or Neumann boundary conditions are applied on both sides of the point or ω > π/2 when one Dirichlet boundary and a Neumann boundary meet at the corner. We now present error bounds for the numerical solution in weighted and unweighted norms. Theorem 5.4 Let Ω be a polygonal domain with one irregular boundary point of interior angle ω and let f ∈ L2 (Ω). Suppose u ∈ V satisfy Lu = f . If uh ∈ V h is chosen to minimize the weighted functional, Gw (uh ; f ) = kLuh − f k20,β = inf kLvh − f k20,β , vh ∈V h
for |1 − β| < α, then the approximation error, u − uh , satisfies the following bounds: ku − uh k1,β ≤ Chα+β−1 kukα+β,β , Gw (u − uh ; 0)
1/2
≤ Chα+β−1 kukα+β,β ,
(5.5) (5.6)
ku − uh k0,β ≤ Chs+β kukα+β,β ,
(5.7)
ku − uh k0 ≤ Chs kukα+β,β ,
(5.8)
where s < α, for α + β ≤ k + 1 with k the degree of the piecewise polynomial elements in V h .
18
LEE, MANTEUFFEL, AND WESTPHAL
Proof. By lemmas 4.8 and 5.2, we have ku − uh k1,β ≤ CkL(u − uh )k0,β ≤ CkL(u − I0h u)k0,β ≤ Cku − I0h uk1,β ≤ Chα+β−1 kukα+β,β ,
which establishes both (5.5) and (5.6) since we may write 1/
kL(u − uh )k0,β = Gw (u − uh ; 0) 2 . Note that lemmas 4.8 and 5.2 are satisfied for |1 − β| < α and α + β ≤ 2. For the weighted L2 -norm, we write X ku − uh k20,β = ku − uh k20,β,τ . τ ∈T h
The Cauchy inequality yields, for any ∈ (0, 1), 1− Z Z Z 1dτ r2β |u − uh |2 dτ ≤ ku − uh k20,β,τ = τ
τ
≤ Ch2−2
Z
τ
r2β |u − uh |2
dτ
2 ·2 β h = Ch2−2 krβ (u − uh )k2 2 r |u − u | dτ 2
τ
1
L (τ )
.
Since H 1 (Ω) is continuously imbedded into Lq (Ω) for all q ∈ [1, ∞), X X ku − uh k20,β = ku − uh k20,β,τ ≤ Ch2−2 krβ (u − uh )k2 2 τ ∈T h 2−2
≤ Ch
krβ (u − uh )k2 2
L (Ω)
τ ∈T h 2−2
≤ Ch
L (τ )
krβ (u − uh )k2H 1 (Ω)
≤ Ch2−2 ku − uh k21,β . Thus, by (5.5), we have ku − uh k0,β ≤ Ch1− ku − uh k1,β ≤ Chs+β kukα+β,β , where any s < α. We now consider the bound on ku − uh k0 . Let K0 = {τ | τ¯ ∩ (0, 0) 6= ∅} and K1 = T h \K0 . First, we consider the case β < 1. If τ ∈ K0 , then r ≤ Ch and r1−β ≤ Ch1−β . Thus, we have ku − uh k20,τ ≤ Ch2(1−β) ku − uh k20,β−1,τ ≤ Ch2(1−β) ku − uh k21,β,τ .
(5.9)
If τ ∈ K1 , we use the technique above to get ku − uh k20,τ ≤ Ch2(1−) ku − uh k2 2
L (τ )
.
Again, since H 1 is continuously imbedded into Lq for all q ∈ [1, ∞), we have X X ku − uh k2 2 ku − uh k20,τ ≤ Ch2−2 τ ∈K1
τ ∈K1
L (τ )
≤ Ch2−2 ku − uh k2 2 ≤ Ch2−2 ku − uh k2H 1 (K1 ) L (K1 ) Z 2−2 r−2β r2β (|u − uh |2 + |∇(u − uh )|2 ) dΩ = Ch ≤ Ch
K1 2(1−β−)
ku − uh k21,β,K1 .
(5.10)
19
FOSLS FOR PROBLEMS WITH SINGULARITIES
Hence by (5.9), (5.10), and (5.5) we have ku − uh k20 =
X
τ ∈K0
ku − uh k20,τ +
≤ Ch2s kuk2α+β,β ,
X
τ ∈K1
ku − uh k20,τ ≤ h2(1−β−) ku − uh k21,β
where s < α. The proof for β ≥ 1 follows analogously. Remark 5.5 For the optimal finite element convergence of O(h) with respect to the weighted functional and H 1 norms, we select β = 2 − α. But theorem 5.4 also requires that β < 1 + α. Thus, when α ∈ [1/2, 1), we may use a weighting with β = 2 − α and expect optimal rates, but when α ∈ (0, 1/2), our theory only guarantees at best O(2α) convergence using β = 1 + α. Numerical results, however, indicate that values of β larger than the theory allows can be used to recover optimal rates. We explore this in the next section. 6. Computational results. In this section, we present some numerical examples of the weighted-norm procedure to validate the error bounds in the previous section. As a test problem, we minimize the weighted functional on the following L-shaped domain: Ω = (−0.5, 0.5)2 \ [0, 0.5) × (−0.5, 0], which yields α = π/ω = 2/3. Function 2 f is chosen so that the solution of this test problem is u = ∇(χ(r)r 3 sin(2θ/3)), where χ(r) = 1 for r < 1/8, χ(r) = 0 for r > 3/8, and χ(r) is C 2 smooth. Again, note that f ∈ L2 (Ω) but u ∈ / H 1 (Ω). Define the following measures of the accuracy of the computed solution, u h : nonweighted functional norm nonweighted L2 norm of the error nonweighted H 1 seminorm of the error weighted functional norm weighted L2 norm of the error weighted H 1 seminorm of the error
G
1/2
1/
= (k∇ · uh − f k20 + k∇ × uh k20 ) 2 ,
0 = ku − uh k0 ,
1 = |u − uh |1 ,
1/
1/
Gw2 = Gw (uh ; f ) 2 , 0w = ku − uh k0,β ,
1w = |u − uh |1,β .
Since α = 2/3, we choose the optimal weight parameter, β = 2 − α = 4/3, for our computations. Table 6.1 summarizes discretization error and convergence rates for β = 4/3. 1 Asymptotic convergence rates in Ω are found to be approximately O(h) for G w/2 2 and 1w , O(h2 ) for 0w and O(h 3 ) for 0 . The approximation does not converge in 1 / H 1 (Ω). either the 1 or G /2 measures since u ∈ To distinguish between behavior near to and away from the singularity, we consider the error of the solution above on a partitioning of Ω. Define Ω0 = Ω ∩ ( 83 , 58 )2 and Ω1 = Ω\Ω0 ; see figure 6.1. Table 6.2 summarizes the asymptotic discretization accuracy obtained at the finest mesh size in subdomains Ω0 and Ω1 . Away from the singularity we observe optimal accuracy in all measures. As expected, near the singularity, the solution fails to converge in the nonweighted functional and H 1 norms. The nonweighted L2 error 2 achieves accuracy of approximately O(h 3 ) near the singularity.
20
LEE, MANTEUFFEL, AND WESTPHAL 1
h 8−1 16−1 32−1 64−1 128−1 256−1 512−1
Gw/2 5.52 4.34 2.34 1.19 5.98e-01 3.00e-01 1.50e-01
Ratio
Rate
1.27 1.85 1.97 1.99 1.99 2.00
0.35 0.89 0.98 0.99 0.99 1.00
h 8−1 16−1 32−1 64−1 128−1 256−1 512−1
0w 3.08e-01 1.35e-01 4.07e-02 1.11e-02 2.98e-03 7.84e-04 2.06e-04
Ratio
Rate
2.28 3.32 3.67 3.72 3.80 3.81
1.19 1.73 1.88 1.90 1.93 1.93
1w 3.81 1.47 6.66e-01 2.97e-01 1.41e-01 6.74e-02 3.31e-02
Ratio
Rate
2.59 2.21 2.24 2.11 2.09 2.04
1.37 1.14 1.16 1.08 1.06 1.03
0 3.72e-01 1.93e-01 8.93e-02 4.93e-02 3.00e-02 1.87e-02 1.18e-02
Ratio
Rate
1.93 2.16 1.81 1.64 1.60 1.58
0.95 1.11 0.86 0.71 0.68 0.66
Table 6.1 Convergence of discretization error for weighted-norm FOSLS.
Ω1
Ω0
Fig. 6.1. L-shaped domain Ω and subdomains Ω0 and Ω1 .
Gw/2
1
G /2
1
1w
1
0w
0
Ω1
O(h)
O(h)
O(h)
O(h)
O(h2 )
O(h2 )
Ω0
O(h)
O(1)
O(h)
O(1)
O(h2 )
O(h 3 )
Ω
O(h)
O(1)
O(h)
O(1)
O(h2 )
O(h 3 )
2 2
Table 6.2 Accuracy in Ω0 , Ω1 , and Ω with β = 2 − α.
Figure 6.2 shows the first component of the exact solution, u1 , and the standard FOSLS approximation uh1 . Figure 6.3 shows the error of the first component of the approximated solution for the standard FOSLS and the weighted-norm FOSLS methods. We see that the error in the approximation in standard FOSLS is highest near the
21
FOSLS FOR PROBLEMS WITH SINGULARITIES
2
2
1
1 1
0
0.4
0.6
−2 0
0.4
0.2
0.8
−1
0.6
−2 0
1
0
0.8
−1
0.4
0.2
0.6 0.8 1
0.4
0.2 0.2
0.6 0.8
0
1
0
(b) Component uh 1 computed by standard FOSLS.
(a) Exact component u1 .
Fig. 6.2. Exact solution component u1 and solution component uh 1 approximated by standard FOSLS on the h = 32−1 mesh.
L2 error in u1
1
L2 error in u1
1
0.9
0.9 2.5
0.8
2.5 0.8
0.7
2
0.6
0.7
2
0.6
0.5
1.5
0.4
0.5
1.5
0.4 1
0.3 0.2
1
0.3 0.2
0.5 0.1
PSfrag replacements
0
0.5 0.1
0
0.2
0.4
PSfrag replacements 0.6 0.8 1
(a) Standard FOSLS.
0
0
0.2
0.4
0.6
0.8
1
(b) Weighted-norm FOSLS, β = 2 − α.
Fig. 6.3. Reduction of the pollution effect by the weighted-norm procedure. Each plot is the −1 mesh. error of solution component uh 1 on the h = 32
singularity, but remains large even away from the corner point. In the weighted-norm FOSLS implementation, the error remains large near the singularity, as we expect, but is now concentrated only near the corner point. The pollution effect is removed by the weighting procedure. There are many boundary value problems not directly covered by the theory presented here that are of interest. For example, Poisson’s equation with mixed boundary conditions on the domain used above has a value of α = 1/3. To recover optimal convergence for this problem, the weighted-norm method requires a value of β larger than theorem 5.4 allows. In other elliptic equations (e.g., Stokes or the linear elasticity equations), the value of α is generally smaller than for Poisson’s equation for the same domain and boundary condition type. In each of these cases, a larger β
22
LEE, MANTEUFFEL, AND WESTPHAL
value is necessary for optimal convergence. This leads us to consider using larger β than the theory allows. Consider the same example problem as above on uniform mesh sizes of h = 1/8, 1/16, ..., 1/512, and values of β ranging from 1/3 to 23/6. Figure 6.4 plots the convergence rate at the finest level for: the weighted functional 1 norm, Gw/2 ; the weighted L2 norm, 0w ; and the L2 norm, 0 . While the functional Convergence Rates vs. Β 2 weighted L2 error
weighted functional norm 1 0.67 L2 error
0 0
0.33
1
1.66
Β
3
4
Fig. 6.4. Convergence rates versus β. The shaded region indicates values of β for which the assumptions of theorem 5.4 are satisfied.
norm retains optimal accuracy for large values of β, the solution fails to converge in the weighted and nonweighted L2 measures for β & 3. This indicates that, although the weighted-norm approach seems to be more robust than the theory allows, large values of β should still be used with caution. The method presented here is applicable to a wide range of problems and provides an efficient alternative to more specialized techniques for treating singularities in boundary value problems. Further numerical results for other problems including systems and in three dimensions can be seen in a companion paper. REFERENCES [1] T. Apel and B. Heinrich, Mesh refinement and windowing near edges for some elliptic problem, SIAM J. Numer. Anal., 31 (1994), pp. 695–708. [2] I. Babuska, R. Kellogg, and J. Pitkaranta, Direct and inverse error estimates for finite elements with mesh refinements, Numer. Math., 33 (1979), pp. 447–471. [3] M. Berndt, T. Manteuffel, S. McCormick, and G. Starke, Analysis of first-order system least squares (fosls) for elliptic problems with discontinuous coefficients: part i, SIAM J. Numer. Anal., 43 (2005), pp. 386–408. [4] , Analysis of first-order system least squares (fosls) for elliptic problems with discontinuous coefficients: part ii, SIAM J. Numer. Anal., 43 (2005), pp. 409–436. [5] P. Bochev and M. Gunzburger, Analysis of least-squares finite element methods for the stokes equations, Math. Comp., 63 (1994), pp. 479–506. [6] , On least-squares finite element methods for the poisson equation and their connection to the dirichlet and kelvin principles, SIAM J. Numer. Anal., 43 (2005), pp. 340–362. [7] D. Braess, Finite Elements: Theory, Fast Solvers and Applications in Solid Mechanics, Cambridge, 2001. [8] J. Bramble, R. Lazarov, and J. Pasciak, A least-squares approach based on a discrete minus one inner product for first order systems, Math. Comp., 66 (1997), pp. 935–955.
FOSLS FOR PROBLEMS WITH SINGULARITIES
23
[9] S. Brenner and L. Scott, The Mathematical Theory of Finite Element Methods, SpringerVerlag, 1994. [10] Z. Cai and S. Kim, A finite element method using singular functions for the poisson equation: Corner singularities, SIAM J. Numer. Anal., 39 (2001), pp. 286–299. [11] Z. Cai, T. Manteuffel, and S. McCormick, First-order system least squares for the stokes equations, with application to linear elasticity, SIAM J. Numer. Anal., 34 (1997), pp. 1727– 1741. , First-order system least squares for velocity-vorticity-pressure form of the stokes [12] equations, with application to linear elasticity, Elect. Trans. Numer. Anal., 3 (1997), pp. 150–159. [13] Z. Cai, T. Manteuffel, S. McCormick, and S. Parter, First-order system least squares (FOSLS) for planar linear elasticity: Pure traction problem, SIAM J. Numer. Anal., 35 (1998), pp. 320–335. [14] Z. Cai, T. Manteuffel, S. McCormick, and J. Ruge, First-order system LL ∗ (FOSLL*): Scalar elliptic partial differential equations, SIAM J. Numer. Anal., 39 (2002), pp. 1418– 1445. [15] C. Cox and G. Fix, On the accuracy of least squares methods in the presence of corner singularities, Comp. and Maths. with Appls., 10 (1984), pp. 463–475. [16] G. Fix, S. Gulati, and G. Wakoff, On the use of singular functions with finite element approximations, J. Comput. Phys., 13 (1973), pp. 209–228. [17] P. Grisvard, Elliptic Problems in Nonsmooth Domains, Pitman, 1985. , Singularities in Boundary Value Problems, Springer-Verlag, 1992. [18] [19] S. Kim, T. Manteuffel, and S. McCormick, First-order system least squares (fosls) for spatial linear elasticity: Pure traction, SIAM J. Numer. Anal., 38 (2001), pp. 1454–1482. [20] V. Kondratiev, Boundary problems for elliptic equations in domains with conical or angular points, Trans. Moscow Math. Soc., 16 (1967), pp. 227–313. [21] V. Kozlov, V. Maz’ya, and J. Rossmann, Spectral Problems Associated with Corner Singularities of Solutions to Elliptic Equations, vol. 85, American Mathematical Society, 2001. [22] T. Manteuffel, S. McCormick, J. Ruge, and J. Schmidt, First-order system ll* (fosll*) for general scalar elliptic problems in the plane, SIAM J. Numer. Anal., 43 (2005), pp. 2098– 2120. [23] V. Maz’ya, S. Nazarov, and B. Plamenevskij, Asymptotic Theory of Elliptic Boundary Vaule Problems in Singularly Perturbed Domains, Volume I, vol. 111, Birkh¨ auser, 2000. Operator Theory Advances and Applications. [24] G. Strang and G. Fix, An analysis of the finite element method, Prentice-Hall, Englewood Cliffs, NJ, 1973.
