A 𝐶 0 INTERIOR PENALTY METHOD FOR 𝑚 TH-LAPLACE EQUATION

. In this paper, we propose a 𝐶 0 interior penalty method for 𝑚 th-Laplace equation on bounded Lipschitz polyhedral domain in R 𝑑 , where 𝑚 and 𝑑 can be any positive integers. The standard 𝐻 1 -conforming piecewise 𝑟 -th order polynomial space is used to approximate the exact solution 𝑢 , where 𝑟 can be any integer greater than or equal to 𝑚 . Unlike the interior penalty method in Gudi and Neilan [ IMA J. Numer. Anal. 31 (2011) 1734–1753], we avoid computing 𝐷 𝑚 of numerical solution on each element and high order normal derivatives of numerical solution along mesh interfaces. Therefore our method can be easily implemented. After proving discrete 𝐻 𝑚 -norm bounded by the natural energy semi-norm associated with our method, we manage to obtain stability and optimal convergence with respect to discrete 𝐻 𝑚 -norm. The error estimate under the low regularity assumption of the exact solution is also obtained. Numerical experiments validate our theoretical estimate.


Introduction
We consider the th-Laplace equation (−1)  ∆   =  in Ω, (1.1a) on simplicial meshes for any  in [16].The finite element space in [16] contains piecewise -th order polynomials with  ≥ 2   + 1.Therefore, the polynomial order of finite element space in [16] is quite big.Though up to this moment they have above mentioned restrictions, conforming   finite element spaces are desirable in both theoretical analysis and practice.In order to simplify the construction of   finite element space, alternative   nonconforming finite element space is introduced in several works.In [22], a   nonconforming finite element space (named Morley-Wang-Xu elements) is introduced for  ≤ .Besides, Hu and Zhang also considered the   nonconforming finite element space in [15] on triangular meshes for  = 2.The finite element space in [22] is generalized for  =  + 1 by Wu and Xu [25].Recently in [24], it is further generalized for arbitrary  and  but with stabilization along mesh interface in order to balance the weak continuity and the penalty terms.In order to obtain stability and optimal convergence in some discrete   -norm, Wang and Xu [22], Wu and Xu [24,25] propose to compute numerical approximation to   , such that their implementation may become quite complicated as  is large.The finite element spaces in [4,14] can be used to solve numerically (1.1) with any source term  ∈  − (Ω).However, the implementation of these conforming and nonconforming finite element spaces can be quite challenging for large .Virtual element methods have been investigated for (1.1).In [1], a conforming   virtual element method is introduced for convex polygonal domain in R 2 .The finite element space in [1] contains piecewise -th order polynomials, where  ≥ 2 − 1.The virtual element method in [1] needs strong assumption on regularity of  ( ∈  −+1 (Ω)) to achieve optimal convergence (see [1], Thm.4.2).In [7], a nonconforming   virtual element method is developed for bounded Lipschitz polyhedral domain in R  , where  can be any positive integer.The design of finite element space in [7], which contains piecewise -th ( ≥ ) order polynomials, is based on a generalized Green's identity for   inner product.It is assumed that  ≤  in [7].In [17], the virtual element method in [7] is extended for  > .Besides above numerical methods based on primary formulations of (1.1), a mixed formulation based on Helmholtz decomposition for tensor valued function is introduced in [19] for two dimensional domain.
We propose a  0 interior penalty method (2.2) for (1.1) for arbitrary positive integers  and .The finite element space of (2.2) is the standard  1 -conforming piecewise -th order polynomials, where  ≥ .The design of (2.2) avoids computing   of numerical solution on each element and high order normal derivatives of numerical solution along mesh interfaces.In fact, equation (2.2) only gets involved with calculation of high order multiplicity of Laplace of numerical solution (∆   ℎ for 1 ≤  ≤ ) and the gradient of high order multiplicity of Laplace of numerical solution (∇∆   ℎ for 0 ≤  ≤  − 1) on both elements and mesh interfaces.Therefore our method (2.2) can be easily implemented, even when  is large and  = 3.After proving (Thm.3.4) that discrete   -norm (see Def. 3.1) is bounded by the natural energy semi-norm associated with (2.2), we manage to show our method (2.2) has stability and optimal convergence on bounded Lipschitz polyhedral domain in R  with respect to the discrete   -norm, for any positive integers  and .Roughly speaking, we have where  ≥ 2 − 1.We refer to Theorems 3.6 and 3.7 for detailed descriptions on stability and optimal convergence.The design and analysis of our method (2.2) can be easily generalized for nonlinear partial differential equations with (−1)  ∆   as their leading term.We would like to point out that our method (2.2) is not a generalization of the interior penalty method for sixth-order elliptic equation ( = 3) in (3.4), (3.5) of [13].Actually, the method in [13] needs to calculate numerical approximation to  3 .
If the exact solution of (1.1) is under the low regularity assumption, Gudi et al. have applied the analysis technique from the a posteriori error analysis to derive the error estimates for the interior penalty methods for the 2nd-order, 4th-order and 6th-order elliptic equations under the low regularity assumptions in [12,13].In this paper, we shall extend the analysis by Gudi et al. for the proposed  0 interior penalty methods for (1.1) when  ≥ 2. Assuming  ∈   0 (Ω) for (1.1), we have where  ℎ is the  0 conforming finite element space and osc  ( ) is the oscillation term defined in (4.6).
The numerical method considered in this paper works for any positive integers  and  from theoretical viewpoint.It can be applied to the practical high order equations.For instance, the modeling for plates in linear elasticity results in consideration of fourth-order partial differential equations [11].Modeling in material science usually applies the fourth-order equation such as the Cahn-Hilliard equation [5,10] and the sixth-order equation such as the thin-film equations [3] and the phase field crystal model [9,21,23].Recently, an eighth-order equation was considered for the nonlinear Schrödinger equation in [18].As mentioned in [22], although there are rare practical applications for general higher order equations, the elliptic equations of order  = /2 in any dimension have been used in differential geometry [6].One can also extends the numerical methods and analysis for the solution of nonlinear Hamilton-Jacobi-Bellman equation and other phase-field models.
In the next section, we present the  0 interior penalty method.In Section 3, we prove stability and optimal convergence with respect to discrete   -norm (see Def. 3.1).In Section 4, we show the error estimates under the low regularity assumption of the exact solution.In Section 5, we provide numerical experiments.

