A new discretization for mth-Laplace equationswith arbitrary polynomial degreesThis work was supported by the Berlin Mathematical School.

# A new discretization for mth-Laplace equations with arbitrary polynomial degrees††thanks: This work was supported by the Berlin Mathematical School.

M. Schedensack Institut für Numerische Simulation, Universität Bonn, Wegelerstraße 6, D-53115 Bonn, Germany
###### Abstract

This paper introduces new mixed formulations and discretizations for th-Laplace equations of the form for arbitrary based on novel Helmholtz-type decompositions for tensor-valued functions. The new discretizations allow for ansatz spaces of arbitrary polynomial degree and the lowest-order choice coincides with the non-conforming FEMs of Crouzeix and Raviart for and of Morley for . Since the derivatives are directly approximated, the lowest-order discretizations consist of piecewise affine and piecewise constant functions for any Moreover, a uniform implementation for arbitrary is possible. Besides the a priori and a posteriori analysis, this paper proves optimal convergence rates for adaptive algorithms for the new discretizations.

Keywords th-Laplace equation, polyharmonic equation, non-conforming FEM, mixed FEM, adaptive FEM, optimality

AMS subject classification 31A30, 35J30, 65N30, 65N12, 74K20

## 1 Introduction

This paper considers th-Laplace equations of the form

 (−1)mΔmu=f (1.1)

for arbitrary Standard conforming FEMs require ansatz spaces in . To circumvent those high regularity requirements and resulting complicated finite elements, non-standard methods are of high interest [Mor68, EGH02, Bre12, GN11]. The novel Helmholtz decomposition of this paper decomposes any (tensor-valued) function in an th derivative and a symmetric part of a Curl. Given a tensor-valued function which satisfies in the weak sense, the projection of to the space of th derivatives then coincides with the th derivative of the exact solution of (1.1) (see Theorem 5.1 below). This results in novel mixed formulations and discretizations for (1.1). This approach generalises the discretizations of [Sch15, Sch16] from to .

The direct approximation of instead of enables low order discretizations; only first derivatives appear in the symmetric part of the Curl and so the lowest order approach only requires piecewise affine functions for any . In contrast to that, even interior penalty methods require piecewise quadratic [Bre12] resp. piecewise cubic [GN11] functions for resp. . Mnemonic diagrams in Figure 1 illustrate lowest-order standard conforming FEMs from [Žen70] and the lowest-order novel FEMs proposed in this work for . Since the proposed new FEMs differ only in the number of components in the ansatz spaces, an implementation of one single program, which runs for arbitrary order, is possible. In particular, the system matrices are obtained by integration of standard FEM basis functions.

For and the lowest polynomial degree in the ansatz spaces, discrete Helmholtz decompositions of [AF89, CGH14] prove that the discrete solutions are piecewise gradients (resp. Hessians) of Crouzeix-Raviart [CR73] (resp. Morley [Mor68]) finite element functions and therefore the new discretizations can be regarded as a generalization of those non-conforming FEMs to higher polynomial degrees and higher-order problems. The generalization of [WX13] of the non-conforming Crouzeix-Raviart and Morley FEMs to is restricted to a space dimension .

In the context of the novel (mixed) formulations, the discretizations appear to be conforming. The new generalization to higher polynomial degrees proposed in this paper appears to be natural in the sense that the inherent properties of the lowest order discretization carry over to higher polynomial ansatz spaces, namely an inf-sup condition, the conformity of the method, and a crucial projection property (also known as integral mean property of the non-conforming interpolation operator).

Besides the a priori and a posteriori error analysis, this paper proves optimal convergence rates for an adaptive algorithm, which are also observed in the numerical experiments from Section 7.

The remaining parts of this paper are organised as follows. Section 2 introduces some notation while some preliminary results are proved in Section 3. The proposed discretization of (1.1) in Section 5 is based on a novel Helmholtz decomposition for higher derivatives which is stated and proved in Section 4. Section 6 introduces an adaptive algorithm and proves optimal convergence rates. Section 7 concludes the paper with numerical experiments on fourth- and sixth-order problems.

Throughout this paper, let be a bounded, polygonal, simply connected Lipschitz domain. Standard notation on Lebesgue and Sobolev spaces and their norms is employed with scalar product . Given a Hilbert space , let resp.  denote the space of functions with values in whose components are in resp. . The space of infinitely differentiable functions reads and the subspace of functions with compact support in is denoted with . The piecewise action of differential operators is denoted with a subscript . The formula represents an inequality for some mesh-size independent, positive generic constant ; abbreviates . By convention, all generic constants do neither depend on the mesh-size nor on the level of a triangulation but may depend on the fixed coarse triangulation and its interior angles.

