A C 0 linear finite element method for sixth order elliptic equations
AA C LINEAR FINITE ELEMENT METHOD FOR SIXTH ORDERELLIPTIC EQUATIONS
HAILONG GUO ∗ , ZHIMIN ZHANG † , AND
QINGSONG ZOU ‡ Abstract.
In this paper, we develop a straightforward C linear finite element method for sixth-order elliptic equations. The basic idea is to use gradient recovery techniques to generate higher-order numerical derivatives from a C linear finite element function. Both theoretical analysis andnumerical experiments show that the proposed method has the optimal convergence rate under theenergy norm. The method avoids complicated construction of conforming C finite element basis ornonconforming penalty terms and has a low computational cost. AMS subject classifications.
Primary 65N30; Secondary 45N08
Key words.
Sixth-order equation, Gradient recovery, Linear finite element.
1. Introduction.
Partial differential equations (PDEs) with order higher than2 have been widely used to describe different physical laws in material sciences [14,15,29, 30, 39], elastic mechanics [42], quantum mechanics [24], plasma physics [10, 11, 22],differential geometry [16, 41], and other areas of science and engineering. Comparingwith the second-order PDEs, higher order PDEs are much less studied, includingsome fundamental theoretical issues such as existence, uniqueness, and regularity ofsolutions.Numerical simulation becomes an important tool to study high order PDEs, andyet the design of efficient and reliable numerical methods is very challenging. As usual,the finite element method (FEM) plays a critical role in the numerical simulation.Both conforming and nonconforming methods have been applied to solve high orderPDEs in the literature. Usually, a conforming method requires higher regularityof the approximating functions (e.g. C functions for a fourth-order PDE and C for a sixth-order PDE), while a nonconforming method avoids the construction ofhigher regularity finite elements by adding some specially designed penalty termsto the scheme. The complicated construction of high regularity finite elements (forconforming methods) or penalty terms (for nonconforming methods) makes thesetwo FEMs hard to be implemented and significantly increases computational cost.Moreover, the analysis of the aforementioned FEMs is often very complicated.In this work, we present a systematic and simple numerical approach to treathigh-order PDEs and shed some light on theoretical analysis for this new method. Tobe more precise, we will develop a gradient recovery technique based C linear finite ∗ Department of Mathematics, University of California Santa Barbara, CA, 93106([email protected]). † Beijing Computational Science Research Center, Beijing 100094 and Department of Mathematics,Wayne State University, Detroit, MI 48202 ([email protected]). The research of this authorwas supported in part by the National Natural Science Foundation of China under grants 11471031,91430216, U1530401, and the U.S. National Science Foundation through grant DMS-1419040. ‡ Corresponding author. School of Data and Computers and Guangdong Province Key Laboratoryof Computational Science, Sun Yat-sen University, Guangzhou 510275 ([email protected]).This author was partially supported by the National Natural Science Foundation of China throughgrants 11571384 and 11428103, by Guangdong Provincial Natural Science Foundation of Chinathrough grant 2014A030313179, and by the Fundamental Research Funds for the Central Universitiesthrough the grant 16lgjc80. 1 a r X i v : . [ m a t h . NA ] A p r lement method for the following sixth-order equation −(cid:52) u = f in Ω (1.1) u = ∂ n u = ∂ nn u = 0 on ∂ Ω , (1.2)where Ω ⊂ R is an open bounded domain, f ∈ L (Ω), n is the outward unit normalof the boundary ∂ Ω. The sixth order derivative is defined as∆ u = ∆(∆(∆ u )) = (cid:88) i,j,k =1 ∂ u∂x i ∂x j ∂x k and the directional derivatives are ∂ n u = ∇ u · n , ∂ nn u = n T D u · n . The (weak)solution of (1.1)-(1.2) is a function u ∈ H (Ω) satisfying a ( u, v ) = ( f, v ) , ∀ v ∈ H (Ω) , (1.3)where the third order derivative tensor is given by D v = ∂ v∂x i ∂x j ∂x k and the bilinearform is a ( v, w ) = (cid:90) Ω D v : D w, ∀ v, w ∈ H (Ω) . Here the Frobenius product “:” for two tensors B = ( b ijk ) , B = ( b ijk ) is defined as B : B = (cid:88) i,j,k =1 b ijk b ijk . Note that the sixth-order elliptic boundary value problem (1.1) arises from manymathematical models including differential geometry( [16, 41]), the thin-film equa-tions [9], and the phase field crystal model [8, 18, 40]. For simplicity, we choose thehomogeneous Dirichlet boundary conditions. The basic principle can be applied toother boundary conditions as well.Let us illustrate the basic idea of the construction of our novel method. Note thatthe bilinear form a ( · , · ) involves the third derivative of the discrete solution, whichis impossible to obtain from a direct calculation of C linear element whose gradientis piecewise constant (w.r.t the underlying mesh) and discontinuous across each ele-ment. To overcome this difficulty, we use the gradient recovery operator G h to “lift”discontinuous piecewise constant Dv h to continuous piecewise linear function G h v h ,see [ ? , 2–5, 7, 20, 45, 47] for the details of different recovery operators. In other words,we use the special difference operator DG h to discretize the third order differential op-erator D . Our algorithm is then designed by applying this special difference operator to the standard Ritz-Galerkin method.From the above construction, our method has some obvious advantages. First,the fact that the recovery operator G h can be defined on a general unstructuredgrid implies that the method is valid for problems on arbitrary domains and meshes.Second, our method only has function value unknowns on nodal points instead of bothfunction value and derivative unknowns, its computational complexity is much lowerthan existing conforming and non-conforming methods in the literature.Naturally, one may question on the consistency, stability, and convergence of theproposed method, which require some more in-depth mathematical analysis. Let us egin with a discussion of consistency. As indicated in [45] (resp. [25]), for reasonablyregular meshes, G h u I (resp. G h u I ) is a second-order finite difference scheme of thegradient Du (resp. D u ), if u is sufficiently smooth. Here u I is the interpolation of u inlinear finite element space. As a consequence, DG h u I is a first-order approximation of D u , provided u is sufficiently smooth. However, for a discrete function v h in the finiteelement space which is not smooth across the element edges, the error (cid:107) Dv h − G h v h (cid:107) is not a small quantity of the high order, sometimes it may not converge to zero atall. Fortunately, an error estimate in [17, 26] set up a consistency in a weak sense,see (3.6) and (3.7) for the details. This weak consistency property of the gradientrecovery operator will play an important role in our error analysis.Next we discuss stability, which in our case can be reduced to verification of the(uniform) coercivity of the bilinear form ( DG h · , DG h · ) in the following sense (cid:107) v h (cid:107) (cid:46) (cid:107) DG h v h (cid:107) , (1.4)for all v h in the finite element space with suitable boundary conditions. Again, sincethe discrete Poincar´e inequality (3.5) has been established in [26], the stability (1.4)is a direct consequence. Note that (1.4) implies that no additional penalty term isneeded in order to guarantee the stability, and this fact makes our method very simple.The convergence properties of our method depends heavily on the aforementionedconsistency, stability, and the nice approximation properties of the recovery operator G h . As usual, the analysis of the error between the exact and approximate solutionscan be decomposed into the analysis of the approximation error and consistency error.Combining the weak consistency error estimates (3.6),(3.7) and approximation errorestimates (3.2)-(3.4) leads to the optimal convergence rate (= 1) under the energynorm ( H norm). This convergence rate is observed numerically. Furthermore, we alsonotice a second-order convergence rate under both H and L norms. However, we areonly able to prove a sub-optimal convergence rate under both the H and L normsat this moment. We would like to emphasize that our analysis here is straightforwardand simpler than the analysis of traditional conforming and nonconforming methodsapplied to sixth-order PDEs.The rest of the paper is organized as follows. We first present our algorithm inSection 2. Several numerical examples are provided in Section 3 to illustrate the effi-ciency and convergence rates of our algorithm. In Section 4, a rigorous mathematicalanalysis of our algorithm is given. Finally, some concluding remarks are presented inthe final section.
2. A recovery based C linear FEM. In this section, we discretize the vari-ational equation (1.3) in the standard C linear finite element space.Let T h be a triangulation of Ω with mesh-size h . We denote by N h and E h theset of vertices and edges of T h , respectively. Let V h be the standard P finite elementspace corresponding to T h . It is well-known that V h = Span { φ p : p ∈ N h } with φ p alinear nodal basis corresponding to each vertex p ∈ N h . Let G h : V h −→ V h × V h bea gradient recovery operator defined as below ( [34, 47]). For each vertex p ∈ N h , wedefine a recovered derivative ( G h v h )( p ) and let the whole recovered gradient functionbe G h v h = (cid:88) p ∈N h ( G h v h )( p ) φ p . For all v h ∈ V h , we have G h v h = ( G x h v h , G x h v h ) ∈ V h × V h . The corresponding ecovered Hessian matrix is defined as follows [25]: G h v h = (cid:18) G x h G x h v h G x h G x h v h G x h G x h v h G x h G x h v h (cid:19) . The derivative of G h v h is a tensor with its component( DG h v h ) ijk = ∂ x i G x j h G x k h v h , i, j, k = 1 , . For all v h , w h ∈ V h , we define a discrete bilinear form a h ( v h , w h ) = (cid:90) Ω D ( G h v h ) : D ( G h w h ) , The gradient recovery linear element scheme for solving (1.1) reads as : Find u h ∈ V h such that a h ( u h , v h ) = ( f, v h ) , v h ∈ V h , (2.1)where the homogenous finite element space V h = { v h ∈ V h | v h = G h v h · n = G h v h · t = n T G h v h n = 0 on ∂ Ω } . Note that here we use an additional condition G h v h · t = 0 since the exact solutionsatisfies ∂u∂ t = 0 on ∂ Ω, where t is the unit tangential vector on ∂ Ω. Remark 2.1.
For sixth order partial differential equation (1.1) - (1.2) , all threeboundary conditions are essential boundary conditions and we should incorporate suchtypes boundary conditions into the discretized linear system instead of the weak form.For partial differential equation (1.1) with nonhomogeneous boundary conditions u | ∂ Ω = g D , ∂ n u | ∂ Ω = g N , ∂ nn u | ∂ Ω = g R . (2.2) The variation form is to find u h ∈ V h with u h | ∂ Ω = g D , ∂ n u h | ∂ Ω = g N , ∂ nn u h | ∂ Ω = g R such that a h ( u h , v h ) = ( f, v h ) , v h ∈ V h , For numerical implement, there is no difference between homogeneous and nonho-mogeneous boundary conditions. In the article, we suppose homogeneous boundaryconditions only for simplifying numerical analysis.
Remark 2.2.
The scheme (2.1) depends on the definition of ( G h v h )( p ) at eachvertex p ∈ N h . In the following, three popular definitions of ( G h v h )( p ) are listed(c.f., [34, 47]).(a) Weighted averaging(WA). For each p ∈ N h , let the element patch ω p = ∪{ τ : p ∈ ¯ τ } and define ( G h v h )( p ) = 1 | ω p | (cid:90) ω p ∇ v h ( x , x ) dx dx . (2.3) (b) Local L -projection. We seek two polynomials P l ∈ P ( ω p ) , ( l = 1 , , such that (cid:90) ω p [ P l ( x , x ) − ∂ x l v ( x , x )] Q ( x , x ) dx dx = 0 , ∀ Q ∈ P ( ω p ) , l = 1 , nd we define ( G h v h )( p ) = ( P ( p ) , P ( p )) . Sometimes, the exact integral in (2.4) is replaced by its discrete counterpart so thatthe two polynomials P l , l = 1 , satisfying the least square fitting equation (SPR) m (cid:88) i =1 [ P l ( x i , x i ) − ∂ x il v ( x i , x i )] Q ( x i , x i ) = 0 , ∀ Q ∈ P ( ω P ) , l = 1 , , (2.5) where ( x i , x i ) , i = 1 , . . . , m are m given points in ω p .(c) The polynomial preserving recovery (PPR). We seek a quadratic function P ∈P ( ω p ) , such that m (cid:88) i =1 [ P ( x i , x i ) − v ( x i , x i )] Q ( x i , x i ) = 0 , ∀ Q ∈ P ( ω P ) . (2.6) Then we can define ( G h v h )( p ) = ( ∂ x P ( p ) , ∂ x P ( p )) .It is known that the above three definitions are equivalent on a mesh of uniformtriangular pattern [45]. Remark 2.3.
Essentially, the operator G h can be regarded as a difference opera-tor defined on unstructured grids. This operator lifts discontinuous gradient generatedfrom a C -FEM to a continuous one, and thereby makes the further calculation ofhigh order derivatives possible. Remark 2.4.
The scheme (2.1) is very simple and straightforward. It avoidsthe complicated construction of conforming C finite element basis (c.f., [28]) or thecomplicated construction of nonconforming penalty terms ( [27]). For
A ⊂
Ω, let V h ( A ) denote the restrictions of functions in V h to A and let V comp h ( A ) denote the set of those functions in V h ( A ) with compact support in theinterior of A [37]. Let Ω ⊂⊂ Ω ⊂⊂ Ω ⊂⊂ Ω be separated by d ≥ c o h and (cid:96) be adirection, i.e., a unit vector in R . Let τ be a parameter, which will typically be amultiply of h . Let T (cid:96)τ denote translation by τ in the direction (cid:96) , i.e., T (cid:96)τ v ( x ) = v ( x + τ (cid:96) ) , (2.7)and for an integer ν T (cid:96)ντ v ( x ) = v ( x + ντ (cid:96) ) . (2.8)Following the definition of [37], the finite element space V h is called translation in-variant by τ in the direction (cid:96) if T (cid:96)ντ v ∈ V comp h (Ω) , ∀ v ∈ V comp h (Ω ) , (2.9)for some integer ν with | ν | < M . Equivalently, T h is called a translation invariantmesh. As illustrated in [25], uniform meshes of regular pattern, chevron pattern,cirsscross patter, and unionjack pattern are all translation invariant.
3. Analysis.
The section is dedicated to a mathematical proof for the conver-gence properties.To this end, we need some properties of G h . For the polynomial preservingrecovery operator G h , there are the following boundedness property (see (2.11) in [34]) (cid:107) G h v h (cid:107) (cid:46) | v h | , v h ∈ V h (3.1) nd the superconvergence approximation properties (cid:107)∇ u − G h u I (cid:107) (cid:46) h | u | , ∞ , u ∈ W , ∞ (Ω) . (3.2)Here u I is the linear interpolation of u in V h . In addition, we will utilize the followingultraconvergence approximation properties of Hessian recovery operator (see Theorem3.5 in [25]) (cid:107) D u − DG h u I (cid:107) (cid:46) h | u | , ∞ , u ∈ W , ∞ (Ω) , (3.3) (cid:107) D u − G h u I (cid:107) (cid:46) h | u | , ∞ , u ∈ W , ∞ (Ω) , (3.4)provided the mesh T h is translation invariant. Remark 3.1.
We would like to comment that the requirement T h translationinvariant for proving approximation properties (3.3) and (3.4) only for theoreticalpurpose. In practical, our method can be applied to and shows optimal convergenceon arbitrary unstructured mesh. To analyze the convergence of the scheme (2.1), we suppose that the T h is sufficientregular such that there holds the following discrete Poincar´e inequality (cf., [26]) (cid:107) v h (cid:107) i (cid:46) (cid:107) G h v h (cid:107) i , ∀ v h ∈ V h , i = 0 , weak approximation properties (cid:12)(cid:12)(cid:12)(cid:12)(cid:90) Ω ∇ v · ( G h v h − ∇ v h ) (cid:12)(cid:12)(cid:12)(cid:12) (cid:46) h (cid:107) v (cid:107) | G h v h | , ∀ v ∈ H , (3.6) (cid:12)(cid:12)(cid:12)(cid:12)(cid:90) Ω ∇ v · ( G h v h − ∇ v h ) (cid:12)(cid:12)(cid:12)(cid:12) (cid:46) h (cid:107) v (cid:107) | G h v h | , ∀ v ∈ H . (3.7)Note that both (3.5) and (3.6) have been discussed in the analysis of a recoveryoperator based linear finite element method for the biharmonic equation by Guo et al.in [26]. In their paper, a counter-example shows that the strong error (cid:107)∇ v h − G h v h (cid:107) is not necessary of O ( h ) for all v h ∈ V h (which means the weak estimate (3.6) mightbe the best estimate of the difference ∇ v h − G h v h ).By (3.5), for all v h ∈ V h , we have (cid:107) v h (cid:107) (cid:46) (cid:107) G h v h (cid:107) (cid:46) (cid:107) G h v h (cid:107) (cid:46) (cid:107) DG h v h (cid:107) . In other words, the semi-norm (cid:107) DG h ·(cid:107) is a norm. Then by the Lax-Milgram theorem,the scheme (2.1) has a unique solution . Moreover, by (2.1), (cid:107) DG h u h (cid:107) = a h ( u h , u h ) = ( f, u h ) (cid:46) (cid:107) f (cid:107) (cid:107) v h (cid:107) . Then (cid:107) DG h u h (cid:107) (cid:46) (cid:107) f (cid:107) (3.8)which implies the stability of our scheme. H error estimate. Theorem 3.1.
Let u h be the solution of (2.1) and u ∈ H the solution of (1.3) . If the mesh T h is translation invariant, G h is properlydefined such that (3.1) - (3.7) hold, then (cid:107) DG h ( u h − u I ) (cid:107) (cid:46) h (cid:107) u (cid:107) , (3.9) here u I is the linear interpolation of u in V h . Consequently, (cid:107) D u − DG h u h (cid:107) (cid:46) h (cid:107) u (cid:107) . (3.10) Proof . Since a weak solution of (1.3) which has regularity u ∈ H is also thestrong solution satisfying (1.1) and u h is a discrete solution satisfying (2.1), we have a h ( u h , v h ) = ( f, v h ) = ( − ∆ u, v h ) , ∀ v h ∈ V h . Using the fact that v h ∈ V h , we have( − ∆ u, v h ) = ( − ( ∇ · ∇ )(∆ u ) , v h ) = ( ∇ (∆ u ) , ∇ v h ) . Therefore, a h ( u h , v h ) = I + ( ∇ (∆ u ) , G h v h ) , with I = ( ∇ (∆ u ) , ∇ v h − G h v h ) . (3.11)Now we deal with the term ( ∇ (∆ u ) , G h v h ). Since G h v h · n = G h v h · t = 0 onthe boundary ∂ Ω, we have that on ∂ Ω, ∂ ∆ u∂ n · G h v h = ∂ ∆ u∂ n ∂ t ( G h v h · t ) + ∂ ∆ u∂ n ( G h v h · n ) = 0 . Then ( ∇ (∆ u ) , G h v h ) = (( ∇ · ∇ )(∆( ∇ u )) , G h v h )= − ( D (∆ u ) , DG h v h ) := − (cid:90) Ω D (∆ u ) : DG h v h . Consequently, ( ∇ (∆ u ) , G h v h ) = I − ( D (∆ u ) , G h v h ) , with I = ( D (∆ u ) , G h v h − DG h v h ) . (3.12)Finally, we deal with the term − ( D (∆ u ) , G h v h ). Writing the gradient as ∇ u = ∂u∂ n n + ∂u∂ t t , we have D u = ∂ u∂ n nn T + ∂ u∂ t tt T + ∂ u∂ n ∂ t ( nt T + tn T ) , and consequently, ∂D u∂ n = ∂ u∂ n nn T + ∂ u∂ n ∂ t tt T + ∂ u∂ n ∂ t ( nt T + tn T ) . oting that u = ∂u∂ n = ∂ u∂ n = 0 on ∂ Ω, we have ∂ u∂ n ∂ t = ∂ u∂ n ∂ t = 0 on Ω . Therefore, on ∂ Ω, ∂D u∂ n : G h v h = ∂ u∂ n nn T : G h v h = ∂ u∂ n n T G h v h n = 0 , where in the last equality, we used the fact that n T G h v h n = 0. By Green’s formulas,we finally obtain − ( D (∆ u ) , G h v h ) = − (∆( D u ) , G h v h ) = ( D u, DG h v h ) . In summary, by letting I = ( D u − DG h u I , DG h v h ) , we obtain a h ( u h − u I , v h ) = I + I + I . (3.13)Next, we estimate I i (i = 1, 2, 3) term by term. By (3.5) and (3.6), we have | I | (cid:46) h (cid:107) u (cid:107) | G h v h | (cid:46) h || u || || DG h v h || . Similarly, by (3.3), | I | (cid:46) h | u | , ∞ || DG h v h || (cid:46) h (cid:107) u (cid:107) || DG h v h || . On the other hand, (3.7) implies | I | (cid:46) h (cid:107) u (cid:107) || DG h v h || . In a conclusion, we obtain that | a h ( u h − u I , v h ) | (cid:46) h (cid:107) u (cid:107) (cid:107) DG h v h (cid:107) , ∀ v h ∈ V h . Choosing v h = u h − u I in the above estimate, we have (3.9).The estimate (3.10) is a direct consequence of (3.9) and (3.3), together with thetriangle inequality. H error estimate. In this section, we use the Aubin-Nitsche technique toestimate the H norm error (cid:107)∇ u − G h u h (cid:107) . To this end, we construct the followingauxiliary problems :1) Find U ∈ H (Ω) such that (cid:90) Ω D U : D v = ( G h ( u h − u I ) , ∇ v ) , ∀ v ∈ H (Ω) . (3.14)2) Find U h ∈ V h such that a h ( U h , v h ) = ( G h ( u h − u I ) , ∇ v h ) , ∀ v h ∈ V h . (3.15) t is easy to deduce from (3.14) that (cid:107) U (cid:107) (cid:46) (cid:107) G h ( u h − u I ) (cid:107) , (3.16)and from (3.15),(3.5) that (cid:107) D ( G h U h ) (cid:107) (cid:46) (cid:107) G h ( u h − u I ) (cid:107) . (3.17) Theorem 3.2.
Let u h be the solution of (2.1) and u ∈ H the solution of (1.3) .If the mesh T h is translation invariant, and G h is properly defined such that (3.1) - (3.7) hold, then (cid:107) G h ( u h − u I ) (cid:107) (cid:46) h (cid:107) u (cid:107) . (3.18) Consequently, (cid:107)∇ u − G h u h (cid:107) (cid:46) h (cid:107) u (cid:107) . (3.19) Proof . First, by the definition of the auxiliary problems, we have( G h ( u h − u I ) , ∇ ( u h − u I )) = a h ( U h , u h − u I )= a h ( u h − u I , U h )= ( f, U h ) − a h ( u I , U h )= ( −(cid:52) u, U h ) − a h ( u I , U h ) . Using the same splitting techniques in the previous theorem, we can write( G h ( u h − u I ) , ∇ ( u h − u I )) = J + J + J + J , where J = ( ∇ ( (cid:52) u ) , ∇ U h − G h U h ) ,J = ( D ( (cid:52) u ) , G h U h − DG h U h ) ,J = ( D u − DG h u I , D U ) ,J = ( D u − DG h u I , DG h U h − D U ) . We first estimate J and J . By (3.7), | J | (cid:46) h (cid:107) u (cid:107) (cid:107) DG h U h (cid:107) (cid:46) h (cid:107) u (cid:107) (cid:107) DG h U h (cid:107) (cid:46) h (cid:107) u (cid:107) (cid:107) G h ( u h − u I ) (cid:107) . and | J | (cid:46) h (cid:107) u (cid:107) (cid:107) DG h U h (cid:107) (cid:46) h (cid:107) u (cid:107) (cid:107) G h ( u h − u I ) (cid:107) . Moreover, using the integration by parts, | J | ≤ (cid:107) D u − G h u I (cid:107) (cid:107) D U (cid:107) (cid:46) h (cid:107) u (cid:107) (cid:107) G h ( u h − u I ) (cid:107) . inally, by Cauchy-Schwartz inequality, | J | (cid:46) h (cid:107) u (cid:107) (cid:107) U (cid:107) (cid:46) h (cid:107) u (cid:107) (cid:107) G h ( u h − u I ) (cid:107) (cid:46) h (cid:107) u (cid:107) . Summarizing all the above estimates, we obtain( G h ( u h − u I ) , ∇ ( u h − u I )) (cid:46) h (cid:107) u (cid:107) (cid:107) G h ( u h − u I ) (cid:107) + h (cid:107) u (cid:107) . Noticing that (cid:107) G h ( u h − u I ) (cid:107) ∼ ( G h ( u h − u I ) , ∇ ( u h − u I )) , we arrive that (cid:107) G h ( u h − u I ) (cid:107) (cid:46) h (cid:107) u (cid:107) (cid:107) G h ( u h − u I ) (cid:107) + h (cid:107) u (cid:107) . Then the estimate (3.18) follows.The H error estimate of (3.19) is a direct consequence of (3.18) and (3.2). L error estimate. Theorem 3.3.
Let u h be the solution of (2.1) and u ∈ H the solution of (1.3) . If the mesh T h is sufficiently regular and uniform, G h is properly defined such that (3.1) - (3.7) hold, then (cid:107) u − u h (cid:107) (cid:46) h (cid:107) u (cid:107) . (3.20) Proof . By (3.5), (cid:107) u I − u h (cid:107) (cid:46) (cid:107) G h ( u I − u h ) (cid:107) (cid:46) h (cid:107) u (cid:107) . Then (cid:107) u − u h (cid:107) (cid:46) (cid:107) u − u I (cid:107) + (cid:107) u I − u h (cid:107) (cid:46) h (cid:107) u (cid:107) . Remark 3.2.
In the second section, we observed the convergence rates O ( h ) both for the errors (cid:107) u − u h (cid:107) and (cid:107)∇ u − G h u h (cid:107) . However, we can only prove theorder O ( h ) from our analysis. Further analysis to the scheme is desired to prove theoptimal convergence rates of (cid:107) u − u h (cid:107) and (cid:107)∇ u − G h u h (cid:107) .
4. Numerical Experiments.
In this section, we present several numerical ex-periments to show the convergence rates and efficiency of our method. In all ournumerical experiments, G h is chosen as the polynomial preserving recovery opera-tor [46]. To present our numerical results, the following notations are used : De := (cid:107) u − u h (cid:107) , D e := (cid:107)∇ u − ∇ u h (cid:107) ,D r e := (cid:107)∇ u − G h u h (cid:107) , D e := (cid:107) D u − DG h u h (cid:107) ,D e := (cid:107) D u − DG h u h (cid:107) . Moreover, the convergence rates are listed with respect to the degree of freedom(Dof).Noticing Dof ≈ h − for a two dimensional grid, the corresponding convergent rateswith respect to the mesh size h are double of what we present in the tables 3.1-3.11. xample 1 . We consider the triharmonic problem (cid:26) − ∆ u = f in Ω = (0 , × (0 , u = ∂ n u = ∂ nn u = 0 on ∂ Ω , (4.1)where f is chosen to fit the exact solution u ( x , x ) = x (1 − x ) x (1 − x ) .First, we apply our scheme (2.1) on regular pattern uniform triangular mesh.The corresponding numerical results are listed in Table 4.1. It shows that numericalsolution u h converges to the exact solution u at a rate of O ( h ) in recovered H norm.Also from this table, we observe that De and D r e converge at a rate of O ( h ) while D e and D e converge at a rate of O ( h ). Note that convergence rate of De and D r e is better than that proved in Theorems 3.2 and 3.3. Table 4.1
Numerical Results Of Example 1 On Regular Pattern Uniform Mesh
Dof De order D e order D r e order D e order D e order1089 5.61e-06 – 5.76e-05 – 2.57e-05 – 4.01e-04 – 4.46e-03 –4225 1.54e-06 0.95 2.20e-05 0.71 7.03e-06 0.96 1.73e-04 0.62 2.06e-03 0.5716641 3.99e-07 0.99 9.65e-06 0.60 1.83e-06 0.98 8.18e-05 0.54 9.99e-04 0.5366049 1.01e-07 0.99 4.62e-06 0.53 4.66e-07 0.99 4.03e-05 0.51 4.93e-04 0.51 Secondly, we test our scheme on uniform triangular meshes of other patterns,including the chevron, Criss-cross, and Union-Jack patterns. Numerical data arelisted in 4.2, Table 4.3, and Table 4.4, respectively. Again, we observed O ( h ) for De , O ( h ) for D e , O ( h ) for D r e , O ( h ) for D e , and O ( h ) for D e , the same as theregular pattern. Table 4.2
Numerical Results Of Example 1 On Chevron Pattern Uniform Mesh
Dof De order D e order D r e order D e order D e order1089 4.48e-06 – 6.00e-05 – 2.02e-05 – 3.81e-04 – 4.26e-03 –4225 1.25e-06 0.94 2.22e-05 0.73 5.65e-06 0.94 1.69e-04 0.60 2.04e-03 0.5416641 3.24e-07 0.99 9.67e-06 0.61 1.47e-06 0.98 8.14e-05 0.53 9.97e-04 0.5266049 8.24e-08 0.99 4.62e-06 0.54 3.75e-07 0.99 4.03e-05 0.51 4.92e-04 0.51 Table 4.3
Numerical Results Of Example 1 On Criss-cross Pattern Uniform Mesh
Dof De order D e order D r e order D e order D e order2113 1.50e-05 – 2.91e-03 – 3.13e-05 – 3.61e-03 – 7.21e-03 –8321 3.83e-06 1.00 1.48e-03 0.50 8.02e-06 0.99 1.83e-03 0.50 3.62e-03 0.5033025 9.46e-07 1.01 7.20e-04 0.52 2.02e-06 1.00 8.98e-04 0.52 1.81e-03 0.50131585 2.39e-07 0.99 3.64e-04 0.49 5.09e-07 1.00 4.55e-04 0.49 9.07e-04 0.50 Finally, we turn to the Delaunay mesh. The first level coarse mesh is generatedby EasyMesh [23] followed by three levels of regular refinement. Table 4.5 presentsthe convergence history for the five different errors. O ( h ) and O ( h ) convergencerates are observed for L and H errors. As for the L error of recovered gradient, O ( h ) superconvergence is observed. Regarding recovered H and H errors, O ( h )convergence are observed .In summary, we see that our method converges with optimal rates on all fourtested uniform meshes as well as the Delaunay mesh. able 4.4 Numerical Results Of Example 1 On Unionjack Pattern Uniform Mesh
Dof De order D e order D r e order D e order D e order1089 2.61e-05 – 3.24e-03 – 6.30e-05 – 4.33e-03 – 9.40e-03 –4225 6.53e-06 1.02 1.61e-03 0.51 1.59e-05 1.02 2.12e-03 0.53 4.52e-03 0.5416641 1.63e-06 1.01 8.03e-04 0.51 4.00e-06 1.01 1.06e-03 0.51 2.22e-03 0.5266049 4.08e-07 1.01 4.01e-04 0.50 1.00e-06 1.00 5.26e-04 0.50 1.10e-03 0.51 To show the efficiency of our method, we make some numerical comparison withthe cubic C interior penalty method [27] on the same Delaunay meshes. Table 4.6shows numerical results of the C interior penalty method in the L norm and theenergy norm defined in [27]. Consisting with the theoretical result established in [27],the error in the energy norm converges linearly and the L error decays at rate O ( h ).Figures 4.1 and 4.2 depict convergent rates of these two methods(i.e. our methodand the C interior penalty method) under the discrete H (the energy) and L norms.The rates are almost the same. However, to achieve the same accuracy, our algorithmuses about one-eighth degrees of freedom of the C interior penalty method. Table 4.5
Numerical Results of Example 1 on Delaunay Triangulation with Regular Refinement
Dof De order D e order D r e order D e order D e order513 1.08e-05 – 9.72e-05 – 4.93e-05 – 6.21e-04 – 6.34e-03 –1969 3.02e-06 0.95 3.04e-05 0.86 1.38e-05 0.95 2.41e-04 0.70 2.90e-03 0.587713 7.94e-07 0.98 1.23e-05 0.66 3.65e-06 0.97 1.08e-04 0.59 1.37e-03 0.5530529 2.03e-07 0.99 5.77e-06 0.55 9.35e-07 0.99 5.20e-05 0.53 6.67e-04 0.52 Table 4.6 C Interior Penalty Method for Example 1 on Delaunay Triangulation with Regular Refinement
Dof De order D h e order4369 7.75e-06 – 6.25e-03 –17233 2.72e-06 0.76 2.92e-03 0.5568449 7.89e-07 0.90 1.37e-03 0.55272833 1.58e-07 1.16 6.63e-04 0.52 Example 2.
In the second example, we will show that our scheme works wellalso for problems nonhomogeneous boundary conditions. We consider the equation − ∆ u = sin(2 πx ) cos(2 πx ) , ( x , x ) ∈ [0 , whose exact solution is u ( x , x ) = 1512 π sin(2 πx ) cos(2 πx ) , It provides nonhomogeneous boundary conditions u | ∂ Ω , ∂ n u | ∂ Ω , ∂ nn u | ∂ Ω .As in Example 1, we first test our algorithm on regular pattern uniform triangularmesh and list the numerical results in Table 4.7. Again, D e decays at rate O ( h ).Asexpected, both De and D r e converges with order O ( h ). Astonishingly, both D e and D e also converge quadratically. Namely, for this example, both D e and D e superconverge. umber of DOF E rr o r -4 -3 -2 ppr based FEMC0IP Fig. 4.1 . Comparison of Discrete H Errorsfor Example 1
Number of DOF E rr o r -7 -6 -5 -4 ppr based FEMC0IP Fig. 4.2 . Comparison of Discrete L Errorsfor Example 1
Table 4.7
Numerical Results Of Example 2 On Regular Pattern Uniform Mesh
Dof De order D e order D r e order D e order D e order1089 2.50e-07 – 3.21e-05 – 1.53e-06 – 1.35e-04 – 2.25e-03 –4225 2.66e-08 1.65 7.09e-06 1.11 1.47e-07 1.73 2.93e-05 1.13 7.65e-04 0.8016641 3.01e-09 1.59 1.46e-06 1.15 1.92e-08 1.48 6.75e-06 1.07 3.24e-04 0.6366049 4.75e-10 1.34 3.32e-07 1.08 3.41e-09 1.25 1.76e-06 0.98 1.37e-04 0.62 We then consider chevron pattern uniform triangular mesh. Table 4.8 clearlyindicates that u h converges to u at a rate of O ( h ) under the L norm, at a rate of O ( h ) under the H norm and the recovered H and H norms. Moreover, the recoverygradient G h u h converges to ∇ u at a rate of O ( h ). We also test our algorithms onDelaunay meshes as in the previous example. The numerical data are demonstratedin Table 4.9. Similar to what we observed in Chevron pattern uniform triangularmesh, the computed error by our method converges to 0 with optimal rates undervarious norms.In addition, we have tested our algorithms on other two types (Criss-cross andUnion-Jack pattern) uniform triangular meshes. Since the numerical results are sim-ilar to the corresponding parts in the previous example, they are not reported here. Table 4.8
Numerical Results of Example 2 on Chevron Pattern Uniform Mesh
Dof De order D e order D r e order D e order D e order1089 4.39e-08 – 1.08e-06 – 4.20e-07 – 1.23e-05 – 2.70e-04 –4225 3.38e-09 1.89 4.47e-07 0.65 3.92e-08 1.75 3.65e-06 0.90 9.95e-05 0.7416641 7.89e-10 1.06 2.22e-07 0.51 7.57e-09 1.20 1.61e-06 0.60 4.13e-05 0.6466049 2.00e-10 1.00 1.11e-07 0.50 1.83e-09 1.03 7.84e-07 0.52 1.78e-05 0.61 Once again, we present a numerical comparison with the C interior penaltymethod. We see from Figure 4.3 that the convergence rates of H error are compa-rable, however, our method requires much less degrees of freedom in order to achievethe same accuracy. Figure 4.4 indicates that our method is slightly better than the C interior penalty method with regard to the L norm which is suboptimal. Here we able 4.9 Numerical Results of Example 2 on Delaunay Triangulation with Regular Refinement
Dof De order D e order D r e order D e order D e order513 5.52e-08 – 1.77e-06 – 3.81e-07 – 1.29e-05 – 1.44e-04 –1969 1.17e-08 1.15 5.43e-07 0.88 8.07e-08 1.15 4.21e-06 0.83 5.78e-05 0.687713 2.79e-09 1.05 2.57e-07 0.55 1.97e-08 1.03 1.90e-06 0.58 1.94e-05 0.8030529 6.86e-10 1.02 1.27e-07 0.51 4.91e-09 1.01 9.40e-07 0.51 8.52e-06 0.60 would like to point out that the error of the C interior penalty method is sensitiveto the penalty parameter. Number of DOF E rr o r -6 -5 -4 -3 ppr based FEMC0IP Fig. 4.3 . Comparison of Discrete H Errorsfor Example 2
Number of DOF E rr o r -10 -9 -8 -7 ppr based FEMC0IP Fig. 4.4 . Comparison of Discrete L Errorsfor Example 2
Example 3.
In previous two examples, we consider sixth order elliptic equationson the unit square. To show the ability of dealing arbitrary complex domain, weconsider the following sixth order partial differential equation − ∆ u = 8 e x + x . on the unit disk, i.e. Ω = { ( x , x ) ∈ R : x + x ≤ } . The exact solution is u ( x , x ) = e x + x . and the corresponding boundary conditions are given by the exact solution. Theinitial mesh is generated by DistMesh [36] as shown in Figure 4.5. The other sevenlevel meshes are obtained by refining the initial mesh using regular refinement. Thenumerical results are reported in Table 4.10. As in two previous examples, O ( h )convergence for D e , D e , and D e are observed and O ( h ) convergence order areobserved for De and D r e . Example 4.
As in [28], we consider the following triharmonic equation − ∆ u = 0 . on the L-shaped domain [ − , \ ([0 , × [ − , u ( x , x ) = x − x . Fig. 4.5 . Initial Mesh On The Unit Disk -1 -0.5 0 0.5 1-1-0.8-0.6-0.4-0.200.20.40.60.81
Fig. 4.6 . Initial Mesh On The Lshape Domain
Table 4.10
Numerical Results of Sixth Order PDE in the Unit Disk
Dof De order D e order D r e order D e order D e order88 3.69e-02 – 4.76e-01 – 1.38e-01 – 1.01e+00 – 3.70e+00 –318 1.44e-02 0.73 1.91e-01 0.71 3.21e-02 1.14 3.48e-01 0.83 1.92e+00 0.511207 1.84e-03 1.54 8.05e-02 0.65 4.76e-03 1.43 1.06e-01 0.89 7.14e-01 0.744701 4.16e-04 1.09 4.00e-02 0.51 1.20e-03 1.01 4.94e-02 0.56 3.66e-01 0.4918553 1.04e-04 1.01 2.00e-02 0.51 3.03e-04 1.00 2.40e-02 0.52 1.99e-01 0.4473713 2.60e-05 1.00 1.00e-02 0.50 7.57e-05 1.01 1.19e-02 0.51 9.70e-02 0.52293857 7.01e-06 0.95 5.00e-03 0.50 1.85e-05 1.02 5.94e-03 0.50 4.62e-02 0.54 Here we use uniform meshes. The initial mesh is plotted in Figure 4.6, while ournumerical results are listed in Table 4.11. As pointed out in [28], the solution u variesfast near the boundary. Even in that case, we observe the optimal convergence ratesunder all the norms. Table 4.11
Numerical Results For Triharmonic Equation On LShape Domain
Dof De order D e order D r e order D e order D e order225 2.04e-01 – 5.78e+00 – 8.10e-01 – 1.42e+01 – 7.58e+01 –833 2.63e-02 1.56 1.69e+00 0.94 1.24e-01 1.43 4.12e+00 0.95 2.96e+01 0.723201 3.54e-03 1.49 4.20e-01 1.03 2.70e-02 1.14 1.42e+00 0.79 1.35e+01 0.5812545 6.66e-04 1.22 1.39e-01 0.81 6.67e-03 1.02 5.71e-01 0.66 6.70e+00 0.5249665 1.44e-04 1.12 5.95e-02 0.62 1.71e-03 0.99 2.62e-01 0.57 3.35e+00 0.50197633 3.35e-05 1.05 2.82e-02 0.54 4.55e-04 0.96 1.27e-01 0.52 1.68e+00 0.50 In summary, our numerical experiments discover that our algorithm convergeswith optimal rates under various norms, for sixth order equations on different kindsof domains, with homogenous or nonhomogeneous boundary conditions. In addition,comparing to some existed algorithm such as C interior penalty method, our algo-rithm has much lower computational cost.
5. Concluding remarks.
In this work, we developed a PPR based discretiza-tion algorithm for a sixth-order PDE. The algorithm has a simple form and is easyto implement. Moreover, it has optimal convergence rates as the existing conforming nd nonconforming FEMs in the literatures for sixth-order PDEs. However, the newmethod seems to be more advantageous with respect to computational complexity.Generally speaking, the recovery operator is a special difference operator onnonuniform grids. It can be used to compute high order derivatives of a functionwhich are piecewise polynomials but only globally in C and thus can be used todiscretize PDEs of higher order. On the other hand, how to choose different recov-ery operators for different PDEs deserves more in-depth mathematical study. Furtherinvestigation is called for to find simple and efficient algorithms for complicated PDEs. REFERENCES[1] A. Adini and R.W. Clough, Analysis of plate bending by the finite element method,
NSF reportG. 7337 , 1961.[2] M. Ainsworth J.T. Oden,
A Posteriori Error Estimation in Finite Element Analysis , WileyInterscience, New York, 2000.[3] I. Babuska and T. Strouboulis,
The Finite Element Method and Its Reliability , Oxford Univer-sity Press, London, 2001.[4] R. E. Bank and A. Weiser, Some a posteriori error estimators for elliptic partial differentialequations.
Math. comp. , 44(1985), 283-301.[5] R. E. Bank and J. Xu, Asymptotically exact a posteriori error estimators, Part I: Grid withsuperconvergence,
SIAM J. Numer. Anal. , 41(2003), 2294-2312.[6] G. A. Baker, Fintie element methods for elliptic equations using nonconforming elements,
Math.Comp. , 31(1977), 45-59.[7] R. E. Bank and J. Xu, Asymptotically exact a posteriori error estimators, Part II: Generalunstructured grids,
SIAM J. Numer. Anal. , 41(2003), 2313-2332.[8] R. Backofen, A. R¨atz, and A. Voigt, Nucleation and growth by a phase crystal (PFC) model,Phil. Mag. Lett., 87(2007):813C-820.[9] J. W. Barrett, S. Langdon, Stephen and R. N¨urnberg, Finite element approximation of a sixthorder nonlinear degenerate parabolic equation, Numer. Math., 96(2004), 401–434.[10] D. Biskamp, E. Schwarz, and J.F. Drake. Ion-controlled collisionless magnetic reconnection.Physical Review Letters, 75(1995):3850–3853.[11] D. Biskamp, E. Schwarz, A. Zeiler, A. Celani, and J.F. Drake. Electron magnetohydrodynamicturbulence. Physics of Plasmas, 6(1999):751-758.[12] S. Brenner and L. Sung, C0 interior penalty methods for fourth order elliptic boundary valueproblems on polygonal domains,
J. Sci. Comput. , 22/23(2005), 83–118.[13] S. Brenner and L.R. Scott,
Mathematical Theory of Finite element Methods , 3rd edition,Spriger-Verlag, New York, 2008.[14] J. W. Cahn, On spinodal decomposition, Acta Metall, 9 (1961), 795–801.[15] G. Caginalp and P. Fife. Higher-order phase field models and detailed anisotropy. PhysicalReview B, 34(1986):4940–4943, 1986.[16] A. S. Chang and W. Chen. A note on a class of higher order comformally covariant equations,Discrete and Continuous Dymanical Systems, 7 (2001), 275C281.[17] H. Chen, H. Guo, Z. Zhang, and Q. Zou, A C0 linear finite element method for two fourth-ordereigenvalue problems. IMA J. Numer. Anal. 37 (2017), no. 4, 2120–2138.[18] M. Cheng and J. A. Warren. An efficient algorithm for solving the phase field crystal model, J.Comput. Phys., 227(2008):6241C6248.[19] P.G. Ciarlet,
The Finite Element Method for Elliptic Problems , Studies in Mathematics andits Applications, Vol.4, North-Holland, Amsterdam, 1978.[20] Introduction to COMSOL Multiphysics Version 5.1, March 2015, page 46.[21] Francoise Chatelin,
Spectral Approximation of Linear Operators , Computer Science and Ap-plied Mathematics, Academic Press Inc., New York, 1983.[22] J.F. Drake, D. Biskamp, and A. Zeiler. Breakup of the electron current layer during 3-d colli-sionless magnetic reconnection. Geophysical Research Letters, 24(1997):2921– 2924.[23]
B. Niceno , EasyMesh Version 1.4: A Two-Dimensional Quality Mesh Generator , .[24] W.I. Fushchych and Z.I. Symenoh, High-order equations of motion in quantum mechanics andgalilean relativity. Journal of Physics A: Mathematical and General, 30(1997):131–135,1997.[25] H. Guo, Z. Zhang, and R. Zhao, Hessian Recovery for Finite Element Methods, Math. Comp.166 (2017), no. 306, 1671–1692.[26] H. Guo, Z. Zhang, and Q. Zou, A C linear finite element method for biharmonic problemsbased on gradient recovery, J. Sci. Comput.
74 (2018), no. 3, 1397–1422.[27] T. Gudi and M. Neilan, An interior penalty method for a sixth-order elliptic equation,
IMA J.Numer. Anal. , 31 (2011), 1734–1753.[28] J. Hu, and S. Zhang, The minimal conforming H k finite element spaces on R n rectangulargrids, Math. Comp. , 84 (2015), 563–579.[29] J.E. Hilliard and J.W. Cahn, Free energy of a non-uniform system. I. Interfacial free energy, J.Chem. Phys., 28 (1958), pp. 258C267.[30] J.E. Hilliard and J.W. Cahn, Free energy of a non-uniform system. III. Nucleation in a twocomponent incompressible fluid, J. Chem. Phys., 31 (1959), pp. 688C699.[31] Mohamed El-Gamel and Mona Sameeh, An efficient technique for finding the eigenvalues offourth-order Sturm-Liouville problems, Applied Mathematics, 3 (2012), 920–925.[32] L. Morley, The triangular equilibrium problem in the solution of plate bending problems.
Aero.Quart. , 19 (1968), 149C-169.[33] A. Naga and Z. Zhang, A posteriori error estimates based on the polynomial preserving recovery,SIAM J. Numer. Anal., 42-4 (2004), 1780–1800.[34] A. Naga and Z. Zhang, The polynomial-preserving recovery for higher order finite elementmethods in 2D and 3D, Discrete and Continuous Dynamical Systems-Series B, 5-3 (2005),769–798.[35] A. Naga and Z. Zhang, Function value recovery and its application in eigenvalue problems,SIAM J. Numer. Anal., 50 (2012), 272–286.[36] P.-O. Persson and G. Strang, A simple mesh generator in Matlab, SIAM Rev. 46 (2004), 329–345.[37] L.B. Wahlbin, Superconvergence in Galerkin finite element methods, Lecture Notes in Mathe-matics, Springer-Verlag, Berkin, 1995.[38] M. Wang and J. Xu, The Morley element for fourth order elliptic equations in any dimensions.Numer. Math. 103 (2006), 155–169.[39] S.M. Wise, J.S. Lowengrub, J.S. Kim, and W.C. Johnson. Efficient phase-field simulation ofquantum dot formation in a strained heteroepitaxial film. Superlattices and Microstruc-tures, 36(2004) : 293–304.[40] S. M. Wise, C. Wang, and J. S. Lowengrub. An energy-stable and convergent finite differencescheme for the phase field crystal equation, SIAM J. Numer. Anal, 47(2009):2269C2288.[41] H. Ugail. Partial Differential Equations for Geometric Design, Springer, NewYork, 2011.[42] E. Ventsel and T. Krauthammer. Thin Plates & Shells: Theory, Analysis, & Applications.CRC, 2001.[43] J. Xu and Z. Zhang, Analysis of recovery type a posteriori error estimators for mildly structuredgrids, Math. Comp., 73 (2004), 1139–1152.[44] S. Zhang and Z. Zhang, Invalidity of decoupling a biharmonic equation to two Poisson equationson non-convex polygons.
Int. J. Numer. Anal. Model.