𝐶 0 interior penalty method
In this section, firstly we give notations to define the  0 interior penalty method for (1.1).Then in Section 2.1, we derive the  0 interior penalty method for any  ≥ 1. Finally in Section 2.2, we provide concrete examples of the method for  = 1, 2, 3, 4.
Let T ℎ be a quasi-uniform conforming simplicial mesh of Ω.Here we define ℎ = max ∈T ℎ ℎ  where ℎ  is the diameter of the element  ∈ T ℎ .We denote by F ℎ , F int ℎ and F  ℎ the collections of all ( − 1)-dimensional faces, interior faces and boundary faces of T ℎ , respectively.Obviously, We introduce some trace operators.For any interior face  ∈ F int ℎ , let  − ,  + ∈ T ℎ be two elements sharing  .We denote by  − and  + the outward unit normal vectors along  − and  + , respectively.For scalar function  : Ω → R and vector field  : Ω → R  , which may be discontinuous across F int ℎ , we define the following quantities.For  − := |  − ,  + := |  + ,  − := |  − and  + := |  + , we define We also define
The  0 interior penalty method is to find  ℎ ∈  ℎ , such that for any  ℎ , where (2.3) Here the parameter  ≥ 1 shall be large enough but independent of ℎ.

Examples of 𝐶 0 interior penalty method
- = 1.The  0 interior penalty method for −∆ =  is to find - = 2.The  0 interior penalty method for ∆ 2  =  is to find Actually, equation (2.5) is the  0 interior penalty method which replaces the discontinuous finite element spaces in [20] with  0 finite element space.
- = 3.The  0 interior penalty method for −∆ 3  =  is to find  ℎ ∈  ℎ satisfying It is easy to see that (2.6) is quite different from the interior penalty method in [13].
The  0 interior penalty method for ∆ 4  =  is to find  ℎ ∈  ℎ satisfying