## 2 Notation

This section introduces notation related to higher-order tensors and tensor-valued functions and triangulations.

Define the set of -tensors over by

 X(ℓ):={R for ℓ=0,∏ℓj=1R2=R2×⋯×R2≅R2ℓ for ℓ≥1

and let denote the symmetric group, i.e., the set of all permutations of . Define the set of symmetric tensors by

 S(ℓ):={A∈X(ℓ)∣∀(j1,…,jℓ)∈{1,2}ℓ∀σ∈Sℓ:Aj1,…,jℓ=Ajσ(1),…,jσ(ℓ)}.

The symmetric part of a tensor is defined by

 (symA)j1,…,jℓ:=(card(Sℓ))−1∑σ∈SℓAjσ(1),…,jσ(ℓ)

for all , where denotes the number of elements in a set . For , the set coincides with the set of symmetric matrices, while for , the tensors consist of the four different components , , , and . Given -tensors and a vector , define the scalar product and the dot product by

 A:B :=∑(j1,…,jℓ)∈{1,2}ℓAj1,…,jℓBj1,…,jℓ, (A⋅q)j1,…,jℓ−1 :=Aj1,…,jℓ−1,1q1+Aj1,…,jℓ−1,2q2

for all . The following definition summarizes some differential operators. Recall that, for a Hilbert space , the space (resp. ) denotes the space of (resp. ) functions with components in .

###### Definition 1 (differential operators).

Let and and define by and . Define the th derivative of , the derivative , the divergence , the Curl, , and the curl, by

 (Dℓv)j1,…,jℓ :=∂ℓv/(∂xj1…∂xjℓ), (Dσ)j1,…,jℓ+1 :=∂σj1,…,jℓ/∂xjℓ+1, (Curlσ)j1,…,jℓ+1 :=(−1)jℓ+1∂σj1,…,jℓ/∂xp(jℓ+1), (divσ)j1,…,jℓ−1 :=∂σj1,…,jℓ−1,1/∂x1+∂σj1,…,jℓ−1,2/∂x2, (curlσ)j1,…,jℓ−1 :=−∂σj1,…,jℓ−1,1/∂x2+∂σj1,…,jℓ−1,2/∂x1

for .

For , these definitions coincide with the row-wise application of , , , and . The scalar product of tensor-valued functions is defined by . Given such that there exists with

 (ψ,Dℓv)L2(Ω)=(−1)ℓ(g,v)L2(Ω)for all v∈Hℓ0(Ω),

define the th order divergence of . The space is defined by

 H(divℓ,Ω):={ψ∈L2(Ω;X(ℓ))∣divℓψ∈L2(Ω)}.

Define furthermore for

 H(divℓ,Ω;X(k)):={ψ∈L2(Ω;X(k))∣∣∣∀(j1,…,jk−ℓ)∈{1,2}k−ℓ:ψj1,…,jk−ℓ,∙∈H(divℓ,Ω)}.
###### Remark 2.1.

Note that the existence of the th weak divergence does not imply the existence of any -th divergence for , e.g., for .

A shape-regular triangulation of a bounded, polygonal, open Lipschitz domain is a set of closed triangles such that and any two distinct triangles are either disjoint or share exactly one common edge or one vertex. Let denote the edges of a triangle and the set of edges in . Any edge is associated with a fixed orientation of the unit normal on (and denotes the unit tangent on ). On the boundary, is the outer unit normal of , while for interior edges , the orientation is fixed through the choice of the triangles and with and is the outer normal of on . In this situation, denotes the jump across . For an edge on the boundary, the jump across reads . For and , let

 Pk(T;X) :={v:T→X∣∣ each component of v is a polynomial of total degree≤k}; Pk(T;X) :={v:Ω→X|∀T∈T:v|T∈Pk(T;X)}

denote the set of piecewise polynomials and . Given a subspace , let denote the projection onto and let abbreviate . Given a triangle , let denote the square root of the area of and let denote the piecewise constant mesh-size with for all . For a set of triangles , let abbreviate

 ∥∙∥M:=√∑T∈M∥∙∥2L2(T).

## 3 Results for tensor-valued functions

The main result of this section is Theorem 3.2, which proves that defines a norm on the space defined in (3.5) below and can, thus, be viewed as a generalized Korn inequality. The following theorem is used in the proof of Theorem 3.2. Recall the definition of the Curl and the symmetric part of a tensor from Section 2.

###### Theorem 3.1.

Any satisfies

 ∥Curlγ∥L2(Ω)≲∥symCurlγ∥L2(Ω)+∥γ∥L2(Ω).
###### Proof.

The proof is subdivided in three steps.

Step 1. Let and with for all and for all , i.e.,

 j(k)=(1,…,1k,2,…,2m−k).

The combination of the definitions of and reads

 (symCurlγ)j(k)=card(Sm)−1∑σ∈Sm(−1)jσ(m)∂∂xp(jσ(m))γjσ(1),…,jσ(m−1). (3.1)

Let be the multi-index with the same number of ones and the number of twos reduced by one and the multi-index with the same number of twos and the number of ones reduced by one, i.e.,

 ¯j(k)=(1,…,1k,2,…,2m−k−1)andj–(k)=(1,…,1k−1,2,…,2m−k).

The symmetry of implies that if and if . Since the number of permutations such that is and the number of permutations such that is and since and , this implies that (3.1) equals

 (symCurlγ)j(k) (3.2) =card(Sm−1)card(Sm)⎛⎝(m−k)∂γ¯j(k)∂x−k∂γj–(k)∂y⎞⎠=m−km∂γ¯j(k)∂x−km∂γj–(k)∂y.

Step 2. This step applies [Neč67, Chap. 3, Thm. 7.6] and [Neč67, Chap. 3, Thm. 7.8] to operators defined below. Step 3 then proves a relation between these operators and the operator .

Define for , , and a multi-index

 ak,s,κ:=⎧⎨⎩−(k−1)/m if s=k−1 and% κ=(0,1),(m−k+1)/m if s=k and κ=(1,0),0 else.

Furthermore, define for

 Nksξ:=∑κ=(1,0),(0,1)ak,s,κξκ=⎧⎨⎩−(k−1)ξ2/m if s=k−1,(m−k+1)ξ1/m if s=k,0else

with the multi-index notation . Then the matrix reads

 1m⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝mξ1000…0−ξ2(m−1)ξ100…00−2ξ2(m−2)ξ10…0⋮⋮⋮⋱⋱⋮⋮⋮0…0−(m−2)ξ22ξ100…0−(m−1)ξ2ξ10…0−mξ2⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠∈R(m+1)×m.

If , the columns of this matrix are linear independent. Define the operators , , by

 Nkv:=m∑s=1∑κ=(1,0),(0,1)ak,s,κDκvs.

Then, the combination of [Neč67, Chap. 3, Thm. 7.6] with [Neč67, Chap. 3, Thm. 7.8] proves

 m∑s=1∥vs∥2H1(Ω)≲m+1∑k=1∥Nkv∥2L2(Ω)+m∑s=1∥vs∥2L2(Ω). (3.3)

Step 3. This step proves a relation between and for a proper choice of .

Define by setting for each the function with and (with for and for ). The symmetry of proves

 ∥Curlγ∥2L2(Ω) ≲m∑s=1∥vs∥2H1(Ω)% andm∑s=1∥vs∥2L2(Ω) ≈∥γ∥2L2(Ω). (3.4)

With the notation from Step 1 it holds that and the definition of from Step 2 and (3.2) reveal

 Nk+1v=(m−k)/m(∂vk+1/∂x)−k/m(∂vk/∂y)=(symCurlγ)j(k).

 m+1∑k=1∥Nkv∥2L2(Ω)≤∥symCurlγ∥2L2(Ω).

This, (3.4), and an application of (3.3) implies the assertion. ∎

Define, for , the spaces

 H(Ω,m−1) :={v∈H1(Ω;S(m−1))∣∣∫Ωvdx=0}, (3.5) Z :={β∈H(Ω,m−1)|symCurlβ=0}, Y :={γ∈H(Ω,m−1)∣∣∀β∈Z:(Curlβ,Curlγ)L2(Ω)=0}.

A computation reveals for , that the spaces and read

 Z ={γ∈H(1)∣∃c1∈R,c2∈R2 with γ(x)=c1x+c2}, (3.6) Y ={γ∈H1(Ω;R2)∣∣∫Ωγdx=0 and ∫Ωdivγdx=0}

 Z=⎧⎪⎨⎪⎩γ∈H(2)∣∣ ∣ ∣∣∃c1,c2,c3∈R,c4∈R2×2 with γ(x,y)=(c1x2+2c2xc1xy+c2y+c3xc1xy+c2y+c3xc1y2+2c3y)+c4⎫⎪⎬⎪⎭. (3.7)

The following theorem generalizes [CGH14, Lemma 3.3] from to higher-order tensors and states that defines a norm on . Note that .

###### Theorem 3.2.

Any satisfies

 ∥Curlγ∥L2(Ω)≲∥symCurlγ∥L2(Ω).
###### Proof.

Assume for contradiction that the statement does not hold. Then there exists a sequence with

 n∥symCurlγn∥L2(Ω)≤∥Curlγn∥L2(Ω)=1.

Since , Poincaré’s inequality implies that all components of are bounded in . Since is reflexive and compactly embedded in , there exists a subsequence (not relabelled) with a limit , in . This and Theorem 3.1 imply

 ∥Curl(γn−γℓ)∥L2(Ω) ≲∥symCurl(γn−γℓ)∥L2(Ω)+∥γn−γℓ∥L2(Ω) ≤1n+1ℓ+∥γn−γℓ∥L2(Ω)→0as n,ℓ→∞.

The Poincaré inequality and the completeness of imply the existence of with in and thus . It holds that and, therefore, defines a bounded functional on . Hence,

 ∥symCurlγ∥L2(Ω)=limn→∞∥symCurlγn∥L2(Ω)=0. (3.8)

Let . Since , the Cauchy inequality reveals

 (Curlβ,Curlγ)L2(Ω) =(Curlβ,Curl(γ−γn))L2(Ω) ≤∥Curlβ∥L2(Ω)∥Curl(γ−γn)∥L2(Ω)→0as n→∞.

This and (3.8) lead to and therefore . This contradicts and, hence, implies the assertion. ∎

###### Remark 3.3 (dependency on the domain).

The proof by contradiction from Theorem 3.2 does not provide information about the dependency on the domain. A scaling argument reveals that it does not depend on the size of the domain, but it may depend on its shape.

## 4 Helmholtz decomposition for higher orders

This section proves a Helmholtz decomposition of tensors into th derivatives and the symmetric part of a Curl in Theorem 4.4. This is a generalization of the Helmholtz decomposition of [BNS07] for fourth-order problems (). The proof is based on Theorem 4.1 below, which characterizes th-divergence-free smooth functions as symmetric parts of Curls.

###### Theorem 4.1.

Let and with . Then there exists with

 τ=symCurlγ.
###### Proof.

The proof is based on mathematical induction.

The base case is a classical result [Rud76]. Assume as induction hypothesis that the statement holds for , i.e., for all with there exists with .

The inductive step is split in five steps. Suppose that with .

Step 1. Then and . Let and . Recall the definition of the divergence from Definition 1. The symmetry of implies

 (divτ)j1,…,jm−1=∂τj1,…,jm−1,1/∂x1+∂τj1,…,jm−1,2/∂x2 =∂τjσ(1),…,jσ(m−1),1/∂x1+∂τjσ(1),…,jσ(m−1),2/∂x2=(divτ)jσ(1),…,jσ(m−1).

Hence, . The induction hypothesis guarantees the existence of with .

Step 2. This step defines some with .

The definitions of and from Section 2 for tensors combine to

 (symCurlβ)j1,…,jm−1 (4.1) =(card(Sm−1))−1∑σ∈Sm−1(−1)jσ(m−1)∂∂xp(jσ(m−1))βjσ(1),…,jσ(m−2).

Define by

 ˆβj1,…,jm:=(−1)p(jm)(card(Sm−1))−1∑σ∈Sm−1jσ(m−1)=p(jm)βjσ(1),…,jσ(m−2). (4.2)

The definition of implies

 (divˆβ)j1,…,jm−1 =(card(Sm−1))−12∑k=1(−1)p(k)∂∂xk∑σ∈Sm−1jσ(m−1)=p(k)βjσ(1),…,jσ(m−2).

Since if and only if , this equals

 (card(Sm−1))−12∑k=1∑σ∈Sm−1p(jσ(m−1))=k(−1)jσ(m−1)∂∂xp(jσ(m−1))βjσ(1),…,jσ(m−2) =(card(Sm−1))−1∑σ∈Sm−1(−1)jσ(m−1)∂∂xp(jσ(m−1))βjσ(1),…,jσ(m−2)

and, hence, the combination of the foregoing two displayed formulae with (4.1) leads to . The combination with Step 1 proves .

Step 3. Since , the base case (applied “row-wise” to ) guarantees the existence of with .

Step 4. This step shows .

Let be fixed and let and be the number of ones and twos. Then

 M1(jm):=N1−(2−jm)andM2(jm):=N2−(jm−1) (4.3)

are the numbers of ones and twos in . Define the index set

 T:={(k1,…,km−2)∈{1,2}m−2∣∣ ∣∣m−2∑ℓ=1kℓ=(