Mixed finite element approximation of the vector Laplacian with Dirichlet boundary conditions
aa r X i v : . [ m a t h . NA ] S e p MIXED FINITE ELEMENT APPROXIMATION OF THE VECTORLAPLACIAN WITH DIRICHLET BOUNDARY CONDITIONS
DOUGLAS N. ARNOLD, RICHARD S. FALK, AND JAY GOPALAKRISHNAN
Abstract.
We consider the finite element solution of the vector Laplace equation on adomain in two dimensions. For various choices of boundary conditions, it is known thata mixed finite element method, in which the rotation of the solution is introduced as asecond unknown, is advantageous, and appropriate choices of mixed finite element spaceslead to a stable, optimally convergent discretization. However, the theory that leads tothese conclusions does not apply to the case of Dirichlet boundary conditions, in which bothcomponents of the solution vanish on the boundary. We show, by computational example,that indeed such mixed finite elements do not perform optimally in this case, and we analyzethe suboptimal convergence that does occur. As we indicate, these results have implicationsfor the solution of the biharmonic equation and of the Stokes equations using a mixedformulation involving the vorticity. Introduction
We consider the vector Laplace equation (Hodge Laplace equation for 1-forms) on a two-dimensional domain Ω. That is, given a vector field f on Ω, we seek a vector field u suchthat(1.1) curl rot u − grad div u = f in Ω . (Notations are detailed at the end of this introduction.) A weak formulation of a boundaryvalue problem for this equation seeks the solution u in a subspace H ⊂ H (rot) ∩ H (div)satisfying(1.2) (rot u , rot v ) + (div u , div v ) = ( f , v ) , v ∈ H. If H is taken to be ˚ H (rot) ∩ H (div), the variational formulation implies the equation (1.1)together with the electric boundary conditions(1.3) u · s = 0 , div u = 0 on ∂ Ω . Magnetic boundary conditions, u · n = 0, rot u = 0, result if instead the subspace H in theweak formulation is taken to be H (rot) ∩ ˚ H (div). (The terms electric and magnetic derivefrom the close relation of the Hodge Laplacian and Maxwell’s equations.) If the domain Mathematics Subject Classification.
Primary: 65N30.
Key words and phrases. vector Laplacian, Hodge Laplacian, mixed finite elements.The work of the first author was supported in part by NSF grant DMS-0713568. The work of the secondauthor was supported in part by NSF grant DMS-0910540. The work of the third author was supported inpart by NSF grant DMS-1014817. This work was primarily carried out at the Institute for Mathematics andits Applications at the University of Minnesota where the authors were visiting.
Ω is simply-connected, both these boundary value problems are well-posed. (Otherwise, H contains a finite dimensional subspace consisting of vector fields which satisfy the boundaryconditions and have vanishing rotation and divergence with dimension equal to the numberof holes in the domain, and each problem can be rendered well-posed by replacing H withthe orthogonal complement of this space.)Even when the domain is simply connected, finite element methods based on (1.2) are prob-lematic. For example, on a non-convex polygon, approximations using continuous piecewiselinear functions converge to a function different from the solution of the boundary value.See [2, § mixed formulation with a stable choice of elements. The mixed formulationfor the electric boundary value problem seeks σ ∈ H , u ∈ H (div) such that( σ, τ ) − ( u , curl τ ) = 0 , τ ∈ H , (curl σ, v ) + (div u , div v ) = ( f , v ) , v ∈ H (div) . On a simply connected domain, this problem has a unique solution for any L vector field f ; u solves (1.1) and (1.3) and σ = rot u . To discretize, we choose finite element spacesΣ h ⊂ H , V h ⊂ H (div), indexed by a sequence of positive numbers h tending to 0, anddetermine σ h ∈ Σ h , u h ∈ V h by( σ h , τ ) − ( u h , curl τ ) = 0 , τ ∈ Σ h , (1.4) (curl σ h , v ) + (div u h , div v ) = ( f , v ) , v ∈ V h . (1.5)In order to obtain a stable numerical method, the finite element spaces Σ h and V h mustbe chosen appropriately. A stable method is obtained by choosing Σ h to be the Lagrangeelements of any degree r ≥ V h to be the Raviart–Thomas elements of the samedegree r (where the case r = 1 refers to the lowest order Raviart–Thomas elements). Inthe notation of [2], Σ h × V h = P r Λ × P − r Λ and the hypotheses required by [2] (the spacesbelong to a subcomplex of the Hilbert complex H −−→ H (div) div −→ L with bounded cochainprojections) are satisfied. From this it follows that the mixed finite element method is stableand convergent. Similar considerations apply to the magnetic boundary value problem,where the finite element spaces are ˚Σ h = Σ h ∩ ˚ H and ˚ V h = V h ∩ ˚ H (div) and the relevantHilbert complex is ˚ H −−→ ˚ H (div) div −→ L . Another possible choice is to take Σ h to beLagrange elements of degree r > V h to be Brezzi–Douglas–Marini elements of degree r − h × V h = P r Λ × P r − Λ ). This case is similar, and will not be discussed furtherhere.We turn now to the main consideration of the current paper, which is the equation (1.1)with Dirichlet boundary conditions u = 0 on ∂ Ω. This problem may of course be treatedin the weak formulation (1.2) with H = ˚ H (Ω; R ). In this case we may integrate by partsand rewrite the bilinear form in terms of the gradient (which, when applied to a vector, ismatrix-valued):(rot u , rot v ) + (div u , div v ) = (grad u , grad v ) , u , v ∈ ˚ H (Ω; R ) . Thus the weak formulation (1.2) is just(1.6) (grad u , grad v ) = ( f , v ) , v ∈ ˚ H (Ω; R ) , IXED FINITE ELEMENT APPROXIMATION OF THE VECTOR LAPLACIAN 3 for which the discretization using Lagrange or similar finite elements is completely standard.However, one might consider using a mixed method analogous to (1.4)–(1.5) for the Dirich-let boundary value problem in the hope of getting a better approximation of σ = rot u , orwhen Dirichlet boundary conditions are imposed on part of the boundary and electric and/ormagnetic boundary conditions are imposed on another part of the boundary. In fact, as wediscuss in Sections 4 and 5, a mixed approach to the vector Laplacian with Dirichlet bound-ary conditions is implicitly used in certain approaches to the solution of the Stokes equationswhich introduce the vorticity, and in certain mixed methods for the biharmonic equation. Inthe mixed formulation of the Dirichlet problem for the vector Laplacian, the vanishing of thenormal component is an essential boundary condition, while the vanishing of the tangentialcomponent arises as a natural boundary condition. No boundary conditions are imposed onthe variable σ . Thus, we define ˚ V h = V h ∩ ˚ H (div), and seek σ h ∈ Σ h , u h ∈ ˚ V h satisfying( σ h , τ ) − ( u h , curl τ ) = 0 , τ ∈ Σ h , (1.7) (curl σ h , v ) + (div u h , div v ) = ( f , v ) , v ∈ ˚ V h . (1.8)Note that curl Σ h * ˚ V h , so there is no Hilbert complex available in this case, and the theory of[2] does not apply. This suggests that there may be difficulties with stability and convergenceof the mixed method (1.7)–(1.8). In the next section, we exhibit computational examplesdemonstrating that this pessimism is well founded. The convergence of the mixed method forthe Dirichlet boundary value problem is severely suboptimal (while it is optimal for electricand magnetic boundary conditions). Thus, the difficulties arising from the loss of the Hilbertcomplex structure are real, not an artifact of the theory.However, the computations indicate that even for Dirichlet boundary conditions, the mixedmethod does converge, albeit in a suboptimal manner. While we do not recommend themixed formulation for the Dirichlet problem, in Section 3 we prove convergence at the sub-optimal rates that are observed and, in so doing, clarify the sources of the suboptimality.Theorem 3.1 summarizes the main results of our analysis, and the remainder of the sectiondevelops the tools needed to establish them.This analysis of the mixed finite element approximation of the vector Laplacian has impli-cations for the analysis of mixed methods for other important problems: for the biharmonicequation using the Ciarlet–Raviart mixed formulation, and for the Stokes equations usinga mixed formulation involving the vorticity, velocity, and pressure, or, equivalently, using astream function-vorticity formulation. As a simple consequence of our analysis of the vectorLaplacian, we are able to analyze mixed methods for these problems, elucidating the sub-optimal rates of convergence observed for them, and establishing convergence at the ratesthat do occur. Some of the estimates we obtain are already known, while others improveon existing estimates. The biharmonic problem is addressed in Section 4 and the Stokesequations in Section 5.We end this introduction with a summary of the main notations used in the paper. Forsufficiently smooth scalar-valued and vector-valued functions σ and u , respectively, we use DOUGLAS N. ARNOLD, RICHARD S. FALK, AND JAY GOPALAKRISHNAN the standard calculus operatorsgrad σ = ( ∂σ∂x , ∂σ∂y ) , curl σ = ( ∂σ∂y , − ∂σ∂x ) , div u = ∂u ∂x + ∂u ∂y , rot u = ∂u ∂x − ∂u ∂y . We use the standard Lebesgue and Sobolev spaces L p (Ω), H l (Ω), W lp (Ω), and also the spaces H (div , Ω) and H (rot , Ω) consisting of L vector fields u with div u in L or rot u ∈ L ,respectively. Since the domain Ω will usually be clear from context, we will abbreviate thesespaces as L p , H l , H (div), etc. For vector-valued functions in a Lebesgue or Sobolev space,we may use notations like H l (Ω; R ), although when there is little chance of confusion we willabbreviate this to simply H l . The closure of C ∞ (Ω) in H , H (div), and H (rot), are denoted˚ H , ˚ H (div), ˚ H (rot). Note that if u ∈ H (div), then the normal trace u · n ∈ H − / ( ∂ Ω)and ˚ H (div) = { u ∈ H (div) | u · n = 0 on ∂ Ω } . Similarly, ˚ H (rot) = { u ∈ H (rot) | u · s =0 on ∂ Ω } . We write ( · , · ) for the L (Ω) inner product (of either scalar- or vector-valuedfunctions) and k · k for the corresponding norm.We shall also need the dual space of ˚ H (div), the space ˚ H (div) ′ , normed by(1.9) k v k ˚ H (div) ′ := sup w ∈ ˚ H (div) ( v , w ) k w k H (div) . Clearly,(1.10) L (Ω; R ) ⊂ ˚ H (div) ′ ⊂ H − (Ω; R )with continuous inclusions. 2. Some numerical results
We begin by considering the solution of the Hodge Laplacian (1.1) with electric boundaryconditions (1.3) using the mixed method (1.4), (1.5). For the space Σ h , we use Lagrangefinite elements of degree r ≥ V h , Raviart–Thomas elements of degree r (consisting locally of certain polynomials of degree ≤ r , including all those of degree ≤ r − k u − u h k = O ( h r ) , k div( u − u h ) k = O ( h r ) , k σ − σ h k = O ( h r +1 ) , k grad( σ − σ h ) k = O ( h r ) . Table 2.1 shows the results of a computation with r = 2. Note that the computed ratesof convergence are precisely as expected. In the test problem displayed, the domain isΩ = (0 , × (0 ,
1) and the exact solution is u = (cos πx sin πy, πx cos πy ). The meshesused for computation were obtained by dividing the square into 2 n × n subsquares, n =1 , , , . . . r ≥ u = (sin πx sin πy, sin πx sin πy ). The finite IXED FINITE ELEMENT APPROXIMATION OF THE VECTOR LAPLACIAN 5
Table 2.1. L errors and convergence rates for degree 2 mixed finite elementapproximation of the vector Laplacian with electric boundary conditions. k u − u h k rate k div( u − u h ) k rate k σ − σ h k rate k curl( σ − σ h ) k rate2.14e-03 1.99 1.17e-02 1.99 2.16e-04 3.03 2.63e-02 1.985.37e-04 1.99 2.93e-03 2.00 2.70e-05 3.00 6.60e-03 1.991.34e-04 2.00 7.33e-04 2.00 3.37e-06 3.00 1.65e-03 2.003.36e-05 2.00 1.83e-04 2.00 4.16e-07 3.02 4.14e-04 2.00 element spaces are as for the computation of Table 2.1, except that the boundary conditionof vanishing normal trace is imposed in the Raviart–Thomas space V h . Note that the L rate of convergence for σ is not the optimal value of 3, but rather roughly 3 /
2. The L rateof convergence of curl σ (i.e., the H rate of convergence of σ ) is also suboptimal by roughly3 /
2: it converges only as h / . For u , the L convergence rate is the optimal 2, but theconvergence rate for div u is suboptimal by 1 / Table 2.2. L errors and convergence rates for degree 2 mixed finite elementapproximation of the vector Laplacian with Dirichlet boundary conditions. k u − u h k rate k div( u − u h ) k rate k σ − σ h k rate k curl( σ − σ h ) k rate1.22e-03 2.01 1.55e-02 1.58 1.90e-02 1.62 2.53e+00 0.633.05e-04 2.00 5.33e-03 1.54 6.36e-03 1.58 1.68e+00 0.607.63e-05 2.00 1.85e-03 1.52 2.18e-03 1.54 1.14e+00 0.561.91e-05 2.00 6.49e-04 1.51 7.58e-04 1.52 7.89e-01 0.53 We have carried out similar computations for r = 3 and 4 and for nonuniform meshes andthe results are all very similar: degradation of the rate of convergence by 3 / σ and curl σ ,and by 1 / u . However the case r = 1 is different. There we saw no degradation ofconvergence rates for uniform meshes, but for nonuniform meshes σ converged in L withrate suboptimal by 1 and curl σ did not converge at all.The moral of this story is that the mixed finite element method using the standard elementsis indeed strongly tied to the underlying Hilbert complex structure which is not present forthe vector Laplacian with Dirichlet boundary conditions, and the method is not appropriatefor this problem. Nonetheless the experiments suggest that the method does converge, albeitat a degraded rate. In the next section, we develop the theory needed to prove that thisis indeed so, and also to indicate where the lack of Hilbert complex structure leads to thesuboptimality of the method. 3. Error analysis
Theorem 3.1, which is the primary result of this section, establishes convergence of themixed method for the Dirichlet problem at the suboptimal rates observed in the previoussection. In it we suppose that Ω is a convex polygon endowed with a shape-regular andquasi-uniform family of triangulations of mesh size h . We continue to denote by Σ h ⊂ H DOUGLAS N. ARNOLD, RICHARD S. FALK, AND JAY GOPALAKRISHNAN and V h ⊂ H (div) the Lagrange and Raviart–Thomas finite element spaces of some fixeddegree r ≥
1, respectively, with ˚ V h = V h ∩ ˚ H (div). Theorem 3.1.
Let u denote the solution of the vector Laplace equation (1.1) with Dirichletboundary condition u = 0 , and let σ = rot u . There exist unique σ h ∈ Σ h , u h ∈ ˚ V h satisfyingthe mixed method (1.7) – (1.8) . If the polynomial degree r ≥ , then the following estimateshold for ≤ l ≤ r (whenever the norms on the right hand side are finite): k u − u h k ≤ Ch l k u k l , k div( u − u h ) k + k σ − σ h k + h k curl( σ − σ h ) k ≤ Ch l − / ( | ln h |k u k W l ∞ + k u k l +1 / ) . If r = 1 , the estimates are: k u − u h k ≤ Ch | ln h | ( | ln h |k u k W ∞ + k u k ) , k div( u − u h ) k + k σ − σ h k + h k curl( σ − σ h ) k ≤ Ch / ( | ln h |k u k W ∞ + h / k u k ) . Note that the error estimate for u is optimal order (modulo the logarithm when r = 1),while (again modulo the logarithm), the estimate for div u is suboptimal by 1 / σ and curl σ are suboptimal by 3 / C to denote a generic constantindependent of h , whose values may differ at different occurrences.The proof of this theorem is rather involved. Without the Hilbert complex structure, thenumerical method is not only less accurate, but also harder to analyze. The analysis willproceed in several steps. First, in Section 3.2, we will establish the well-posedness of thecontinuous problem, not in the space H × ˚ H (div), but rather using a larger space than H with weaker norm for σ . Next, in Section 3.3, we mimic the well-posedness proof on thediscrete level to obtain stability of the discrete problem, but with a mesh-dependent normon Σ h . This norm is even weaker than the norm used for the continuous problem, whichmay be seen as the cause of the loss of accuracy. To continue the analysis, we then introduceprojection operators into ˚ V h and Σ h and develop bounds and error estimates for them inSection 3.4. In Section 3.5 we combine these with the stability result to obtain basic errorestimates for the scheme, and we improve the error estimate for u h in Section 3.6 usingduality.3.1. Preliminaries.
First we recall two forms of the Poincar´e–Friedrichs inequality:(3.1) k τ k ≤ C P k curl τ k , τ ∈ ˚ H , k ψ k ≤ C P k grad ψ k , ψ ∈ ˆ H . Here ˆ H denotes the subspace of functions in H with zero mean. Similarly, we will use ˆ L to denote the zero mean subspace of L .Next we recall the Hodge decomposition. The space L (Ω; R ) admits a decomposition intothe orthogonal closed subspaces curl H and grad ˚ H , or, alternatively, into the subspacescurl ˚ H and grad H . The decomposition of a given v ∈ L according to either of these IXED FINITE ELEMENT APPROXIMATION OF THE VECTOR LAPLACIAN 7 may be computed by solving appropriate boundary value problems. For example, we maycompute the unique ρ ∈ ˚ H and φ ∈ ˆ H such that(3.2) v = curl ρ + grad φ, by a Dirichlet problem and a Neumann problem for the scalar Poisson equation, respectively:(3.3) (curl ρ, curl τ ) = ( v , curl τ ) , τ ∈ ˚ H , (grad φ, grad ψ ) = ( v , grad ψ ) , ψ ∈ ˆ H . Clearly, k grad φ k ≤ k v k . If v ∈ ˚ H (div), then φ satisfies the Neumann problem∆ φ = div v in Ω , ∂φ∂n = 0 on ∂ Ω , Z Ω φ dx = 0 , so, by elliptic regularity, k φ k ≤ C k div v k if the domain is convex, and k φ k ≤ C k div v k for any domain.We shall need analogous results on the discrete level. To this end, let S h denote the spaceof piecewise polynomials of degree at most r −
1, with no imposed interelement continuity.Then the divergence operator maps V h onto S h and also maps ˚ V h onto ˆ S h , the codimensionone subspace consisting of functions with mean value zero. The former pair of spaces is usedto solve the Dirichlet problem for the Poisson equation, and the later is used to solve theNeumann problem. Each pair forms part of a short exact sequence:(3.4) 0 → ˆΣ h curl −−→ V h div −→ S h → → ˚Σ h curl −−→ ˚ V h div −→ ˆ S h → , respectively.The usual Raviart–Thomas approximate solution to Poisson equation ∆ φ = g with Dirich-let boundary condition φ = 0 is then: find v h ∈ V h , φ h ∈ S h such that( v h , w ) + (div w , φ h ) = 0 , w ∈ V h , (div v h , ψ ) = ( g, ψ ) , ψ ∈ S h . Define the operator grad h : S h → V h by(grad h φ, w ) = − ( φ, div w ) , φ ∈ S h , w ∈ V h . From the stability of the mixed method, we obtain the discrete Poincar´e inequality k φ k ≤ ¯ C P k grad h φ k , φ ∈ S h , with ¯ C P independent of h . The solution ( v h , φ h ) ∈ V h × S h of themixed method may be characterized by(grad h φ h , grad h ψ ) = − ( g, ψ ) , ψ ∈ S h , and v h = grad h φ h .Corresponding to the first sequence in (3.4), we have the discrete Hodge decomposition(3.5) V h = curl Σ h + grad h S h , and corresponding to the second, the alternate discrete Hodge decomposition(3.6) ˚ V h = curl ˚Σ h + grad ◦ h S h , where grad ◦ h : S h → ˚ V h is defined by(grad ◦ h φ, w ) = − ( φ, div w ) , φ ∈ S h , w ∈ ˚ V h . DOUGLAS N. ARNOLD, RICHARD S. FALK, AND JAY GOPALAKRISHNAN
Both of the discrete Hodge decompositions can be characterized by finite element compu-tations. For example, in analogy to (3.3), for given v ∈ ˚ V h we may compute the unique ρ h ∈ ˚Σ h and φ h ∈ ˆ S h such that v = curl ρ h + grad ◦ h φ h from the following finite elementsystems (one primal, one mixed):(curl ρ h , curl τ ) = ( v , curl τ ) , τ ∈ ˚Σ h , (grad ◦ h φ h , grad ◦ h ψ ) = ( v , grad ◦ h ψ ) , ψ ∈ ˆ S h . Well-posedness of the continuous formulation.
As a first step towards analyz-ing the mixed method, we obtain well-posedness of a mixed formulation of the continuousboundary value problem for the vector Laplacian. To do so, we need to introduce a largerspace than H for the scalar variable, namelyΣ = { τ ∈ L : curl τ ∈ ˚ H (div) ′ } , with norm k τ k = k τ k + k curl τ k H (div) ′ (see (1.9)). The space Σ has appeared before instudies of the vorticity-velocity-pressure and stream function-vorticity formulations of theStokes problem [10], and an equivalent space (at least for domains with C , boundary) isused in [4]. The bilinear form for the mixed formulation is B ( ρ, w ; τ, v ) = ( ρ, τ ) − h curl τ, w i + h curl ρ, v i + (div w , div v ) , where h· , ·i denotes the pairing between ˚ H (div) ′ and ˚ H (div) (or more generally between aHilbert space and its dual.) Often, we will tacitly use the fact that if τ is in H , then h curl τ, w i = (curl τ, w ). Clearly, | B ( ρ, w ; τ, v ) | ≤ k ρ k + k w k H (div) ) / ( k τ k + k v k H (div) ) / , ρ, τ ∈ Σ , w , v ∈ ˚ H (div) , so B is bounded on (Σ × ˚ H (div)) × (Σ × ˚ H (div)). For τ ∈ Σ, we define τ ∈ ˚ H by(curl τ , curl ψ ) = h curl τ, curl ψ i , ψ ∈ ˚ H , Taking ψ = τ shows that(3.7) k curl τ k ≤ k curl τ k ˚ H (div) ′ ≤ k τ k Σ , τ ∈ Σ . It is also true that(3.8) k τ k Σ ≤ C ( k τ k + k curl τ k ) , τ ∈ Σ . To see this, define φ ∈ ˆ L by(3.9) ( φ, div v ) = h curl τ, v i − (curl τ , v ) , v ∈ ˚ H (div) . This is well-defined, since div ˚ H (div) = ˆ L , and, if div v vanishes, then v = curl ψ for some ψ ∈ ˚ H , so the right-hand side vanishes as well. Clearly, h curl τ, v i = (curl τ , v ) + ( φ, div v ) ≤ ( k curl τ k + k φ k ) k v k H (div) , v ∈ ˚ H (div) . Choosing v ∈ ˚ H in (3.9) with div v = φ and k v k ≤ C k φ k , we get k φ k ≤ C k curl( τ − τ ) k − .This implies k curl τ k ˚ H (div) ′ ≤ C ( k τ k + k curl τ k ), thus establishing (3.8). We concludefrom (3.7) and (3.8) that the norm τ
7→ k τ k + k curl τ k is an equivalent norm on Σ. IXED FINITE ELEMENT APPROXIMATION OF THE VECTOR LAPLACIAN 9
Assuming that f ∈ L (or even ˚ H (div) ′ ), we now give a mixed variational formulation ofthe continuous problem. We seek σ ∈ Σ, u ∈ ˚ H (div), such that( σ, τ ) − h curl τ, u i = 0 , τ ∈ Σ , h curl σ, v i + (div u , div v ) = ( f , v ) , v ∈ ˚ H (div) . We note that, if u ∈ ˚ H (Ω; R ) is the solution of the standard variational formulation (1.6)and σ = rot u , then σ , u solve this mixed variational formulation. Indeed, u ∈ ˚ H (Ω; R ) ⊂ ˚ H (div), σ ∈ L , and, for v ∈ ˚ H (Ω; R ) h curl σ, v i = ( σ, rot v ) = (rot u , rot v ) = ( f , v ) − (div u , div v ) . This implies that curl σ ∈ ˚ H (div) ′ , so σ ∈ Σ, and, extending to v ∈ ˚ H (div) by density, thatthe second equation above holds. Finally( σ, τ ) = (rot u , τ ) = h u , curl τ i for all τ ∈ L , so the first equation holds.In the next theorem, we establish well-posedness of the mixed variational problem byproving the inf-sup condition for B , following the approach of [1]. Note that the theoremestablishes well-posedness of the more general problem where the zero on the right hand sideof the first equation is replaced by the linear functional h g, τ i , where g ∈ Σ ′ , and we allow f ∈ ˚ H (div) ′ . Theorem 3.2.
There exist constants c > , C < ∞ such that, for any ( ρ, w ) ∈ Σ × ˚ H (div) ,there exists ( τ, v ) ∈ Σ × ˚ H (div) with B ( ρ, w ; τ, v ) ≥ c ( k ρ k + k w k H (div) ) , (3.10) k τ k Σ + k v k H (div) ≤ C ( k ρ k Σ + k w k H (div) ) . (3.11) Moreover, if w ∈ curl ˚ H , then we may choose v ∈ curl ˚ H .Proof. Define ρ ∈ ˚ H by (curl ρ , curl ψ ) = h curl ρ, curl ψ i , ψ ∈ ˚ H . Next, use the Hodgedecomposition to write w in the form w = curl µ + grad φ , with µ ∈ ˚ H and φ ∈ ˆ H , andrecall that(3.12) k grad φ k ≤ C k div w k . We then choose τ = ρ − δµ, v = w + curl ρ , where δ is a constant to be chosen. Hence, B ( ρ, w ; τ, v ) = k ρ k − δ ( ρ, µ ) − h curl ρ, w i + δ (curl µ, w )+ h curl ρ, w i + h curl ρ, curl ρ i + k div w k = k ρ k + δ k curl µ k − δ ( ρ, µ ) + k curl ρ k + k div w k . Recalling the constant C P in the Poincar´e inequality (3.1) and choosing δ sufficiently small,we obtain B ( ρ, w ; τ, v ) ≥ k ρ k + ( δ − δ C P / k curl µ k + k curl ρ k + k div w k ≥ c (cid:0) k ρ k + k w k H (div) (cid:1) , where we have used the facts that k w k = k curl µ k + k grad φ k , (3.12), and (3.8) in thelast step. This establishes (3.10).To establish (3.11), we observe that k v k H (div) ≤ k w k H (div) + k curl ρ k ≤ k w k H (div) + k ρ k Σ by (3.7), while k τ k Σ ≤ k ρ k Σ + δ k µ k Σ ≤ k ρ k Σ + δ k µ k ≤ C ( k ρ k Σ + k w k ) , since k µ k ≤ C k curl µ k ≤ C k w k .To establish the final claim, we observe that if w ∈ curl ˚ H , then obviously v = w +curl ρ ∈ curl ˚ H . (cid:3) Remark.
Had we posed the weak formulation using the space H × ˚ H (div) instead of Σ × ˚ H (div), we would not have obtained a well-posed problem.3.3. Stability of the discrete formulation.
In this section, we establish the stability ofthe mixed method (1.7)–(1.8), guided by the arguments used for the continuous problem inthe preceding subsection. Analogous to the norm on Σ, we begin by defining a norm on Σ h by k τ k h = k τ k + k curl τ k V ′ h , τ ∈ Σ h , where k v k ˚ V ′ h := sup w ∈ ˚ V h ( v , w ) k w k H (div) . The bilinear form is bounded on the finite element spaces in this norm: | B ( ρ, w ; τ, v ) | ≤ k ρ k h + k w k H (div) ) / ( k τ k h + k v k H (div) ) / , ρ, τ ∈ Σ h , w , v ∈ ˚ V h . For τ ∈ Σ h , we define τ ∈ ˚Σ h by(curl τ , curl ψ ) = (curl τ, curl ψ ) , ψ ∈ ˚Σ h . The discrete analogue of (3.7) again follows by choosing ψ = curl τ : k curl τ k ≤ k curl τ k ˚ V ′ h ≤ k τ k Σ h , τ ∈ Σ h . Next we establish discrete analogue of (3.8), that is,(3.13) k τ k Σ h ≤ C ( k τ k + k curl τ k ) , τ ∈ Σ h . To see this, define φ ∈ ˆ S h by( φ, div v ) = (curl τ, v ) − (curl τ , v ) , v ∈ ˚ V h . This is well-defined, since div ˚ V h = ˆ S h , and, if div v vanishes, then v = curl ψ for some ψ ∈ ˚Σ h , so the right-hand side vanishes as well. It follows that k curl τ k ˚ V ′ h ≤ k curl τ k + k φ k .To bound k φ k , as in the continuous case, we choose v ∈ ˚ H with div v = φ and k v k ≤ C k φ k . IXED FINITE ELEMENT APPROXIMATION OF THE VECTOR LAPLACIAN 11
In the discrete case, we also introduce Π Vh v , the canonical projection of v into the Raviart–Thomas space ˚ V h (see (3.22)), so div Π Vh v = P S h div v = φ and k v − Π Vh v k ≤ Ch k v k .Then k φ k = ( φ, div Π Vh v ) = (curl τ, Π Vh v ) − (curl τ , Π Vh v )= (curl τ, Π Vh v − v ) + (curl τ, v ) − (curl τ , Π Vh v ) ≤ Ch ( k τ k + k curl τ k − + k curl τ k ) k v k . Using the inverse inequality k τ k ≤ Ch − k τ k and the fact that k v k ≤ k φ k , gives the bound k φ k ≤ C ( k τ k + k curl τ k ), and implies (3.13).With this choice of norm, stability of the finite element approximation scheme is estab-lished by an argument precisely analogous to that used in the proof of Theorem 3.2, simplyusing the Σ h norm, the discrete gradient operator grad ◦ h , the discrete Hodge decomposition(3.6), the estimate (3.13), and the discrete Poincar´e inequality, instead of their continuouscounterparts. Theorem 3.3.
There exists constants c > , C < ∞ , independent of h , such that, for any ( ρ, w ) ∈ Σ h × ˚ V h , there exists ( τ, v ) ∈ Σ h × ˚ V h with B ( ρ, w ; τ, v ) ≥ c ( k ρ k h + k w k H (div) ) , k τ k Σ h + k v k H (div) ≤ C ( k ρ k Σ h + k w k H (div) ) . Moreover, if w ∈ curl ˚Σ h , then we may choose v ∈ curl ˚Σ h .Remark. Note that k τ k Σ h ≤ k τ k Σ for τ ∈ Σ h , but, in general equality does not hold. Hadwe used the Σ norm instead of the Σ h norm on the discrete level, we would not have beenable to establish stability.3.4. Projectors.
Our error analysis will be based on the approximation and orthogonalityproperties of certain projection operators into the finite element spaces: P S h : L → S h , P Σ h : H → Σ h , P ˚Σ h : ˚ H → ˚Σ h , P ˚ V h : ˚ H (div) → ˚ V h . For P S h , we simply take the L projection. By standard approximation theory, k s − P S h s k L p ≤ Ch l k s k W lp , ≤ l ≤ r, ≤ p ≤ ∞ . For P Σ h and P ˚Σ h , we use elliptic projections. Namely, for any τ ∈ H ,(curl P Σ h τ , curl ρ ) = (curl τ, curl ρ ) , ρ ∈ Σ h , ( P Σ h τ ,
1) = ( τ, , and, for any τ ∈ ˚ H (curl P ˚Σ h τ , curl ρ ) = (curl τ, curl ρ ) , ρ ∈ ˚Σ h . Then, by standard estimates,(3.14) k σ − P Σ h σ k + h k σ − P Σ h σ k ≤ Ch l k σ k l , ≤ l ≤ r + 1 . Moreover,(3.15) (curl[ σ − P Σ h σ ] , v ) ≤ Ch k curl( σ − P Σ h σ ) kk div v k , v ∈ V h , σ ∈ H . To prove this last estimate, we use the discrete Hodge decomposition (3.5) to write v =curl γ h +grad h ψ h , with γ h ∈ ˆΣ h and ψ h ∈ S h . As explained in § h ψ h , ψ h ) ∈ V h × S h is the mixed approximation of (grad ψ, ψ ) where ψ ∈ ˚ H solves ∆ ψ = div v in Ω.Since Ω is convex, k ψ k ≤ C k div v k . Therefore,(curl[ σ − P Σ h σ ] , v ) = (curl[ σ − P Σ h σ ] , curl γ h + grad h ψ h ) = (curl[ σ − P Σ h σ ] , grad h ψ h )= (curl[ σ − P Σ h σ ] , grad h ψ h − grad ψ ) ≤ Ch k curl( σ − P Σ h σ ) kk ψ k ≤ Ch k curl( σ − P Σ h σ ) kk div v k . For P ˚Σ h τ , τ ∈ ˚ H , we will use the W p estimate (due to Nitsche [15] for r ≥ r = 1; cf. also [5, Theorem 8.5.3]):(3.16) k τ − P ˚Σ h τ k W p ≤ Ch l − k τ k W lp , ≤ l ≤ r + 1 , ≤ p ≤ ∞ , which holds with constant C independent of p as well as h .We define the fourth projection operator, P ˚ V h : ˚ H (div) → ˚ V h , by the equations( P ˚ V h v , curl τ + grad ◦ h s ) = ( v , curl τ ) − (div v , s ) , τ ∈ ˚Σ h , s ∈ S h . In view of the discrete Hodge decomposition (3.6), P ˚ V h v ∈ ˚ V h is well defined for any v ∈ ˚ H (div). It may be characterized as well by the equations( v − P ˚ V h v , curl τ ) = 0 , τ ∈ ˚Σ h , (3.17) (div[ v − P ˚ V h v ] , s ) = 0 , s ∈ S h . (3.18)Similar projectors have been used elsewhere, e.g., [7, eq. (2.6)]. The properties of P ˚ V h aresummarized in the following theorem. Theorem 3.4.
For v ∈ ˚ H (div) and U ∈ ˚ H , (3.19) div P ˚ V h v = P S h div v , P ˚ V h curl U = curl P ˚Σ h U. Moreover, the following estimates hold k v − P ˚ V h v k L p ≤ Cph l k v k W lp , ≤ l ≤ r, ≤ p < ∞ , (3.20) k div( v − P ˚ V h v ) k ≤ Ch l k div v k l , ≤ l ≤ r, (3.21) whenever the norm on the right hand side is finite.Proof. The first commutativity property in (3.19) is immediate from (3.18), and the diver-gence estimate (3.21) follows immediately. For the second commutativity property, we notethat curl P ˚Σ h U ∈ ˚ V h and that, if we set v = curl U and replace P ˚ V h v by curl P ˚Σ h U , then thedefining equations (3.17), (3.18) are satisfied.To prove the L p estimate (3.20), we follow the proof of corresponding results for mixedfinite element approximation of second order elliptic problems given in [11]. First, we in-troduce the canonical interpolant Π Vh : H (Ω; R ) → V h into the Raviart–Thomas space,defined through the degrees of freedom(3.22) v Z e v · n w ds, w ∈ P r − ( e ) , v Z T v · w dx, w ∈ P r − ( T ) , IXED FINITE ELEMENT APPROXIMATION OF THE VECTOR LAPLACIAN 13 where e ranges over the edges of the mesh and T over the triangles. Then(3.23) k v − Π Vh v k L p ≤ Ch l k v k W lp , ≤ l ≤ r, ≤ p ≤ ∞ , and, since div Π Vh v = P S h div v ,(3.24) k div( v − Π Vh v ) k L p ≤ Ch l k div v k W lp , ≤ l ≤ r, ≤ p ≤ ∞ . Writing v − P ˚ V h v = ( v − Π Vh v ) + ( Π Vh v − P ˚ V h v ), it thus remains to bound the second term.From (3.19), div( P ˚ V h v − Π Vh v ) = 0, so P ˚ V h v − Π Vh v = curl ρ h for some ρ h ∈ ˚Σ h . Applyingthe decomposition (3.2), we have v − Π Vh v = curl ρ + grad ψ for some ρ ∈ ˚ H and ψ ∈ ˆ H .From (3.17), (curl ρ h , curl τ ) = (curl ρ, curl τ ) , τ ∈ ˚Σ h . Thus, ρ h = P ˚Σ h ρ and so satisfies the bound k curl ρ h k L p ≤ C k curl ρ k L p given above in (3.16).Since (curl ρ, curl τ ) = ( v − Π Vh v , curl τ ) = (rot( v − Π Vh v ) , τ ) , τ ∈ ˚ H ,ρ ∈ ˚ H satisfies − ∆ ρ = rot( v − Π Vh v ). Using the elliptic regularity result of [13, Corollary1], we have for 1 < p < ∞ that k ρ k W p ≤ C p k rot( v − Π Vh v ) k W − ,p ≤ C p k v − Π Vh v k L p . Following the proof of that result, the dependence of the constant C p on p arises from theuse of the Marcinkiewicz interpolation theorem for interpolating between a weak L and an L estimate. Using the explicit bound on the constant in this theorem found in [8, TheoremVIII.9.2], it follows directly that C p ≤ Cp , where C is a constant independent of p . Weremark that this regularity result requires the assumed convexity of Ω, and does not hold forall 1 < p < ∞ if Ω is only Lipschitz (c.f. [14]). Estimate (3.20) follows by combining theseresults and applying (3.23). (cid:3) Theorem 3.6 below gives one more property of P ˚ V h , inspired by an idea in [17]. To proveit we need a simple lemma. Lemma 3.5.
Let ρ be a piecewise polynomial function with respect to some triangulationwhich is nonzero only on triangles meeting ∂ Ω . Then for any ≤ q ≤ , k ρ k L q ≤ Ch /q − / k ρ k L , where the constant C depends only on the polynomial degree and the shape regularity of thetriangulation.Proof. By scaling and equivalence of norms on a finite dimensional space, we have k ρ k L q ( T ) ≤ Ch /q − k ρ k L ( T ) , ρ ∈ P r ( T ) , where the constant C depends only on the polynomial degree r and the shape constant forthe triangle T . Now, let T ∂h denote the set of triangles meeting ∂ Ω. Then k ρ k qL q (Ω) = X T ∈T ∂h k ρ k qL q ( T ) ≤ Ch − q X T ∈T ∂h k ρ k qL ( T ) . Applying H¨older’s inequality we have X T ∈T ∂h k ρ k qL ( T ) ≤ ( T ∂h ) (2 − q ) / (cid:0) X T ∈T ∂h k ρ k L ( T ) (cid:1) q/ , and T ∂h ≤ Ch − by the assumption of shape regularity. Combining these results gives thelemma. (cid:3) Theorem 3.6.
Let ≤ p ≤ ∞ . Then (3.25) ( v − P ˚ V h v , curl τ ) ≤ Ch − / − /p k v − P ˚ V h v k L p k τ k , τ ∈ Σ h , v ∈ ˚ H (div) ∩ L p . Proof.
Define ˚ τ ∈ ˚Σ h by taking the Lagrange degrees of freedom to be the same as those for τ , except setting equal to zero those associated to vertices or edges in ∂ Ω. Then k ˚ τ k ≤ C k τ k and τ − ˚ τ is nonzero only on triangles meeting ∂ Ω. By (3.17),( v − P ˚ V h v , curl τ ) = ( v − P ˚ V h v , curl[ τ − ˚ τ ]) . Let q = p/ ( p − ≤ q ≤
2. Applying H¨older’s inequality, the lemma, and an inverseinequality, we obtain( v − P ˚ V h v , curl( τ − ˚ τ )) ≤ k v − P ˚ V h v k L p k curl( τ − ˚ τ ) k L q ≤ C k v − P ˚ V h v k L p h / − /p k curl( τ − ˚ τ ) k L ≤ C k v − P ˚ V h v k L p h − / − /p k τ − ˚ τ k L , from which the result follows. (cid:3) Error estimates by an energy argument.
Using the projection operators definedin the last subsection and the stability result of the preceding section, we now obtain a basicerror estimate (which is not, however, of optimal order).
Theorem 3.7.
Let r ≥ denote the polynomial degree. There exists a constant C indepen-dent of the mesh size h and of p ∈ [2 , ∞ ) , for which k σ − σ h k + h k σ − σ h k + k u − u h k H (div) ≤ C h l − / − /p (cid:16) p k u k W lp + k u k l +1 / − /p (cid:17) , ≤ l ≤ r, if r ≥ , h / − /p (cid:16) p k u k W p + h / /p k u k (cid:17) , if r = 1 , whenever the norms on the right hand side are finite.Proof. We divide the errors into the projection and the remainder: σ − σ h = ( σ − P Σ h σ ) + ( P Σ h σ − σ h ) , u − u h = ( u − P ˚ V h u ) + ( P ˚ V h u − u h ) . Since, k σ − P Σ h σ k + h k σ − P Σ h σ k ≤ Ch t k σ k t ≤ Ch t k u k t +1 , ≤ t ≤ r + 1 , and, by Theorem 3.4, k u − P ˚ V h u k H (div) ≤ Ch t k u k t +1 , ≤ t ≤ r, IXED FINITE ELEMENT APPROXIMATION OF THE VECTOR LAPLACIAN 15 the projection error satisfies the necessary bounds (without the p k u k W lp term on the right-hand side).Therefore, setting ρ = σ h − P Σ h σ and w = u h − P ˚ V h u , it suffices to show that for 2 ≤ p < ∞ ,(3.26) k ρ k + k w k H (div) ≤ C (cid:0) k σ − P Σ h σ k + h k σ − P Σ h σ k + h − / − /p k u − P ˚ V h u k L p (cid:1) . Indeed, both cases of the theorem follow from (3.26), Theorem 3.4, and the inverse inequality Ch k ρ k ≤ k ρ k . By the stability result of Theorem 3.3, there exists ( τ, v ) ∈ Σ h × ˚ V h satisfying B ( ρ, w ; τ, v ) ≥ c (cid:0) k ρ k h + k w k H (div) (cid:1) , k τ k Σ h + k v k H (div) ≤ C (cid:0) k ρ k Σ h + k w k H (div) (cid:1) . By Galerkin orthogonality, B ( ρ, w ; τ, v ) = B ( σ − P Σ h σ, u − P ˚ V h u ; τ, v )= ( σ − P Σ h σ, τ ) − ( u − P ˚ V h u , curl τ ) + (curl( σ − P Σ h σ ) , v ) , where we used the definition of B and (3.19) in the last step. Applying the Cauchy-Schwarzinequality, Theorem 3.6, and (3.15), we then obtain B ( ρ, w ; τ, v ) ≤ C (cid:0) k σ − P Σ h σ k + h k curl( σ − P Σ h σ ) k + h − / − /p ) k u − P ˚ V h u k L p (cid:1) / × (cid:0) k τ k + k v k H (div) (cid:1) / . Together, these imply (3.26) and so complete the proof of the theorem. (cid:3)
Choosing p = | ln h | in the theorem gives a limiting estimate. Corollary 3.8.
The following estimates hold whenever the right hand side norm is finite: k σ − σ h k + h k σ − σ h k + k u − u h k H (div) ≤ C ( h l − / (cid:0) | ln h |k u k W l ∞ + k u k l +1 / (cid:1) , ≤ l ≤ r, if r ≥ , h / (cid:0) | ln h |k u k W ∞ + h / k u k (cid:1) , if r = 1 . For smooth solutions, choosing the maximum value of l = r in the corollary gives subop-timal approximation of σ by order h / , and suboptimal approximation of u and div u byorder h / (ignoring logarithms). In the next section, we show how to improve the L errorestimate for u to optimal order. The other estimates are essentially sharp, as demonstratedby the numerical experiments already presented.3.6. Improved estimates for u − u h . Using duality, we can prove the following estimatefor u − u h in L , which is of optimal order (modulo logarithms for r = 1). Theorem 3.9.
These estimates hold whenever the right hand side norm is finite: k u − u h k ≤ C ( h l k u k l , ≤ l ≤ r, if r ≥ ,h (cid:0) | ln h | / k u k W ∞ + k u k (cid:1) , if r = 1 . Proof.
Define φ ∈ Σ, w ∈ ˚ H (div) by B ( τ, v ; φ, w ) = ( v , u − u h ) , τ ∈ Σ , v ∈ ˚ H (div) . Thus w solves the Poisson equation − ∆ w = u − u h in Ω with homogeneous Dirichletboundary conditions, and φ = − rot w . Under our assumption that Ω is a convex polygon,we know that w ∈ H , φ ∈ H , and k φ k + k w k ≤ C k u − u h k .Choosing τ = σ − σ h and v = u − u h and then using Galerkin orthogonality, we obtain k u − u h k = B ( σ − σ h , u − u h ; φ, w ) = B ( σ − σ h , u − u h ; φ − P Σ h φ, w − P ˚ V h w ) . The right hand side is the sum of following four terms: T = ( σ − σ h , φ − P Σ h φ ) , T = − ( u − u h , curl[ φ − P Σ h φ ]) ,T = (curl[ σ − σ h ] , w − P ˚ V h w ) , T = (div[ u − u h ] , div[ w − P ˚ V h w ]) . We have replaced h· , ·i by the L -inner products because φ ∈ H and σ = rot u is in H whenever the right hand side norm in the theorem is finite. For T , we use the Cauchy–Schwarz inequality, the bound k φ − P Σ h φ k ≤ Ch k φ k ≤ Ch k u − u h k for the elliptic projection,and the estimate of Theorem 3.7 with p = 2 to obtain | T | ≤ C ( h l k u k l k u − u h k , ≤ l ≤ r, if r ≥ h ( k u k + h k u k ) k u − u h k , if r = 1.Similar considerations give the same bound for T .To bound T , we split it as ( P ˚ V h u − u , curl[ φ − P Σ h φ ]) and T ′ = ( u h − P ˚ V h u , curl[ φ − P Σ h φ ]).The first term is clearly bounded by Ch l k u k l k u − u h k , while, for the second, we use (3.15)to find that | T ′ | ≤ Ch k curl( φ − P Σ h φ ) kk div( u h − P ˚ V h u ) k . Bounding div( u h − P ˚ V h u ) via Theorem 3.7 and (3.21), we get | T ′ | ≤ Ch l k u k l k u − u h k , ≤ l ≤ r, | T ′ | ≤ Ch ( k u k + h k u k ) k u − u h k , r = 1 . Finally, we bound T . If r ≥
2, then we simply use the Cauchy–Schwarz inequality, thebound(3.27) k w − P ˚ V h w k ≤ Ch k w k ≤ Ch k u − u h k , and the p = 2 case of Theorem 3.7 to obtain | T | ≤ Ch l k u k l k u − u h k , ≤ l ≤ r. If r = 1, then (3.27) does not hold. Instead we split T as (curl[ σ − P Σ h σ ] , w − P ˚ V h w ) +(curl[ P Σ h σ − σ h ] , w − P ˚ V h w ). Since k w − P ˚ V h w k ≤ Ch k w k ≤ Ch k u − u h k , the first termis bounded by Ch k σ k k u − u h k ≤ Ch k u k k u − u h k . For the second, we apply Theorem 3.6and (3.20) to obtain | (curl[ P Σ h σ − σ h ] , w − P ˚ V h w ) | ≤ Ch − / − /p k w − P ˚ V h w k L p k P Σ h σ − σ h k≤ Ch / − /p p k w k W p k P Σ h σ − σ h k , ≤ p < ∞ . IXED FINITE ELEMENT APPROXIMATION OF THE VECTOR LAPLACIAN 17
By the Sobolev inequality, k w k W p ≤ K p k w k W q , where q = 2 p/ (2 + p ) <
2. Moreover, from[18] and a simple extension argument the constant K p ≤ Cp / . Since k w k W q ≤ C k w k with C depending only on the area of the domain, we obtain | (curl[ P Σ h σ − σ h ] , w − P ˚ V h w ) | ≤ Ch / − /p p / k P Σ h σ − σ h kk u − u h k , ≤ p < ∞ . By (3.14) and Theorem 3.7 with r = 1, k P Σ h σ − σ h k ≤ k σ − P Σ h σ k + k σ − σ h k ≤ C ( h / − /p p k u k W p + h k u k ) . Thus we obtain | T | ≤ C (cid:16) h − /p p / k u k W p + h / − /p p / k u k (cid:17) k u − u h k , ≤ p < ∞ , and, by choosing p = | ln h | and noting that h / | ln h | / is bounded, | T | ≤ Ch (cid:0) | ln h | / k u k W ∞ + k u k (cid:1) k u − u h k . The theorem follows easily from these estimates. (cid:3) The Ciarlet-Raviart Mixed Method for the Biharmonic
In this section, we show that the above analysis immediately gives estimates for theCiarlet–Raviart mixed method for the biharmonic, including some new estimates whichimprove on those available in the literature.Given g ∈ H − (Ω) = ( ˚ H (Ω)) ′ , the standard weak formulation of the Dirichlet problemfor the biharmonic seeks U ∈ ˚ H such that(∆ U, ∆ V ) = ( g, V ) , V ∈ ˚ H . Letting σ := − ∆ U ∈ L , we have ∆ σ = − g . Assuming that g ∈ H − (Ω), as we henceforthshall, for Ω a convex polygon, we have that U ∈ H (Ω), σ ∈ H (Ω) and k U k + k σ k ≤ C k g k − . Hence ( σ, U ) ∈ H × ˚ H satisfy( σ, τ ) − (curl U, curl τ ) = 0 , τ ∈ H , (curl σ, curl V ) = ( g, V ) , V ∈ ˚ H . We note that a mixed formulation in these variables, but with spaces that are less regular,can also be given for this problem (c.f. [4]), but we shall not pursue this approach here.The Ciarlet–Raviart mixed method [6] for the approximation of the Dirichlet problem forthe biharmonic equation using Lagrange elements of degree r , seeks σ h ∈ Σ h , U h ∈ ˚Σ h suchthat ( σ h , τ ) − (curl U h , curl τ ) = 0 , τ ∈ Σ h , (curl σ h , curl V ) = ( g, V ) , V ∈ ˚Σ h . This discretization has been analyzed in many papers under the assumption that Ω is aconvex polygon. It is proved in [12] and [3] that for r ≥ k U − U h k ≤ Ch r k U k r +1 , k σ − σ h k ≤ Ch r − k U k r +1 . The former estimate is optimal, while the estimate for k σ − σ h k is two orders suboptimal.The case r = 1 was analyzed in [17], where it was proven that k U − U h k ≤ Ch / | ln h | / k U k , k σ − σ h k ≤ Ch / | ln h |k U k . These estimates are suboptimal by 1 / / H regularity of U . (As noted in [17], the same technique could be applied for r ≥ / k σ − σ h k .) Below we improve the estimateon k U − U h k for r = 1 to an optimal order estimate (modulo logarithms), with decreasedassumptions on the regularity of the solution U .We now show how to obtain all of these results from the analysis of the previous section,with only minor modifications. Let u = curl U . Then B ( σ, u ; τ, curl V ) = ( g, V ) , ( τ, V ) ∈ H × ˚ H . Similarly, with u h = curl U h , B ( σ h , u h ; τ, curl V ) = ( g, V ) , ( τ, V ) ∈ Σ h × ˚Σ h . As above, set ρ = σ h − P Σ h σ ∈ Σ h , w = u h − P ˚ V h u ∈ ˚ V h . Note that w = curl U h − curl P ˚Σ h U ∈ curl ˚Σ h . Subtracting the above equations and writing v for curl V , we have B ( ρ, w ; τ, v ) = B ( σ − P Σ h σ, u − P ˚ V h u ; τ, v ) , ( τ, v ) ∈ Σ h × curl ˚Σ h . Since the stability result of Theorem 3.3 holds over the space Σ h × curl ˚Σ h , as stated in thelast sentence of the theorem, we can argue exactly as in proof of Theorem 3.7 and concludethat the estimates proved in that theorem for the Hodge Laplacian hold as well in thiscontext with one improvement. To estimate the term k u − P ˚ V h u k L p in (3.26), instead ofusing (3.20), we note that k u − P ˚ V h u k L p = k curl( U − P ˚Σ h U ) k L p and invoke (3.16). In thisway we avoid a factor of p . The improved estimates of Theorem 3.9 also translate to thisproblem, with essentially the same proof and a similar improvement. The dual problem is,of course, now taken to be: Find φ ∈ Σ , w ∈ curl ˚ H such that B ( τ, v ; φ, w ) = ( v , u − u h ) , τ ∈ Σ , v ∈ curl ˚ H . Thus w = curl W , where W solves the biharmonic problem ∆ W = rot( u − u h ) ∈ H − with Dirichlet boundary conditions, and φ = ∆ W . The relevant regularity result, valid ona convex domain, is k w k + k φ k ≤ C k W k ≤ C k rot( u − u h ) k − ≤ C k u − u h k . The remainder of the proof goes through as before, with the simplification that now theterms T and T ′ are zero, and the term k w − P ˚ V h w k L p can be bounded without introducinga factor of p as just described. The suppressed factors of p lead to fewer logarithms in thefinal result. Stating this result in terms of the original variable U instead of u = curl U , wehave the following theorem. IXED FINITE ELEMENT APPROXIMATION OF THE VECTOR LAPLACIAN 19
Theorem 4.1.
Let U solve the Dirichlet problem for the biharmonic equation, σ = − ∆ U ,and let U h ∈ ˚Σ h , σ h ∈ Σ h denote the discrete solution obtained by the Ciarlet–Raviart mixedmethod with Lagrange elements of degree r ≥ . If r ≥ and ≤ l ≤ r , then the followingestimates, requiring differing amounts of regularity, hold whenever the norms on the righthand side are finite: k σ − σ h k + h k σ − σ h k ≤ C (cid:26) h l − k U k l +1 h l − / (cid:0) k U k W l +1 ∞ + k U k l +3 / (cid:1) , k U − U h k ≤ Ch l k U k l +1 . If r = 1 , the estimates are: k σ − σ h k + h k σ − σ h k ≤ C (cid:26) ( k U k + h k U k ) ,h / (cid:0) k U k W ∞ + h / k U k (cid:1) , k U − U h k ≤ Ch ( | ln h | / k U k W ∞ + k U k ) . Stationary Stokes equations
Another application in which the vector Laplacian with Dirichlet boundary conditionsarises is the stationary Stokes equations, in which the vector field represents the velocity,subject to no-slip conditions on the boundary. A standard weak formulation (with viscosityequal to one) seeks u ∈ ˚ H (Ω , R ) and p ∈ ˆ L such that(grad u , grad v ) − ( p, div v ) = ( f , v ) , v ∈ ˚ H (Ω , R ) , (div u , q ) = 0 , q ∈ L . Mixed methods, such as we have discussed, have been used to approximate this problem,based on the vorticity-velocity-pressure formulation. For example, using the spaces definedin Section 3, the following weak formulation is discussed in [10]. Find σ ∈ Σ, u ∈ ˚ H (div), p ∈ ˆ L such that ( σ, τ ) − h curl τ, u i = 0 , τ ∈ Σ , h curl σ, v i − ( p, div v ) = ( f , v ) , v ∈ ˚ H (div) . (div u , q ) = 0 , q ∈ L . This formulation is obtained just as for the vector Laplacian, by writing(grad u , grad v ) = (rot u , rot v ) + (div u , div v )and introducing the variable σ = rot u . When f ∈ L (Ω; R ) and Ω is a convex polygon, u ∈ H (Ω; R ), p ∈ ˆ H (Ω), and σ = rot u ∈ H (Ω). Assuming this extra regularity, andsetting u = curl U , and v = curl V , ( σ, U ) ∈ H × ˚ H satisfy the stream function-vorticityequations: ( σ, τ ) − (curl U, curl τ ) = 0 , τ ∈ H (curl σ, curl V ) = ( f , curl V ) , V ∈ ˚ H . Taking g = rot f , this formulation coincides with the mixed formulation of the biharmonicproblem discussed in the previous section.We consider here the finite element approximation which seeks σ h ∈ Σ h , u h ∈ ˚ V h , p h ∈ ˆ S h such that ( σ h , τ ) − ( u h , curl τ ) = 0 , τ ∈ Σ h , (curl σ h , v ) − ( p h , div v ) = ( f , v ) , v ∈ ˚ V h . (div u , q ) = 0 , q ∈ ˆ S h . where the spaces Σ h , ˚ V h , and ˆ S h are defined as above. The existence and uniqueness of thesolution is easily established by standard methods. When f = 0, we get by choosing τ = σ h , v = u h , q = p h and adding the equations that σ h = 0 and div u h = 0. Hence u h = curl U h , U h ∈ ˚Σ h , and choosing τ = U h , we see that curl U h = 0. Since div ˚ V h = ˆ S h , we also get p h = 0.Error estimates for k u − u h k and k σ − σ h k are easily obtained by reducing the problem toits stream function-vorticity form and using the estimates obtained in the previous section.Letting u h = curl U h , and choosing v = curl V , V ∈ ˚Σ h , we see that ( σ h , U h ) is the uniquesolution of the Ciarlet-Raviart formulation of the biharmonic with g = rot f . Hence, theestimates for σ − σ h in Theorem 4.1 remain unchanged, except that we can replace k U k s by k u k s − . In particular, we have the following theorem. Theorem 5.1.
Let ( u , p ) solve the Dirichlet problem for the Stokes equation, σ = rot u , andlet u h ∈ ˚ V h , σ h ∈ Σ h , and p h ∈ ˆ S h denote the discrete solution obtained by the vorticity-velocity-pressure mixed method with r ≥ the polynomial degree. If r ≥ and ≤ l ≤ r ,then the following estimates, requiring differing amounts of regularity, hold whenever thenorms on the right hand side are finite: k σ − σ h k + h k σ − σ h k ≤ C (cid:26) h l − k u k l ,h l − / (cid:0) k u k W l ∞ + k u k l +1 / (cid:1) , k u − u h k ≤ Ch l k u k l . If r = 1 , the estimates are: k σ − σ h k + h k σ − σ h k ≤ C (cid:26) k u k + h k u k ,h / (cid:0) k u k W ∞ + h / k u k (cid:1) , k u − u h k ≤ Ch ( | ln h | / k u k W ∞ + k u k ) . The only item remaining is to derive error bounds for the approximation of the pressure.We obtain the following result, which gives error bounds that are suboptimal by O ( h / ). Theorem 5.2. If r ≥ and ≤ l ≤ r , then k p − p h k ≤ C ( h l − ( k u k l + k p k l − ) ,h l − / (cid:0) k u k W l ∞ + k u k l +1 / + k p k l − / (cid:1) . IXED FINITE ELEMENT APPROXIMATION OF THE VECTOR LAPLACIAN 21 If r = 1 , the estimates are k p − p h k ≤ C ( k u k + h k u k + k p k ,h / (cid:0) k u k W ∞ + h / k u k + k p k / (cid:1) . Proof.
From the variational formulation, we get the error equation( p h − P S h p, div v h ) = ( p − P S h p, div v h ) + (curl[ σ h − σ ] , v h ) , v h ∈ ˚ V h . We choose v ∈ ˚ H (Ω; R ) such that div v = p h − P S h p and k v k ≤ C k p h − P S h p k , and take v h = Π Vh v . We have that div v = div Π Vh v and k Π Vh v k H (div) ≤ C k v k ≤ C k p h − P S h p k , so k p h − P S h p k = ( p h − P S h p, div Π Vh v ) = ( p − P S h p, div Π Vh v ) + (curl[ σ h − σ ] , Π Vh v ) , = ( p − P S h p, p h − P S h p ) + (curl[ σ h − σ ] , Π Vh v − v ) + ( σ h − σ, rot v ) ≤ C ( k p − P S h p k + h k curl( σ h − σ ) k + k σ h − σ k ) k p h − P S h p k . It easily follows that k p − p h k ≤ C ( k p − P S h p k + k σ h − σ k + h k curl( σ h − σ ) k ) . The theorem follows directly by applying Theorem 5.1 and standard estimates for the errorin the L projection. (cid:3) A number of papers have been devoted to finite element approximation schemes of eitherthe vorticity-velocity-pressure or stream-function-vorticity formulation of the Stokes prob-lem. In particular, the lowest order ( r = 1) case of the method analyzed here was discussedin [9], (in which additional references can also be found). In the case of the magnetic bound-ary conditions, σ = u · n = 0, the authors established stability and first-order convergence,which is optimal, for all variables. But for the no-slip boundary conditions u = 0, withwhich we are concerned and which arise much more commonly in Stokes flow, they observein numerical experiments stability problems and reduced rates of convergence which are inagreement with the theory presented above.We close with a simple numerical example in the case r = 2 that demonstrates thatthe suboptimal convergence orders obtained above are sharp even for very smooth solutions.Our discretization of the vorticity-velocity-pressure mixed formulation of the Stokes problemthen approximates the velocity u by the second lowest order Raviart–Thomas elements, thevorticity σ by continuous piecewise quadratic functions, and the pressure p by discontinuouspiecewise linear functions. We take Ω to be the unit square and compute f corresponding tothe polynomial solution velocity field u = ( − x ( x − y (2 y − y − , y ( y − x (2 x − x − p = ( x − / + ( y − / . The computations, summarized inTable 5.1, indeed confirm the convergence rates established above, i.e., u h converges withoptimal order 2 to u in L , while the approximations to σ and curl σ are both suboptimalby 3 / p is suboptimal by 1 / Table 5.1. L errors and convergence rates for the mixed finite element ap-proximation of the Stokes problem for the vector Laplacian with boundaryconditions u · n = 0, u · s = 0 on the unit square. k u − u h k rate k p − p h k rate k σ − σ h k rate k curl( σ − σ h ) k rate3.26e-04 1.9 2.34e-03 1.3 2.70e-03 1.3 1.67e-01 0.28.35e-05 2.0 8.05e-04 1.5 9.70e-04 1.5 1.24e-01 0.42.10e-05 2.0 2.74e-04 1.6 3.47e-04 1.5 8.96e-02 0.55.27e-06 2.0 9.39e-05 1.6 1.24e-04 1.5 6.42e-02 0.5 References
1. Douglas N. Arnold, Richard S. Falk, and Ragnar Winther,
Finite element exterior calculus, homologicaltechniques, and applications , Acta Numer. (2006), 1–155. MR 2269741 (2007j:58002)2. , Finite element exterior calculus: from Hodge theory to numerical stability , Bull. Amer. Math.Soc. (N.S.) (2010), no. 2, 281–354. MR 2594630 (2011f:58005)3. I. Babuˇska, J. Osborn, and J. Pitk¨aranta, Analysis of mixed methods using mesh dependent norms , Math.Comp. (1980), no. 152, 1039–1062. MR 583486 (81m:65166)4. Christine Bernardi, Vivette Girault, and Yvon Maday, Mixed spectral element approximation of theNavier-Stokes equations in the stream-function and vorticity formulation , IMA J. Numer. Anal. (1992), no. 4, 565–608. MR 1186736 (93i:65113)5. Susanne C. Brenner and L. Ridgway Scott, The mathematical theory of finite element methods , thirded., Texts in Applied Mathematics, vol. 15, Springer, New York, 2008. MR 2373954 (2008m:65001)6. P. G. Ciarlet and P.-A. Raviart,
A mixed finite element method for the biharmonic equation , Mathemat-ical aspects of finite elements in partial differential equations (Proc. Sympos., Math. Res. Center, Univ.Wisconsin, Madison, Wis., 1974), Math. Res. Center, Univ. of Wisconsin-Madison, Academic Press, NewYork, 1974, pp. 125–145. Publication No. 33. MR 0657977 (58
Multigrid in a weighted space arisingfrom axisymmetric electromagnetics , Math. Comp. (2010), no. 272, 2033–2058. MR 26843548. Emmanuele DiBenedetto, Real analysis , Birkh¨auser Advanced Texts: Basler Lehrb¨ucher. [Birkh¨auser Ad-vanced Texts: Basel Textbooks], Birkh¨auser Boston Inc., Boston, MA, 2002. MR 1897317 (2003d:00001)9. F. Dubois, M. Sala¨un, and S. Salmon,
First vorticity-velocity-pressure numerical scheme for the Stokesproblem , Comput. Methods Appl. Mech. Engrg. (2003), no. 44-46, 4877–4907. MR 2013871(2004i:76141)10. Fran¸cois Dubois, Michel Sala¨un, and St´ephanie Salmon,
Vorticity-velocity-pressure and stream function-vorticity formulations for the Stokes problem , J. Math. Pures Appl. (9) (2003), no. 11, 1395–1451.MR 2020806 (2004i:35247)11. Ricardo G. Dur´an, Error analysis in L p , ≤ p ≤ ∞ , for mixed finite element methods for linear andquasi-linear elliptic problems , RAIRO Mod´el. Math. Anal. Num´er. (1988), no. 3, 371–387. MR 958875(89i:65114)12. R. S. Falk and J. E. Osborn, Error estimates for mixed methods , RAIRO Anal. Num´er. (1980), no. 3,249–277. MR 592753 (82j:65076)13. Stephen J. Fromm, Potential space estimates for Green potentials in convex domains , Proc. Amer. Math.Soc. (1993), no. 1, 225–233. MR 1156467 (93k:35076)14. David Jerison and Carlos E. Kenig,
The inhomogeneous Dirichlet problem in Lipschitz domains , J. Funct.Anal. (1995), no. 1, 161–219. MR 1331981 (96b:35042)15. J. A. Nitsche, L ∞ -convergence of finite element approximation , Journ´ees “´El´ements Finis” (Rennes,1975), Univ. Rennes, Rennes, 1975, p. 18. MR 568857 (81e:65058)16. Rolf Rannacher and Ridgway Scott, Some optimal error estimates for piecewise linear finite elementapproximations , Math. Comp. (1982), no. 158, 437–445. MR 645661 (83e:65180) IXED FINITE ELEMENT APPROXIMATION OF THE VECTOR LAPLACIAN 23
17. Reinhard Scholz,
A mixed method for 4th order problems using linear finite elements , RAIRO Anal.Num´er. (1978), no. 1, 85–90, iii. MR 0483557 (58 Best constant in Sobolev inequality , Ann. Mat. Pura Appl. (4) (1976), 353–372.MR 0463908 (57
Institute for Mathematics and its Applications and School of Mathematics, Universityof Minnesota, Minneapolis, MN 55455
E-mail address : [email protected] URL : Department of Mathematics, Rutgers University, Piscataway, NJ 08854
E-mail address : [email protected] URL : Portland State University, PO Box 751, Portland, OR 97207-0751
E-mail address : [email protected] URL ::