Analysis
In this section, firstly we prove Theorem 3.4, which states the discrete   -norm (see Def. 3.1) bounded by the natural energy semi-norm associated with the  0 interior penalty method (2.2).Then we prove Theorem 3.6, which shows the energy estimate of (2.2).Finally, we prove Theorem 3.7, which gives optimal convergence of numerical approximation to  in the discrete   -norm.Throughout this paper,  with or without a subscript denotes a positive constant depending only on the property of Ω, the shape regularity of the meshes and the degree of polynomial spaces.The constant  can take on different values in different occurrences.Definition 3.1.For any integers  ≥ 2, we define the discrete   -norm ‖‖ ,ℎ by For any  ∈ F int ℎ , there are two elements  − ,  + ∈ T ℎ sharing the common face  .We denote by  − := |  − and  + := |  + .We define .
For any  ∈ F  ℎ , we define .

Discrete 𝐻 𝑚 -norm bounded by natural energy semi-norm
The main result of Section 3.1 is Theorem 3.4, which shows that the discrete   -norm (see Def. 3.1) is bounded by the natural energy semi-norm associated with the  0 interior penalty method (2.2).The proof of Theorem 3.4 is based on Lemmas 3.2 and 3.3.Lemma 3.2.For any integers  ≥  ≥ 2, there is a constant  > 0 such that Proof.We choose  ∈ F ℎ arbitrarily.There is an orthonormal coordinate system {  }  =1 such that the   -axis is parallel to normal vector along  .Therefore  1 -axis, • • • ,  −1 -axis are all parallel to  .
We claim that for any 1 ≤  ≤ , there is a positive integer  ′ such that We prove (3.2) by induction.When  = 1, it is easy to see By discrete inverse inequality and the fact that  1 -axis, • • • ,  −1 -axis are all parallel to  , we have that Therefore we have )︁ .
Thus (3.2) holds when  = 1.We assume that (3.2) holds for any  < .Then by discrete inverse inequality and the fact that )︁ .
We have applied (3.6) for ∆ l−  ℎ to obtain last inequality.We also notice that for any 0 .
Therefore we have that (3.8) implies where  = 2 l + 1.Now we assume that 1 ≤  <  is an even number,  = 2 l and Then by similar argument in last paragraph, we have that (3.10) implies where  = 2 l.According to (3.3), (3.7)-(3.11)and the fact that  ∈ F ℎ is chosen arbitrarily, we can conclude that the proof is complete.
According to (3.1c) in Theorem 3.1 of [8], there is a constant  > 0 such that where r ≥ 2 is a positive integer.
Lemma 3.3.We define 2 m + 1 =  if  is an odd number, while 2 m =  if  is an even number.Then there is a positive constant  such that for any  ℎ ∈  ℎ .
It is easy to see Applying (3.12) to each     ℎ , we have )︁ .
According to Lemmas 3.2, 3.3 and the discrete Poincaré inequality, we immediately have the following Theorem 3.4.Theorem 3.4.For any integers  ≥  ≥ 1, there is a constant  > 0 such that for any  ℎ ∈  ℎ .‖ ℎ ‖ ,ℎ is introduced in Definition 3.1.We point out that the right hand side of the above inequality is the natural energy semi-norm associated with the method (2.2).

Energy estimate of 𝐶 0 interior penalty method
We provide Theorem 3.6, which shows energy estimate of  0 interior penalty method (2.2) with respect to the discrete   -norm (see Def. 3.1).Before we prove Theorem 3.6, we introduce Lemma 3.5.Lemma 3.5.For any integers  ≥  ≥ 2 and any spatial dimension  ≥ 1, there is a positive number  0 ≥ 1 such that for any  ℎ ∈  ℎ , (3.15) Proof.We prove (3.15) for  = 2 m ( is an even number) in the following.It is similar to prove (3.15) for  which is an odd number.
By Definition 2.2, it is easy to see that (3.15) holds.Therefore the proof is complete.
Theorem 3.6.When  = 1, the method (2.2) is well-posed such that For any  ≥ 2, there is a positive number  0 ≥ 1 which is the same as Lemma 3.5, such that if  ≥  0 , then the method (2.2) is well-posed such that (3.17) Here  ℎ ∈  ℎ is the numerical solution of the method (2.2).
Proof.By (2.4), the method is the standard finite element method for Poisson equation when  = 1.Therefore, the method (2.2) is well-posed and (3.16) holds, when  = 1.Now we consider  ≥ 2. By the definition of the bilinear form  ℎ (•, •), Theorem 3.4 and Lemma 3.5, the coercivity and the continuity of  ℎ (•, •) are obtained which imply the well-posedness of the method (2.2).We assume  = 2 m to be an even number.By taking  ℎ =  ℎ in the method (2.2), we have We choose  0 the same as Lemma 3.5.Then Lemma 3.5 implies Then by Theorem 3.4 and the above inequality, we obtain (3.17) when  is an even number.
It is similar to show that (3.17) holds when  is an odd number.Thus we can conclude that the proof is complete.

Error analysis of 𝐶 0 interior penalty method
We provide Theorem 3.7, which gives error analysis of  0 interior penalty method (2.2) with respect to the discrete   -norm (see Def. 3.1).
Proof.When  = 1, the method (2.2) is the standard finite element method (2.4) for Poisson equation with homogeneous Dirichlet boundary condition.So it is easy to see that (3.18) holds.In the following, we assume  ≥ 2. By Theorem 3.6, the method (2.2) has the unique numerical solution  ℎ ∈  ℎ .Since  ∈   0 (Ω), it is easy to see that for any 0 ≤  ≤  − 1, every component of    is continuous across any face  ∈ F int ℎ and is equal to zero along Ω.Therefore by Definitions 2.1 and 2.2, we have We denote by Π ℎ  ∈  ℎ the standard  2 -orthogonal projection of  into  ℎ .We define We assume  = 2 m to be an even number.By (2.1), (3.20) and the method (2.2), we have It is easy to see that We have used trace inequality and (3.21) to obtain the last inequality above.By trace inequality and (3.21) again, we have By inverse inequality and discrete trace inequality, Combing (3.23) with above estimates, we have We have used the fact that  is independent of ℎ to obtain the above inequality.Then by Theorem 3.4, we have Now we obtain the error estimate (3.19) when  ≥ 2 is an even number.It is similar to show that (3.19) holds when  ≥ 2 is an odd number.Therefore we can conclude that the proof is complete.

Error analysis under the low regularity assumption
In the above section, we assume the exact solution  ∈   0 (Ω) ∩   (Ω) with  ≥ 2 − 1.Since this regularity assumption may be high for the realistic problems, we further deduce the error analysis in this section for the exact solution under the low regularity assumption, i.e.,  ∈   0 (Ω).In this case, the Galerkin orthogonality does not hold true for the  0 interior penalty method if  ≥ 2. We derive the error analysis by the technique developed by Gudi in [12] which utilizes the analysis idea from the a posteriori error analysis.
Let  :=   0 (Ω) and   ℎ be the   -conforming finite element space in  .One can refer to the construction of   -conforming finite element space in  in any dimension according to a recent work in [16].For any ,  ∈  , let (, ) = (∆ m, ∆ m) Ω if  = 2 m and (, ) = (∇∆ m, ∇∆ m) Ω if  = 2 m + 1.As the three abstract assumptions in [12], firstly we assume there exists an enriching operator  ℎ :  ℎ →   ℎ such that Actually, for the cases of  = 2 and 3, this enriching operator  ℎ has been constructed by averaging technique and the above estimate has been derived in [12,13].Secondly, by the definition of  ℎ (•, •) in (2.3) and Lemma 3.5, choosing  as in Theorem 3.6, we easily have that Thirdly, we have the following estimate: for any  ∈  ,  ∈   ℎ and  ℎ ∈  ℎ , it holds that Actually, for  = 2 m, due to the fact that  ∈   ℎ and  ∈  , we can derive By the trace inequality and inverse estimate, we have which, together with (4.4), yields the estimate (4.3).For the case of  = 2 m + 1, one can similarly derive (4.3) and we omit the details here.By the estimates (4.1)-( 4.3) and following Lemma 2.1 in [12], we have In order to get the upper bound for the second term on the right-hand side of (4.5), we first provide two lemmas.
Lemma 4.1.Let  ∈  ℎ .There exists a positive constant  independent of mesh size such that where Proof.We provide the proof for the case of  = 2 m, and the case of  = 2 m + 1 can be similarly deduced.
By the above estimate and the triangular inequality, we directly obtain )︁ , which yields the desired estimate.
Lemma 4.2.Let  ∈  ℎ .For  = 2 m, there exists a positive constant  independent of mesh size such that, for For  = 2 m + 1, there exists a positive constant  independent of mesh size such that Proof.For brevity, we only provide the proof for the case of  = 2 m.The estimates (4.10) and (4.11) for the case of  = 2 m + 1 can be similarly deduced.The proof is based on the induction approach.Now we prove (4.8) with  = 0.For any  ∈ F int ℎ , we denote   =  − ∪  + where  − ∩  + =  .Let   be the unit normal vector along  pointing from  − to  + .Let  1 ∈  −1 (  ) be defined by For the construction of  1 , we can firstly assume  = .Let  + +1 and  − +1 be the linear basis functions at the nodes opposite to the face  on  + and  − respectively.We choose , where  + and  − are constants.It is obviously that  1 satisfies (4.12) and (4.13).One can easily choose  + and  − such that For  > , one can similarly construct ) such that  1 satisfies (4.12)-(4.14).Let  2 ∈   0 (  ) be a piecewise polynomial bubble function such that  2 (  ) = 1, where   is the barycenter of  .Denote  =  1  2 on   and extend it by zero on Ω ∖   .It follows from the definitions of  1 ,  2 and integration by parts that By the scaling argument, for  ∈   , we have which directly yields that By the inverse estimate, we have Combining (4.15) and (4.16) yields The above estimate (4.17) can be similarly deduced for the case of  ∈ Ω.Now combining (4.17), Lemma 4.1 and summing over all  ∈ F ℎ , we get the estimate (4.8) with  = 0. Next we prove (4.9) with  = 0.For any Here  1 can be similarly constructed as  1 .Let  2 ∈   0 (  ) be a piecewise polynomial bubble function such that  2 (  ) = 1, where   is the barycenter of  .Denote  =  1  2 on   and extend it by zero on Ω ∖   .By integration by parts, we have .
By the scaling argument, for  ∈   , we have which yields By the inverse estimate and the trace inequality, we get The above estimate can be similarly deduced for the case of  ∈ Ω.Combining the above estimate with (4.18), (4.17), Lemma 4.1 and summing over all  ∈ F ℎ , we obtain the estimate (4.9) with  = 0. We assume (4.8) and (4.9) hold true for 0 <  =  ≤ m − 2, we would like to prove that (4.8) and (4.9) hold true with  =  + 1.Since the derivations are similar, we only show the proof for (4.8) with  =  + 1.
For any  ∈ F int ℎ , let  1 ∈  −4−5 (  ) be defined by Here  1 can be similarly constructed as  1 .Let  2 ∈   0 (  ) be a piecewise polynomial bubble function such that  2 (  ) = 1, where   is the barycenter of  .Denote  =  1  2 on   and extend it by zero on Ω ∖   .By integration by parts, we have By the definition of , integration by parts yields We also have By the scaling argument, for  ∈   , we have which yields By the trace inequality, inverse estimate and (4.22), we obtain Similarly, by the trace inequality, inverse estimate and (4.22), we have For the estimate of  0 , following integration by parts, the trace inequality, inverse estimate and (4.22) yields Combining (4.19)-(4.21)and the upper bounds for  0 ,  1 and  2 , we have We note that For any  ∈ Ω, we can derive the similar estimate as (4.23).Now combining (4.23), the Cauchy-Schwarz inequality, the estimates of (4.8) and (4.9) with  ≤ , the inverse estimate and (4.22) and summing over all  ∈ F ℎ , we can finally obtain (4.8) with  =  + 1.
Now we can start to derive the upper bound for the second term on the right-hand side of (4.5).We provide the estimate for the case of  = 2 m, and the estimate for the case of  = 2 m + 1 can be similarly obtained.For any  ∈  ℎ ∖ {0}, let  =  −  ℎ .By the definition of  ℎ (•, •) in (2.3) and integration by parts, we have By the definition of  ℎ (•, •), the trace inequality, inverse estimate and (4.1), we further have Combining the above estimate, Lemmas 4.1 and 4.2, we directly have

Numerical experiments and discussions
In this section, we provide several numerical experiments to verify the theoretical prediction of the  0 interior penalty finite element method proposed in the previous sections in two and three dimensions.We calculate the rate of convergence of ‖ −  ℎ ‖ ,ℎ in various discrete   norms and compare each computed rate with its theoretical estimate.It is pointed out that the estimated convergence rates have very little dependency on the particular value when  = (1), so we choose  = 1 in the following tests.All the numerical experiments are carried out in C, and the resulting linear algebraic systems are solved using GMRES solvers from the PETSc package [2].
We list the errors along with their estimated rates of convergence in Tables 1 and 2 when  =  and  = +1, respectively.It is remarked that long double in C99 standard is used to represent extended precision floating point value for the 4th-order Laplacian operator, which is accurate up to 10 −20 .The tables indicate the following rates of convergence: ‖ −  ℎ ‖ ,ℎ = (ℎ), when  = , ‖ −  ℎ ‖ ,ℎ = (ℎ 2 ), when  =  + 1. Example 5.2.In the second example, we test the proposed method in which the solutions have partial regularity on a convex domain [13] and a non-convex one [25], respectively.To this end, we solve the third-Laplace equation (−∆) 3  = .
where (, ) are polar coordinates.Here  = 0 and  2 ∈  3+1/2 (Ω 2 ) due to the singularity at the origin.In both cases, the observed errors of the proposed method converge asymptotically with the optimal order ℎ and ℎ 1/2 , respectively, in the discrete  3 norm, as shown in Table 3.
We list the errors and rates of convergence in Table 4, which indicates that the computed solution converges asymptotically linearly to the exact solution in the discrete  3 norm.The observed rate is in agreement with Theorem 3.7.

Conclusion
A  0 interior penalty method is considered for th-Laplace equation on bounded Lipschitz polyhedral domain in R  in this paper.In order to avoid computing   of numerical solution on each element, we reformulate the  0 interior penalty method for the odd and even  respectively, and only the gradient and Laplace operators are used in the new method.A rigorous and detailed analysis is given for the key estimate that the discrete   -norm of the solution can be bounded by the natural energy semi-norm associated with our method.Then the stability estimate and the optimal error estimates with respect to discrete   -norm are achieved.The error estimate under the low regularity assumption of the exact solution is also provided.We believe that the proposed  0 interior penalty method for th-Laplace equation can be applied for the nonlinear high order partial differential equations which will be our consideration in future.