A Priori and A Posteriori Error Control of Discontinuous Galerkin Finite Element Methods for the Von Kármán Equations
Carsten Carstensen Gouranga Mallik Neela Nataraj

	
Semantic Foggy Scene Understanding with Synthetic Data



    Christos Sakaridis Dengxin Dai Luc Van Gool 

    
=========================================================



 
This paper analyses discontinuous Galerkin finite element methods (DGFEM) to approximate a regular solution to the von Kármán equations defined on a polygonal
domain. A discrete inf-sup condition sufficient for the stability of the discontinuous
Galerkin discretization of a well-posed linear problem is established and this allows
the proof of local existence and uniqueness of a discrete solution to the non-linear problem
with a Banach fixed point theorem. The Newton scheme is locally second-order convergent and appears to be  a robust solution strategy up to machine precision. A comprehensive a priori and a posteriori energy-norm error analysis relies on one sufficiently large stabilization parameter and sufficiently fine triangulations. In case the other stabilization parameter
degenerates towards infinity, the DGFEM reduces to a novel C^0 interior penalty method (IPDG). In contrast to the known C^0-IPDG due to Brenner et al <cit.>, the overall
discrete formulation maintains symmetry of the trilinear form in the first two components – despite the general non-symmetry of the discrete nonlinear problems. Moreover, a reliable and efficient a posteriori error analysis immediately follows for the DGFEM of this paper, while the
different norms in the known C^0-IPDG lead to complications with some non-residual type remaining terms. Numerical experiments confirm the best-approximation results and the
equivalence of the error and the error estimators. A related adaptive mesh-refining algorithm leads to optimal empirical convergence rates for a non convex domain.





§ INTRODUCTION

The discontinuous Galerkin finite element methods (DGFEM) have become popular for the numerical solution of a large range of problems in partial differential equations, which include linear and nonlinear problems, convected-dominated diffusion for second- and fourth-order elliptic problems.  Their advantages are well-known;  the flexibility offered by the discontinuous basis functions eases the global finite element assembly and the hanging nodes in mesh generation helps to handle complicated geometry. The continuity restriction for conforming FEM is relaxed, thereby making it an interesting choice for adaptive mesh refinements. On the other hand, conforming finite element methods for plate problems demand C^1 continuity and involve complicated higher-order finite elements. The simplest examples are Argyris finite element with 21 degrees of freedom in a triangle and element with 16 degrees of freedom in a rectangle. 

Nonconforming <cit.>, mixed and hybrid <cit.> finite element methods are also alternative approaches that have been used to relax the C^1-continuity. Discontinuous Galerkin methods are well studied for linear fourth-order elliptic problems, e.g., the hp-version of the nonsymmetric interior penalty DGFEM (NIPG) <cit.>, the hp-version of the symmetric interior penalty DGFEM (SIPG) <cit.> and a combined analysis of NIPG and SIPG in <cit.>. The literature on a posteriori error analysis for biharmonic problems with DGFEM include <cit.> and a quadratic C^0-interior penalty method <cit.>. The medius analysis in <cit.> combines ideas of a priori and a posteriori analysis to establish error estimates for DGFEM under minimal regularity assumptions on the exact solution.  


This paper concerns discontinuous Galerkin finite element methods for the approximation of regular solution to the von Kármán equations, which describe the deflection of very thin elastic plates. Those plates are  modeled by a semi-linear system of fourth-order partial differential equations and can be described as follows. For a given load function f∈, seek u,v ∈ H^2_0(Ω) such that

    Δ^2 u =[u,v]+f  and Δ^2 v =-[u,u]    in Ω

with the biharmonic operator Δ^2 and the von Kármán bracket [∙,∙],
Δ^2φ:=φ_xxxx+2φ_xxyy+φ_yyyy, and [η,χ]:=η_xxχ_yy+η_yyχ_xx-2η_xyχ_xy=(D^2η):D^2χ for the co-factor matrix (D^2η) of D^2η. The colon : denotes the scalar product of two 2× 2 matrices.

In <cit.>, conforming finite element approximations for the von Kármán equations are analyzed and an error estimate in energy norm is derived for approximations of regular solutions. Mixed and hybrid methods reduce the system of fourth-order equations into a system of second-order equations <cit.>. Conforming finite element methods for the canonical von Kármán equations have been proposed and error estimates in energy, H^1 and L^2 norms are established in <cit.> under realistic regularity assumption on the exact solution. Nonconforming FEMs have also been analyzed for this problem <cit.>. An a priori error analysis for a C^0 interior penalty method of this problem is studied in <cit.>. Recently, an abstract framework for nonconforming discretization of a class of semilinear elliptic problems which include is analyzed in <cit.>. 

In this paper, discontinuous Galerkin finite element methods are applied to approximate the regular solutions of the von Kármán equations. To highlight the contribution, under  minimal regularity assumption of the exact solution, optimal order a priori error estimates are obtained and a reliable and efficient a posteriori error estimator is designed. Moreover, a priori and a posteriori error estimates for a C^0 interior penalty method for the von Kármán equations are recovered as a special case. The comprehensive a priori analysis in <cit.> controls the error in the stronger norm ·_h≡∙_ and therefore requires a more involved mathematics and a trilinear form b_IP without symmetry in the first two variables, cf. Remark <ref>.  




The remaining parts of the paper are organized as follows. Section 2 describes some preliminary results and introduces  discontinuous Galerkin finite element methods for von Kármán equations. Section 3 discusses some auxiliary results required for a priori and a posteriori error analysis. In Section 4, a discrete inf-sup condition is established for a linearized problem for the proof of the existence, local uniqueness and error estimates of the discrete solution of the non-linear problem. In Section 5, a reliable and efficient a posteriori error estimator is derived. Section 6 derives a priori and a posteriori error estimates for a C^0 interior penalty method.  Section 7 confirms the theoretical results in various numerical experiments and establishes an adaptive mesh-refining algorithm. 

Throughout the paper, standard notation on Lebesgue and Sobolev spaces and their norms are employed.

The standard semi-norm and norm on H^s(Ω) (resp. W^s,p (Ω)) for s>0 are denoted by |∙|_s and ∙_s (resp. |∙|_s,p and ∙_s,p ).  Bold letters refer to vector valued functions and spaces, e.g. = X× X. The positive constants C appearing in the inequalities denote generic constants which do not depend on the mesh-size. The notation A≲ B means that there exists a generic constant C independent of the mesh parameters and independent of the stabilization parameters σ_1 and σ_2≥ 1 such that A ≤ CB; A≈ B abbreviates A≲ B≲ A.



§ PRELIMINARIES

This section introduces weak and discontinuous Galerkin (dG) formulations for the von Kármán equations.



 §.§ Weak formulation

The weak formulation of von Kármán equations (<ref>) reads: Given f∈, seek  u,v∈  X:= such that

	
    a(u,φ_1)+ b(u,v,φ_1)+b(v,u,φ_1)=l(φ_1)   φ_1∈ X
        a(v,φ_2)-b(u,u,φ_2)   =0            φ_2 ∈ X.


Here and throughout the paper, for all η,χ,φ∈ X, 





    a(η,χ):= D^2 η:D^2χ,   b(η,χ,φ):=- [η,χ]φ,  and  l(φ):=∫_Ωfφ.

Given F=(f,0)∈ L^2(Ω)× L^2(Ω), the combined vector form seeks Ψ=(u,v)∈:=X× X≡× such that

    N(Ψ;Φ):=A(Ψ,Φ)+B(Ψ,Ψ,Φ)-L(Φ)=0Φ∈,

where, for all Ξ=(ξ_1,ξ_2),Θ=(θ_1,θ_2), and Φ=(φ_1,φ_2)∈,

    A(Θ,Φ):=a(θ_1,φ_1)+a(θ_2,φ_2),
       B(Ξ,Θ,Φ):=b(ξ_1,θ_2,φ_1)+b(ξ_2,θ_1,φ_1)-b(ξ_1,θ_1,φ_2) and  L(Φ):=l(φ_1).

Let ∙_2 denote the product norm on   defined by Φ_2:=(|φ_1|_2,Ω^2+|φ_2|_2,Ω^2)^1/2 for all Φ=(φ_1,φ_2)∈. It is easy to verify that the following boundedness and ellipticity properties hold

    A(Θ,Φ)≤Θ_2  Φ_2, A(Θ,Θ) ≥Θ_2^2,
          B(Ξ, Θ, Φ) ≤  C Ξ_2  Θ_2  Φ_2.

Since b(∙,∙,∙) is symmetric in first two variables, the trilinear form B(∙,∙,∙) is symmetric in first two variables.

For results regarding the existence of solution to (<ref>), regularity and bifurcation phenomena, we refer  to <cit.>. It is well known <cit.> that on a polygonal domain Ω, for given f∈ H^-1(Ω), the solutions u,v belong to ∩ H^2+α(Ω), for the index of elliptic regularity α∈ (,1] determined by the interior angles of Ω. Note that when Ω is convex; α=1; that is, the solution belongs to ∩ H^3(Ω). Unless specified otherwise, the parameter α is supposed to satisfy 1/2<α≤ 1. 

Throughout the paper, we consider the approximation of a regular solution <cit.> Ψ to the non-linear operator N(Ψ; Φ)=0 for all Φ∈  of (<ref>) in the sense that the bounded derivative DN(Ψ) of the operator N at the solution Ψ is an isomorphism in the Banach space; this is equivalent to an inf-sup condition 

    0<β:=inf_Θ∈
    Θ_2=1sup_Φ∈
    Φ_2=1(A(Θ,Φ)+2B(Ψ,Θ,Φ)).




 §.§ Triangulations

Let   be a shape-regular <cit.> triangulation of Ω into closed triangles. 

The set of all internal vertices (resp. boundary vertices) and  interior edges (resp.  boundary edges)  of the triangulation  are denoted by (Ω) (resp.  (∂Ω)) and (Ω) (resp. (∂Ω)).
Define a piecewise constant mesh function h_(x)=h_K= diam (K) for all x ∈ K, K∈, and set h:=max_K∈h_K. Also define a piecewise constant edge-function on :=(Ω)∪(∂Ω) by h_|_E=h_E= diam(E) for any E∈. Set of all edges of K is denoted by (K). Note that for a shape-regular family, there exists a positive constant C independent of h such that any K∈ and any E∈∂ K satisfy 

    Ch_K≤ h_E≤ h_K.

Let P_r( K) denote the set of all polynomials of degree less than or equal to r and P_r():={φ∈ L^2(Ω): ∀ K∈,φ|_K∈ P_r(K)} and write P_r():=P_r()× P_r() for pairs of piecewise polynomials.




For a nonnegative integer s, define the broken Sobolev space for the subdivision   as

    H^s()={φ∈: φ|_K∈ H^s(K) K∈}

with the broken Sobolev semi-norm  |∙|_H^s() and norm ∙_H^s() defined by

    |φ|_H^s()=(∑_K∈ |φ|_H^s(K)^2)^1/2 and φ_H^s()=(∑_K∈φ_H^s(K)^2)^1/2.

Define the jump [φ]_E=φ|_K_+-φ|_K_- and the average ⟨φ⟩_E=(φ|_K_++φ|_K_-) across the interior edge E of φ∈ H^1() of the adjacent triangles  K_+ and K_-. Extend the definition of the jump and the average to an edge lying in boundary by [φ]_E=φ|_E and ⟨φ⟩_E=φ|_E for E∈(∂Ω) owing to the homogeneous boundary conditions.  For any vector function, jump and average are understood componentwise.


Set Γ≡⋃_E∈E. 








 §.§ Discrete norms and bilinear forms

For 1/2<α≤ 1, abbreviate Y_h:=(X∩ H^2+α(Ω))+ and :=Y_h× Y_h. For all η,χ∈ Y_h, φ∈ X+, introduce the bilinear, trilinear and linear forms by

    a_(η,χ):=∑_K∈∫_K D^2η:D^2χ-(J(η,χ)+ J(χ,η)) +J_σ_1,σ_2(η,χ),
        b_(η,χ,φ):=-∑_K∈∫_K [η,χ]φ,     l_(φ):= ∑_K∈∫_K  fφ,
       J(η,χ)=∑_E∈∫_E [∇χ]_E·⟨ D^2η ν_E⟩_E ,

with σ_1>0 and σ_2>0 to be suitably chosen in the jump terms across any edge E∈ with unit normal vector ν_E and

    J_σ_1,σ_2(η,χ):=∑_E∈σ_1/h_E^3∫_E[η]_E[χ]_E+∑_E∈σ_2/h_E∫_E[∇η·ν_E]_E[∇χ·ν_E]_E.






The discontinuous Galerkin (dG) finite element formulation of (<ref>) seeks (u_,v_)∈:=× such that, for all (φ_1, φ_2)∈,

    a_(u_,φ_1)+b_(u_,v_,φ_1)+b_(v_,u_,φ_1)=l_(φ_1), 
        a_(v_,φ_2)-b_(u_,u_,φ_2)=0.

The combined vector form seeks Ψ_≡ (u_,v_)∈ such that, for all Φ_∈,

    N_h(Ψ_;Φ_):=A_(Ψ_,Φ_)+B_(Ψ_,Ψ_,Φ_)-L_(Φ_)=0,

where, Ξ_=(ξ_1, ξ_2),Θ_=(θ_1,θ_2),Φ_=(φ_1,φ_2)∈,

    A_(Θ_,Φ_):=a_(θ_1,φ_1)+a_(θ_2,φ_2),
       B_(Ξ_,Θ_,Φ_):=b_(ξ_1,θ_2,φ_1)+b_(ξ_2,θ_1,φ_1)-b_(ξ_1,θ_1,φ_2), 
       L_(Φ_):=l_(φ_1).

Note that b_(∙,∙,∙) is symmetric in the first and second variables, and so is B_(∙,∙,∙).

For φ∈ H^2() and Φ=(φ_1,φ_2)∈^2()≡ H^2()× H^2(), define the mesh dependent  norms ∙_ and ∙_ by 

    φ_^2:=|φ|_ H^2()^2+∑_E∈σ_1/h_E^3[φ]_E_L^2(E)^2+∑_E∈σ_2/h_E[∇φ·ν_E]_E_L^2(E)^2,
       Φ_^2:=φ_1_^2+φ_2_^2.
 

For ξ∈ Y_h≡ (X∩ H^2+α(Ω))+ and Ξ=(ξ_1,ξ_2)∈≡ Y_h× Y_h, define the auxiliary norms ∙_h and ∙_h  by

    ξ_h^2:=ξ_^2+∑_E∈∑_j,k=1^2h_E^1/2⟨∂^2ξ/∂ x_j ∂ x_k⟩_E  _L^2(E)^2 and Ξ_h^2:=ξ_1_h^2+ξ_2_h^2.




§ AUXILIARY RESULTS 

This section discusses some auxiliary results and establishes the boundedness and ellipticity results required for the analysis.


 §.§ Some known operator bounds

This subsection recalls a few standard results. Throughout this subsection, the generic multiplicative constant C≈ 1 hidden in the brief notation ≲ depends on the shape regularity of the triangulation  and arising parameters like the polynomial degree r∈_0 or the Lebesgue index p and the Sobolev indices ℓ,s>1/2 and 1/2<α≤ 1; C is independent of the mesh-size.

 <cit.> There exists a positive constant C that depends on ℓ∈ℕ, p∈ and the shape regularity of the triangulation such that, for 1≤ℓ,  2≤ p≤∞, for any ξ∈ P_r(K) satisfies
	
    ξ_L^p(K)   ≲ h_K^(2-p)/pξ_L^2(K) and |ξ|_H^ℓ(K)≲ h_K^-1|ξ|_H^ℓ-1(K),
	
	for any K ∈ with E⊂(K), where
	
    ξ_L^p(E)≲ h_E^1/p-1/2ξ_L^2(E).


 The following trace inequalities hold for K∈ and s>1/2.
	
		
  * <cit.> ξ_L^2(∂ K)≲ h_K^-1/2ξ_L^2(K)ξ∈ P_r(K);
		
  * <cit.>
       ξ_L^2(∂ K)≲h_K^s-1/2ξ_H^s(K)+h_K^-1/2ξ_L^2(K) for all ξ∈ H^s(K).	
		


 <cit.> There exists a linear operator
Π_h: H^s()→ P_r(), such that, for 0≤ q≤ s, m=min(r+1,s) and 1/2<α≤ 1,

    φ-Π_hφ_H^q(K)   ≲ h_K^m-qφ_H^s(K) K∈ andφ∈ H^s(),
    φ-Π_hφ_   ≤φ-Π_hφ_h≲ h^αφ_2+αφ∈ H^2+α(Ω).



The proof of (<ref>) follows from Lemma <ref>.b, (<ref>), and an interpolation of Sobolev spaces <cit.>.















For an easy notation of vectors, we denote the componentwise interpolation of ζ∈H^s():=H^s()× H^s() by Π_hζ.





	<cit.> For K∈, a macro-element of degree 4 is a nodal finite element (K,P̃_4,Ñ), consisting of sub-triangles K_j, j=1,2,3 (see Figure <ref>). The local element space P̃_4 is defined by 
	
    P̃_4:={φ∈ C^1(K): φ|_K_j∈ P_4(K_j), j=1,2,3}.
	
	The degrees of freedom Ñ are defined as (a) the value and the first (partial) derivatives at the vertices of K; (b) the value at the midpoint of each edge of K; (c) the normal derivative at two distinct points in the interior of each edge of K; (d) the value and the first (partial) derivatives at the common vertex of  K_1, K_2 and K_3. The corresponding finite element space consisting of the above macro-elements will be denoted by S_4()⊂.







The enrichment operator of <cit.> is outlined in the sequel for a convenient reading. For each nodal point p of the C^1-conforming finite element space S_4(), define (p) to be the set of K∈ which shares the nodal point p and let |(p)| denotes its cardinality. Define the operator E_h:→ S_4() for any nodal variable N_p at p by

    N_p(E_h(φ_)):=1/|(p)|∑_K∈(p) N_p(φ_|_K)  if  p∈(Ω),
    
    0  if  p∈(∂Ω).




	
<cit.> The enrichment operator E_h: → S_4() satisfies, for m=0,1,2,

    ∑_K∈|φ_-E_hφ_|_H^m(K)^2    ≲h_^1/2-m[φ_]__L^2(Γ)^2+h_^3/2-m[∇φ_]__L^2(Γ)^2 ≲ h^2-mφ__.

		Moreover, for some positive constant Λ≈ 1,
		
    φ_-E_hφ__≤Λinf_φ∈ Xφ_-φ_.

	
	
		See <cit.> for a proof of (<ref>). For the proof of (<ref>), choose m≤ 2 in (<ref>) and obtain (with h_≲ h≲ 1) that
		
    φ_-E_hφ__ H^2()^2≲h_^-3/2[φ_]__L^2(Γ)^2+h_^-1/2[∇φ_]__L^2(Γ)^2.

		Since [φ_-E_hφ_]_E=[φ_]_E and [∇(φ_-E_hφ_)]_E=[∇φ_]_E,  those edge terms in both sides of the above inequality lead (in the definition of ·_) to
		
    φ_-E_hφ__^2≲h_^-3/2[φ_]__L^2(Γ)^2+h_^-1/2[∇φ_]__L^2(Γ)^2.

		Furthermore, any φ∈ X satisfies (with (<ref>) for m=2 in the end) that
		
    φ_-E_hφ__^2≲h_^-3/2[φ_-φ]__L^2(Γ)^2+h_^-1/2[∇ (φ_-φ)]__L^2(Γ)^2≲φ_-φ_.

		This completes the proof of (<ref>) for some h-independent positive constant Λ.
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	

		It holds
		
    h_∇φ_L^∞(Ω)   ≲φ_L^∞(Ω)φ∈+S_4(),
    φ_W^1,4()+φ_L^∞(Ω)   ≲φ_φ∈+X.

	

This follows with the arguments of <cit.> on the enrichment and interpolation operator. Further details are omitted for brevity.

















































 §.§ Continuity and ellipticity




This subsection is devoted to the boundedness and ellipticity results for the bilinear form a_(∙, ∙) and boundedness results for b_(∙, ∙,∙).











 Any θ_,φ_∈+S_4() satisfies
	
    a_(θ_,φ_)≲θ__φ__.

	



	Given any θ_,φ_∈, recall the definition of a_(∙,∙)
	
    a_(θ_,φ_)    = D^2θ_: D^2φ_-(J(θ_,φ_)+ J(φ_,θ_)) 
            +J_σ_1,σ_2(θ_,φ_).

	The definition of J(∙, ∙), the Cauchy-Schwarz inequality, and Lemma <ref> imply
	
    J(θ_,φ_)= [∇φ_]_E·⟨ D^2θ_ν_E⟩_E 
       ≤σ_2^-1/2(∑_E∈σ_2/h_E[∇φ_]_E_L^2(E)^2)^1/2(∑_E∈h_E^1/2⟨ D^2θ_⟩_E  _L^2(E)^2)^1/2
       ≲σ_2^-1/2φ__|θ_|_H^2()≤σ_2^-1/2θ__φ__.

	The same arguments show, J(φ_,θ_)≲σ_2^-1/2θ__φ__. The definitions of J_σ_1,σ_2(∙,∙) and ∙_, and the Cauchy-Schwarz inequality lead to 
	
    J_σ_1,σ_2(θ_,φ_)=∑_E∈σ_1/h_E^3∫_E[θ_]_E[φ_]_E+∑_E∈σ_2/h_E∫_E[∇θ_·ν_E]_E[∇φ_·ν_E]_E
       ≤(∑_E∈σ_1/h_E^3[θ_]_E_L^2(E)^2)^1/2(∑_E∈σ_1/h_E^3[θ_]_E_L^2(E)^2)^1/2
         +(∑_E∈σ_2/h_E[∇θ_·ν_E]_E_L^2(E)^2)^1/2(∑_E∈σ_2/h_E[∇φ_·ν_E]_E_L^2(E)^2)^1/2≲θ_φ_.

	The combination of all displayed formulas and σ_2≥ 1 conclude the proof.





	The definitions of a_(∙,∙), the auxiliary norm ∙_h, and the estimate (<ref>) imply (since σ_2≥ 1)
	
    a_(θ,φ)≲θ_hφ_hθ,φ∈ Y_h≡ (X∩ H^2+α(Ω))+.




	
	The trace inequality Lemma <ref>.a implies that ∙_≈∙_h are equivalent norms on +S_4() with equivalence constants, which do neither depend on the mesh-size nor on σ
	_1,σ_2>0.

 
	For any σ_1>0 and for a sufficiently large parameter σ_2, there exists some h-independent positive constant β_0 (which depends on σ_2) such that
	
    β_0 θ__^2 ≤ a_(θ_,θ_) θ_∈.



	For  θ_∈, the definition of a_(∙,∙) leads to
	
    a_(θ_,θ_)=θ__^2-2J(θ_,θ_).

Recall (<ref>) in the form 
J(θ_,θ_)≤ C_0 σ_2^-1/2θ__^2
with some constant C_0≈ 1. For any 0<β_0<1 and any choice of σ_2≥ 4C_0(1-β_0)^-2, the combination of the previous estimates concludes the proof.































Recall that h denotes the maximal mesh-size of the underlying triangulation .

	Any ξ∈ H^2+α(Ω) with 1/2<α≤ 1 and φ_∈ satisfy
	
    a_(ξ,φ_-E_hφ_)≲ h^αξ_2+αφ__.

	Consequently, for ξ∈^2+α(Ω) and Φ_∈,
	
    A_(ξ,Φ_-E_hΦ_)≲ h^αξ_2+αΦ__.




	The method of real interpolation <cit.> of Sobolev spaces defines H^2+α(Ω):=[ H^2(Ω),H^3(Ω)]_α,2 from H^2(Ω) and H^3(Ω).
	Given any ξ∈ H^2(Ω) and φ_∈, the boundedness of a_(∙,∙) and Lemma <ref> imply
	
    a_(Π_hξ,φ_-E_hφ_)≲Π_hξ_φ_-E_hφ__≲ξ_2φ__.

	
	
	
	
	Given any ξ∈ H^3(Ω) and φ_∈, the definition of a_(∙,∙), an integration by parts, and Lemma <ref> result in
	
    a_(ξ,φ_-E_hφ_)
    	   =∑_K∈∫_K D^2ξ:D^2(φ_-E_hφ_)- J(ξ,φ_-E_hφ_)
       =∑_K∈∫_K(D^2ξ)·∇(φ_-E_hφ_)≤ Chξ_3φ__.

	This,  Remark <ref>, (<ref>) with α=1, and Lemma <ref> lead to
	
    a_(Π_hξ,φ_-E_hφ_)   =a_(Π_hξ-ξ,φ_-E_hφ_)+a_(ξ,φ_-E_hφ_) 
       ≲ hξ_3φ__.

	The estimates (<ref>)-(<ref>) and the interpolation between Sobolev spaces show
	
    a_(Π_hξ,φ_-E_hφ_)≲ h^αξ_2+αφ__.

	For any α>1/2, the combination with Remark <ref>, (<ref>), and Lemma <ref> results in
	
    a_(ξ,φ_-E_hφ_)    =a_(ξ-Π_hξ,φ_-E_hφ_)+a_(Π_hξ,φ_-E_hφ_) 
       ≲ h^αξ_2+αφ__.

	This concludes the proof of (<ref>). The estimate (<ref>) follows from (<ref>) and the definition of A_(∙,∙).











 
	
		
  * Any  η,χ, φ∈ X+ satisfy 
    b_(η,χ,φ)≲η_χ_φ_.

		
  * Given any α>1/2, any η∈ X∩ H^2+α(Ω) and χ∈ X+ satisfy
		
    b_(η,χ,φ)≲η_2+αχ_φ_1φ∈ H^1_0(Ω),
    η_2+αχ_φ_L^4(Ω)φ∈ X+.

		



	(a).
	For η,χ,φ∈ X+, the definition of b_(∙,∙,∙) and Lemma <ref> lead to
	
    |2b_(η,χ,φ)|   =| [η,χ]φ|≲ |η|_ H^2()|χ|_ H^2()η_L^∞()
       ≲η_χ_φ_.

	
	
	(b). For η∈ X∩ H^2+α(Ω), χ∈ X+, and φ∈ H^1_0(Ω)∪ (X+), the generalized inequality and the continuous imbedding H^2+α(Ω)↪ W^2,4(Ω) yields
	
    |2b_(η,χ,φ)|   =| [η,χ]φ|≲η_W^2,4(Ω)χ_H^2()φ_L^4(Ω)
       ≲η_2+αχ_φ_L^4(Ω).

	This implies second part of (b). For φ∈ H^1_0(Ω)↪ L^4(Ω) this proves the first.




























































































§ A PRIORI ERROR CONTROL


This section first establishes the discrete inf-sup condition for the linearized problem, then the existence of a discrete solution to the nonlinear problem (<ref>), and finally the convergence of a Newton method. 



 §.§ Discrete inf-sup condition

This subsection is devoted to discrete inf-sup condition. Throughout the paper, the statement “there exists σ_2 such that for all σ_2≥σ_2 as in Lemma <ref> on ellipticity, there exists h(σ_2) such that for all h≤ h(σ_2) ...” is abbreviated by the phrase “for sufficiently large σ_2 and sufficiently small h ...”.

 Let Ψ∈^2+α(Ω)∩^2_0(Ω) be a regular solution to (<ref>). For sufficiently large σ_2 and sufficiently small h, there exists β such that the following discrete inf-sup condition holds
	
    0<β≤inf_Θ_∈
    Θ__=1sup_Φ_∈
    Φ__=1(A_(Θ_,Φ_)+2B_(Ψ,Θ_,Φ_)).



	Given any Θ_∈ with Θ__=1, let ξ∈ and η∈ solve the biharmonic problems
	
    A(ξ,Φ)   =2B_(Ψ,Θ_,Φ)Φ∈,
    
    	A(η,Φ)   =2B(Ψ,E_hΘ_,Φ)Φ∈.

	Lemma <ref>.b implies that B_(Ψ,Θ̃,∙) and B (Ψ,Θ̃,∙) belong to ^-1(Ω) for Θ̃∈+. The reduced elliptic regularity for the biharmonic problem <cit.> yields ξ,η∈^2+α(Ω)∩. Since Ψ is a regular solution to (<ref>), there exists β from (<ref>) and Φ∈ with Φ_2=1 such that
	
    β E_hΘ__2≤ A(E_hΘ_,Φ)+2B(Ψ,E_hΘ_,Φ).

	The solution property in (<ref>), the boundedness of A(∙,∙), and the triangle inequality in the above result imply
	
    β E_hΘ__2    ≤ A(E_hΘ_+η,Φ)≤ E_hΘ_+η_2 
       ≤  E_hΘ_-Θ__+Θ_+ξ_+η-ξ_2.

	The definition of ξ,η in (<ref>)-(<ref>) and Lemma <ref>.a lead to
	
    η-ξ_2≤ 2C_bΨ_2 E_hΘ_-Θ__

	for some positive constant C_b≈ 1.
	The combination of the previous two displayed inequalities reads
	
    β E_hΘ__2≤Θ_+ξ_+(1+2C_bΨ_2) E_hΘ_-Θ__.

	This and (<ref>) result in 
	
    E_hΘ__2≤1/β(1+Λ(1+2C_bΨ_2))Θ_+ξ_.

	The triangle inequality, (<ref>), and (<ref>) lead to
	
    1   =Θ__≤ E_hΘ_-Θ__+ E_hΘ__2
       ≤(Λ+1/β(1+Λ(1+2C_bΨ_2)))Θ_+ξ_.

	In other words, β_1:=β/(1+Λ(1+β+2C_bΨ_2)) satisfies
	
    β_1≤Θ_+ξ_≤Θ_+Π_hξ_+ξ-Π_hξ_.

	For any given Θ_+Π_hξ∈, the ellipticity of A_(∙,∙) from Lemma <ref> implies the existence of some Φ_∈ with Φ__=1 and
	
    β_0Θ_+Π_hξ_≤ A_(Θ_+Π_hξ,Φ_) 
    
    	        =A_(Θ_,Φ_)+A_(Π_hξ-ξ,Φ_)+A_(ξ,Φ_-E_hΦ_)+A(ξ,E_hΦ_).

	The choice of Φ=E_hΦ_ in (<ref>) plus straightforward calculations result in
	
    β_0Θ_+Π_hξ_
    	
    	   ≤ A_(Θ_,Φ_)+ A_(Π_hξ-ξ,Φ_) 
    	+A_(ξ,Φ_-E_hΦ_) 
          + 2B_(Ψ,Θ_,Φ_)
    	+2B_(Ψ,Θ_,E_hΦ_-Φ_).

	
	Remark <ref>, Lemma <ref> and <ref> for the second and third term, plus Lemma <ref>.b and <ref> for the last term in (<ref>) lead to 	
	
    β_0Θ_+Π_hξ_≤ A_(Θ_,Φ_)+2B_(Ψ,Θ_,Φ_)+Ch^α.

	The combination of (<ref>) and (<ref>) with Lemma <ref> shows
	
    β_0β_1-C_*h^α≤ A_(Θ_,Φ_)+2B_(Ψ,Θ_,Φ_)

	for some h-independent positive constant C_*. Hence,  for all  h≤ h_0:=(β_0β_1/2C_*)^1/α, the discrete inf-sup condition (<ref>) follows.


The next lemma establishes that the perturbed bilinear form  

    _(Θ_,Φ_):= A_(Θ_,Φ_)+2B_(Π_hΨ,Θ_,Φ_)

satisfies a discrete inf-sup condition.

 
    Let Π_hΨ be the interpolation of Ψ from Lemma <ref>.
	Then, for sufficiently large σ_2 and sufficiently small h, the perturbed bilinear form (<ref>) satisfies the discrete inf-sup condition
	
    β/2≤inf_Θ_∈
    Θ__=1sup_Φ_∈
    Φ__=1_(Θ_,Φ_).

	
	
	

Proof.
	
	
	
	
	
	Lemma <ref> leads to the existence of h_1>0, such that  Ψ-Π_hΨ_≤β/4C_b holds for h≤ h_1. Given any Θ_∈,
	Theorem <ref> and  Lemma <ref>.a lead to 
	
    β/2Θ__≤βΘ__-2C_bΨ-Π_hΨ_Θ_
       ≤sup_Φ_∈
    Φ__=1(A_(Θ_,Φ_)+2B_(Ψ,Θ_,Φ_)) 
                   -2sup_Φ_∈
    Φ__=1 B_(Ψ-Π_hΨ,Θ_,Φ_)
       ≤sup_Φ_∈
    Φ__=1(A_(Θ_,Φ_)+2B_(Π_hΨ,Θ_,Φ_))=sup_Φ_∈
    Φ__=1_(Θ_,Φ_).














 §.§ Existence, uniqueness and error estimate

The discrete inf-sup condition is employed to define a nonlinear map μ:→ which enables to analyze the existence and uniqueness of solution of (<ref>). 
For any Θ_∈, define μ(Θ_)∈ as the solution to the discrete fourth-order problem

    _(μ(Θ_),Φ_)=L_(Φ_)+2B_(Π_hΨ,Θ_,Φ_)-B_(Θ_,Θ_,Φ_)
 
for all Φ_∈.
Lemma <ref> guarantees that the mapping μ is well-defined and continuous. Also, any fixed point of μ is a solution to (<ref>) and vice-versa. In order to show that the mapping μ has a fixed point, define the ball 
    _R(Π_hΨ):={Φ_∈: Φ_-Π_hΨ_≤ R }.



	For sufficiently large σ_2 and sufficiently small h, there exists a positive constant R(h) such that μ maps the ball _R(h)(Π_hΨ) to itself; μ(Θ_ )-Π_hΨ_≤ R(h) holds for any Θ_∈_R(h)(Π_hΨ).



	The discrete inf-sup condition of _(∙,∙) in Lemma <ref> implies the existence of Φ_∈ with  Φ__=1 and 
	
    β/2μ(Θ_)-Π_hΨ_≤_(μ(Θ_)-Π_hΨ, Φ_).

	Let E_hΦ_ be the enrichment of Φ_ from Lemma <ref>.   The definition of _(∙,∙), the symmetry of B_(∙,∙,∙) in first and second variables, (<ref>), and (<ref>) lead to 
	
    _(μ(Θ_)-Π_hΨ, Φ_)
    	=_(μ(Θ_), Φ_)-_(Π_hΨ, Φ_)
         =L_(Φ_)+2B_(Π_hΨ,Θ_,Φ_)-B_(Θ_,Θ_,Φ_)-A_(Π_hΨ,Φ_)
             -2B_(Π_hΨ,Π_hΨ,Φ_)
         =L_(Φ_-E_hΦ_) +(A(Ψ, E_hΦ_)-A_(Π_hΨ,Φ_)) 
             +(B(Ψ,Ψ, E_hΦ_)-B_(Π_hΨ,Π_hΨ,Φ_))
    	
    	+B_(Π_hΨ-Θ_,Θ_-Π_hΨ,Φ_)
          =:T_1+T_2+T_3+T_4.

	The term T_1 can be estimated using the continuity of L_ and Lemma <ref>.
	The continuity of A_(∙, ∙), Lemma <ref> and <ref> with Φ__=1  lead to
	
    T_2:   =A(Ψ, E_hΦ_)-A_(Π_hΨ,Φ_)
        = A_(Ψ,E_hΦ_-Φ_)
    	+A_(Ψ-Π_hΨ,Φ_) ≲ h^αΨ_2+α.

	Lemma <ref>, <ref>, and <ref> result in 
	
    T_3   := B(Ψ,Ψ,E_hΦ_)-B_(Π_hΨ,Π_hΨ,Φ_)
       = B_(Ψ,Ψ-Π_hΨ,E_hΦ_)+B_(Ψ-Π_hΨ,Π_hΨ,Φ_)
    	
         +B_(Ψ,Π_hΨ,E_hΦ_-Φ_) ≲ h^αΨ_2+αΨ_2.
 
	
	
	
	
	
	
	
	
	Lemma <ref> implies
	
    T_4:=B_(Π_hΨ-Θ_,Θ_-Π_hΨ,Φ_)≲Θ_-Π_hΨ_^2.

	A substitution of the estimates for T_1, T_2, T_3, and T_4 in (<ref>) and Ψ_2+α≈ 1≈Ψ_2 lead to C_1≈ 1 with
	
    μ(Θ_)-Π_hΨ_≤ C_1( h^α+Θ_-Π_hΨ_^2).

	Then h≤ h_2:=(2C_1)^-2/α and Θ_-Π_hΨ_≤ R(h):=2C_1 h^α lead to
	
    μ(Θ_)-Π_hΨ_≤ C_1 h^α(1+4C_1^2 h^α)≤ R(h).

	This concludes the proof.









For sufficiently large σ_2 and sufficiently small h, there exists a unique solution Ψ_ to the discrete problem (<ref>) in _R(h)(Π_hΨ).




	First we prove the contraction result of the nonlinear map μ in the ball _R(h)(Π_hΨ) of Theorem <ref>. Given any Θ_,Θ̃_∈_R(h)(Π_hΨ) and for all Φ_∈, the solutions μ(Θ_) and  μ(Θ̃_) satisfy
	
    _(μ(Θ_),Φ_)   =L_(Φ_)+2B_(Π_hΨ,Θ_,Φ_)-B_(Θ_,Θ_,Φ_), 
    _(μ(Θ̃_),Φ_)   =L_(Φ_)+2B_(Π_hΨ,Θ̃_,Φ_)-B_(Θ̃_,Θ̃_,Φ_).

	The discrete inf-sup of _(∙,∙) from Lemma <ref> guarantees the existence of Φ_ with Φ__=1 below.  With (<ref>)-(<ref>) and Lemma <ref>, it follows that
	
    β/2μ(Θ_)-μ(Θ̃_)_≤_(μ(Θ_)-μ(Θ̃_),Φ_)
       =2B_(Π_hΨ,Θ_-Θ̃_, Φ_)+B_(Θ̃_,Θ̃_, Φ_)-B_(Θ_,Θ_, Φ_)
       =B_(Θ̃_-Θ_,Θ_-Π_hΨ,Φ_)+B_(Θ̃_-Π_hΨ,Θ̃_-Θ_,Φ_)
       ≲Θ̃_-Θ__(Θ_-Π_hΨ_+
    	Θ̃_-Π_hΨ_).

	Since Θ_,Θ̃_∈_R(h)(Π_hΨ), for a choice of R(h) as in the proof of  Theorem <ref>, for sufficiently large σ_2 and h≤min{h_0,h_1,h_2}, 
	
    μ(Θ_)-μ(Θ̃_)_≲ h^αΘ̃_-Θ__.

	Hence, there exists positive constant h_3, such that for h≤ h_3 the contraction result holds.
	
	





	For h≤h:=min{h_0,h_1,h_2,h_3}, Lemma <ref> and Theorem <ref> lead to the fact that μ is a contraction map that maps the ball _R(h)(Π_hΨ) into itself. An application of the Banach fixed point theorem yields that the mapping μ has a unique fixed point in the ball _R(h)(Π_hΨ), say Ψ_, which solves (<ref>) with Ψ_-Π_hΨ_≤ R(h).




	Let  Ψ be a regular solution to (<ref>)  and let Ψ_ be the solution to (<ref>). For sufficiently large σ_2 and sufficiently small h, it holds
	
    Ψ-Ψ__≤ C h^α.




	A triangle inequality yields
	
    Ψ-Ψ__≤Ψ-Π_hΨ_+Π_hΨ-Ψ__.

	For h≤h and sufficiently large σ_2, Theorem <ref> leads to
	
    Π_hΨ-Ψ__≤ Ch^α.

	This, Lemma <ref>, (<ref>),  and (<ref>) conclude the proof.






 §.§ Convergence of the Newton method


The discrete solution Ψ_ of (<ref>) is characterized as the fixed point of (<ref>) and so depends on the unknown Π_hΨ. The approximate solution to (<ref>) is computed with the Newton method, where the iterates Ψ_^j solve 

    A_(Ψ_^j,Φ_)+2B_(Ψ_^j-1,Ψ_^j,Φ_)=B_(Ψ_^j-1,Ψ_^j-1,Φ_)+L_(Φ_)
 
for all Φ_∈. The Newton method has locally quadratic convergence.

	Let Ψ be a regular solution to (<ref>) and let  Ψ_ solve (<ref>). There exists a positive constant R independent of h, such that for any initial guess Ψ_^0 with
	Ψ_- Ψ_^0_≤  R, it follows Ψ_- Ψ_^j_≤ R j=0,1,2,… and the iterates of the Newton method  in (<ref>) are well defined and converges quadratically to Ψ_.


	Following the proof of Lemma <ref>, there exists a positive constant ϵ (sufficiently small) independent of h such that for each Z_∈ with Z_-Π_hΨ_≤ϵ, the bilinear form
	
    A_(∙,∙)+2B_(Z_,∙,∙)

	satisfies discrete inf-sup condition in ×.
	 For sufficiently large σ_2 and sufficiently small h, the equation (<ref>) implies Π_hΨ-Ψ__≤ Ch^α. Thus h can be chosen sufficiently small so that Π_hΨ-Ψ__≤ϵ/2. Recall β from (<ref>). Lemma <ref>.a implies that there exists a positive constant C_b̃ independent of h such that  B_(Ξ_,Θ_,Φ_)≤ C_b̃Ξ__Θ__Φ__. Set
	
    R:=min{ϵ/2,β/8C_b̃}.

	Assume that the initial guess Ψ_^0 satisfies Ψ_-Ψ_^0_≤  R. Then,
	
    Π_hΨ-Ψ_^0_≤Π_hΨ-Ψ__+Ψ_-Ψ_^0_≤ϵ.

	This implies Ψ_-Ψ_^j-1_≤  R for j=1 and suppose for mathematical induction that this holds for some j∈. Then Z_:=Ψ_^j-1 in (<ref>) leads to an discrete inf-sup condition of A_(∙,∙)+2B_(Ψ^j-1_,∙,∙) and so to an unique solution Ψ_^j in step j of the Newton scheme. The discrete inf-sup condition (<ref>) implies the existence of Φ_∈ with Φ__=1 and
	
    β/4Ψ_-Ψ_^j_≤ A_(Ψ_-Ψ_^j, Φ_)+2B_(Ψ_^j-1,Ψ_-Ψ_^j, Φ_).

	The application of (<ref>), (<ref>), and Lemma <ref> result in 
	
    A_(Ψ_-Ψ_^j, Φ_)+2B_(Ψ_^j-1,Ψ_-Ψ_^j, Φ_)
       =A_(Ψ_,  Φ_)+2B_(Ψ_^j-1,Ψ_,  Φ_)-B_(Ψ_^j-1,Ψ_^j-1,  Φ_)-L_(  Φ_)
       =-B_(Ψ_,Ψ_,Φ_)+2B_(Ψ_^j-1,Ψ_,  Φ_)-B_(Ψ_^j-1,Ψ_^j-1,  Φ_)
       =B_(Ψ_-Ψ_^j-1,Ψ_^j-1-Ψ_,  Φ_)≤ C_b̃Ψ_-Ψ_^j-1_^2.

	This implies 
	
    Ψ_-Ψ_^j_≤(4C_b̃/β)Ψ_-Ψ_^j-1_^2

	and establishes the quadratic convergence of the Newton method  to Ψ_. The definition of R and (<ref>) guarantee Ψ_-Ψ_^j_≤Ψ_-Ψ_^j-1_<R to allow an induction step j→ j+1 to conclude the proof.










§ A POSTERIORI ERROR CONTROL


This section establishes a reliable and efficient error estimator for the DGFEM. For K ∈ and E ∈(Ω), define the volume and edge estimators η_K and η_E by 

    η_K^2   := h_K^4(f+[u_,v_]_L^2(K)^2+[u_,u_]_L^2(K)^2),
    η_E^2   :=h_E([D^2 u_ ν_E]_E_L^2(E)^2+[ D^2 v_ν_E]_E_L^2(E)^2)
    
         +h_E^-3([u_]_E_L^2(E)^2+[v_]_E_L^2(E)^2)+h_E^-1([∇ u_]_E_L^2(E)^2+[∇ v_]_E_L^2(E)^2).










 Let Ψ=(u,v)∈ be a regular solution to (<ref>) and let Ψ_=(u_,v_)∈ solve (<ref>). For sufficiently large σ_2 and sufficiently small h, there exists an h-independent positive constant C_ rel such that
	
    Ψ-Ψ__^2   ≤ C_ rel^2(∑_K∈η_K^2+∑_E∈(Ω)η_E^2).




	The Fréchet derivative of N(Ψ) at Ψ in the direction Θ∈ reads
	
    DN(Ψ;Θ,Φ):=A(Θ,Φ)+2B(Ψ,Θ,Φ)Φ∈.

	Since Ψ is a regular solution, for any 0<ϵ<β with β from (<ref>), there exists some Φ∈ with Φ_2=1 and
	
    (β-ϵ)Ψ-E_hΨ__2≤ DN(Ψ;Ψ-E_hΨ_,Φ).

	Since N is quadratic, the finite Taylor series is exact and shows
	
    N(E_hΨ_;Φ)   =N(Ψ;Φ)+DN(Ψ;E_hΨ_-Ψ,Φ)
         + D^2N(Ψ;E_hΨ_-Ψ, E_hΨ_-Ψ,Φ).

	Since N(Ψ;Φ)=0 and D^2N(Ψ;Θ,Θ,Φ)=2B(Θ,Θ,Φ) for Θ=Ψ-E_hΨ_, (<ref>)-(<ref>) plus Lemma <ref>.a with boundedness constant C_b lead to
	
    (β-ϵ)Ψ-E_hΨ__2   ≤ -N(E_hΨ_;Φ)+B(Ψ-E_hΨ_, Ψ-E_hΨ_,Φ)
       ≤ |N(E_hΨ_;Φ)|+C_bΨ-E_hΨ__2^2.

	The triangle inequality, (<ref>), and Theorem <ref> imply 
	
    Ψ-E_hΨ__≤Ψ-Ψ__+Ψ_-E_hΨ__
    	≤ C(1+Λ)h^α.

	With ϵ↘ 0, (<ref>)-(<ref>) verify
	
    (β-(1+Λ)CC_bh^α)Ψ-E_hΨ__2≤ |N(E_hΨ_;Φ)|.

	There exists positive constant h_4, such that h≤ h_4 implies β-(1+Λ)CC_bh^α>0. Hence, for h≤ h_4, the above equation and triangle inequality lead to
	
    Ψ-Ψ__   ≤Ψ-E_hΨ__+ E_hΨ_-Ψ__
       ≲ |N(E_hΨ_;Φ)|+ E_hΨ_-Ψ__.

	
	
	
	
	The definitions of N and  N_h, N_h(Ψ_;Π_hΦ)=0, the boundedness of A_(∙, ∙), B_(∙, ∙, ∙), and Theorem <ref> result in
	
    N(E_hΨ_;Φ)=A(E_hΨ_,Φ)+B(E_hΨ_,E_hΨ_,Φ)-L(Φ)
       =N_h(Ψ_;Φ)+A_(E_hΨ_-Ψ_,Φ)+B(E_hΨ_,E_hΨ_,Φ)-B_(Ψ_,Ψ_,Φ)
       =N_h(Ψ_;Φ-Π_hΦ)+A_(E_hΨ_-Ψ_,Φ)
           +B_(E_hΨ_-Ψ_,E_hΨ_,Φ)+B_(Ψ_,E_hΨ_-Ψ_,Φ)
       ≲ |N_h(Ψ_;Φ-Π_hΦ)|+ E_hΨ_-Ψ__.

	The combination of the previous displayed estimates proves 
	
    Ψ-Ψ__≲ |N_h(Ψ_;Φ-Π_hΦ)|+ E_hΨ_-Ψ__.
	
	Abbreviate Φ-Π_hΦ=:χ=(χ_1,χ_2) and recall
	
    N_h(Ψ_;χ)=A_(Ψ_,χ)+B_(Ψ_,Ψ_,χ)-L_(χ).

	The definition of A_(∙,∙) and an integration by parts with the facts that  u_ and v_ are piecewise quadratic polynomials lead to
	
    A_(Ψ_,χ)=(D^2 u_: D^2χ_1+D^2 v_: D^2χ_2)-J(u_,χ_1)- J(χ_1,u_)
               -J(v_,χ_2)- J(χ_2,v_) +∑_E∈σ_1/h_E^3∫_E( [u_]_E[χ_1]_E+[v_]_E[χ_2]_E)
                 +∑_E∈σ_2/h_E∫_E([∇ u_·ν_E]_E[∇χ_1·ν_E]_E+[∇ v_·ν_E]_E[∇χ_2·ν_E]_E)
       =([D^2 u_ ν_E]_E·⟨∇χ_1⟩_E  +[D^2 v_ ν_E ]_E·⟨∇χ_2⟩_E  )
           -(J(χ_1,u_)+J(χ_2,v_))
    	+∑_E∈σ_1/h_E^3∫_E([u_]_E[χ_1]_E+[v_]_E[χ_2]_E)
             +∑_E∈σ_2/h_E∫_E([∇ u_·ν_E]_E[∇χ_1·ν_E]_E+[∇ v_·ν_E]_E[∇χ_2·ν_E]_E).

	The trace inequality of Lemma <ref> and interpolation estimates χ_1_H^m(K)=φ_1-Π_hφ_1_H^m(K)≤ h_K^2-mφ_1_H^2(K) for m=1,2, from Lemma <ref> result in
	
    h_^-1/2⟨∇χ_1⟩__L^2(Γ)^2≲∑_K∈h_K^-1∇χ_1_L^2(∂ K)^2
       ≲∑_K∈h_K^-1(h_K^-1χ_1_H^1(K)^2+h_Kχ_1_H^2(K)^2)≲Φ_2^2=1.

	This and the Cauchy-Schwarz inequality imply 
	
    ([D^2 u_ ν_E]_E·⟨∇χ_1⟩_E  +[D^2 v_ ν_E ]_E·⟨∇χ_2⟩_E  )
    
    	     ≲(∑_E∈(Ω)h_E^1/2[D^2 u_ ν_E]_E_L^2(E)^2)^1/2+(∑_E∈(Ω)h_E^1/2[D^2 v_ ν_E]_E_L^2(E)^2)^1/2.

	The Cauchy-Schwarz inequality, the trace inequality of Lemma <ref>, and the interpolation of Lemma <ref> lead to
	
    J(χ_1,u_)+J(χ_2,v_)≲(h_^-1/2[∇ u_]_E_Γ^2+h_^-1/2[∇ v_]_E_Γ^2)^1/2Φ_2.

	Similar arguments for the penalty terms result in
	
    ∑_E∈σ_1/h_E^3∫_E( [u_]_E[χ_1]_E+[v_]_E[χ_2]_E)≲(h_^-3/2[u_]__L^2(Γ)^2+h_ ^-3/2[v_]__L^2(Γ)^2)^1/2,
       ∑_E∈σ_2/h_E∫_E([∇ u_·ν_E]_E[∇χ_1·ν_E]_E+[∇ v_·ν_E]_E[∇χ_2·ν_E]_E)
                       ≲(h_^-1[∇ u_·ν_E]__L^2(Γ)^2+h_ ^-1[∇ v_·ν_E]__L^2(Γ)^2)^1/2.

	
	
	
	
	
	
	
	
	
	
	The definition of B_(∙,∙,∙) yields
	
    B_(Ψ_,Ψ_,χ)
    	   = -[u_,v_]χ_1+[u_,u_]χ_2.

	
	
	
	
    The Cauchy-Schwarz inequality and Lemma <ref> show
	
    B_(Ψ_,Ψ_,χ)-L_(χ)≲∑_K∈ h_K^4(f+[u_,v_]_L^2(K)^2+[u_,u_]_L^2(K)^2).

	
	
	A substitution of the previous six estimates in (<ref>) and Lemma <ref> establish (<ref>).


The efficiency of the error estimator obtained in Theorem <ref> can be established using the standard bubble function techniques <cit.>; further details are therefore omitted.





	Let Ψ=(u,v)∈ be a regular solution to (<ref>) and let Ψ_=(u_,v_)∈ be the local solution to (<ref>). There exists a positive constant C_ eff independent of h but dependent on Ψ such that
	
    ∑_K∈η_K^2+∑_E∈(Ω)η_E^2≤ C_ eff^2(Ψ-Ψ__^2+ osc^2(f)),

where osc^2(f):=∑_K∈h_K^4f-f_h_L^2(K)^2 and f_h denotes the piecewise average of f.




§ A C^0 INTERIOR PENALTY METHOD


The analysis of this paper extends to a C^0 interior penalty method for the von Kármán equations formally for σ_1→∞ when σ_1 disappears but the trial and test functions become continuous. The novel scheme is the above dG method but with ansatz test function restricted to ∩^1_0(Ω)=:^2_0()≡ S^2_0()× S^2_0() and the norm ∙_ is ∙_ with restriction to ^2_0() excludes σ_1 (which has no meaning as it is multiplied by zero) and ∙_ is ∙_h with restriction to ^2_0().

Since the discrete functions are globally continuous for this case, the bilinear form  a_(∙,∙) simplifies for some positive penalty parameter σ_2, for η_,χ_∈ S^2_0(),  to

    a_(η_,χ_):= D^2η_:D^2χ_-⟨ D^2η_ν_E⟩_E  ·[∇χ_]_E
           -⟨ D^2χ_ν_E⟩_E  ·[∇η_]_E+∑_E∈σ_2/h_E∫_E[∇η_·ν_E]_E[∇χ_·ν_E]_E.


This novel C^0 interior penalty (C^0-IP) method for the von Kármán equations seeks u_,v_∈ S^2_0() such that

    a_(u_,φ_1)+b_(u_,v_,φ_1)+b_(v_,u_,φ_1)=l_(φ_1)φ_1∈ S^2_0(), 
        a_(v_,φ_2)-b_(u_,u_,φ_2)=0φ_2∈ S^2_0().

The term related to jump which is of the form [η_]_E for each η_∈ S^2_0() vanishes in the C^0-IP method and this simplifies the analysis. 

	Let  Ψ be a regular solution to (<ref>)  and let Ψ_=(u_,v_) be the solution to (<ref>)-(<ref>). For sufficiently large σ_2 and sufficiently small h, it holds
	
    Ψ-Ψ__≤ C h^α.




The Lemmas <ref>-<ref> hold as it is and the boundedness results in Lemma <ref> for b_(∙,∙,∙) can be modified to

    b_(η,χ,φ)≲η_χ_φ_η,χ,φ∈ X+S^2_0(),
    η_2+αχ_φ_1η∈ X∩ H^2+α(Ω) and χ,φ∈ X+S^2_0().

Theorems <ref>-<ref> follow in the same lines and hence, a priori error estimates in energy norm can be established without any additional difficulty.	

For K ∈ and E ∈(Ω), a posteriori error estimates for the C^0-interior penalty method (<ref>)-(<ref>) lead to
the volume estimator η_K and the edge estimator η_E defined by

    η_K^2   := h_K^4(f+[u_,v_]_L^2(K)^2+[u_,u_]_L^2(K)^2),
    η_E^2   :=h_E([D^2 u_ ν_E]_E_L^2(E)^2+[ D^2 v_ν_E]_E_L^2(E)^2)
    
           +h_E^-1([∇ u_]_L^2(E)^2+[∇ v_]_L^2(E)^2).

 Let Ψ=(u,v)∈ be a regular solution to (<ref>) and Ψ_=(u_,v_)∈ S^2_0()× S^2_0() be the solution to (<ref>)-(<ref>). For sufficiently large σ_2 and sufficiently small h, there exist h-independent positive constants C_ rel and C_ eff such that
	
    C_ rel^-2Ψ-Ψ__^2   ≤∑_K∈η_K^2+∑_E∈(Ω)η_E^2≤ C_ eff^2Ψ-Ψ__^2+ osc^2(f).



	This follows as in  Theorems <ref>-<ref>, and hence the details are omitted for brevity.


	The C^0-IP formulation of <cit.> chooses the trilinear form b_(∙,∙,∙) as
	
    b_(η_,χ_,φ_)   :=- [η_,χ_]φ_
           +∑_E∈(Ω)∫_E
    	[⟨(D^2η_)⟩_E  ∇χ_·ν_E]_Eφ_

		for all η_,χ_,φ_∈ S^2_0().
	For the C^0-IP formulation (<ref>)-(<ref>) with b_(∙,∙,∙) replacing b_(∙,∙,∙) and ∙_h≡∙_, the efficiency of the estimator is still open, for difficulties caused by the the non-residual type average term ⟨(D^2η_)⟩_E.




§ NUMERICAL EXPERIMENTS

This section is devoted to numerical experiments to investigate the practical parameter choice and adaptive mesh-refinements.



 §.§ Preliminaries

The discrete solution to (<ref>) is obtained using the Newton method defined in (<ref>) with initial guess Ψ_^0∈ computed as the solution of the biharmonic part of the von Kármán equations, i.e., Ψ_^0∈ solves

    A_(Ψ_^0,Φ_)=L(Φ_)Φ_∈.

Let the ℓ-th level error (for example, in the norm Ψ-Ψ__) and the number of degrees of freedom (ndof) be denoted by e_ℓ and (ℓ), respectively. The ℓ-th level empirical rate of convergence is defined by

    (ℓ):=2×log(e_ℓ-1/e_ℓ)/log((ℓ)/(ℓ-1) ) for ℓ=1,2,3,…

In all the numerical tests, the Newton iterates converge within 4 steps with the stopping criteria Ψ_^5-Ψ_^j-1_<10^-8 for j∈, where Ψ_^5 denotes the discrete solution generated by Newton iterates at 5-th iteration. The penalty parameters for the DGFEM and C^0-IP are consistently chosen 
as σ_1=σ_2=20 in all numerical examples and appear as sensitive as in the case of the linear biharmonic equations. 



 §.§ Example on unit square domain

The exact solution to (<ref>)  is u(x,y)=x^2y^2(1-x)^2(1-y)^2 and 
v(x,y)=sin^2(π x)sin^2(π y) on the unit square Ω
with elliptic regularity index α=1 and corresponding data f and g. 
Figure <ref> displays  the initial mesh and its successive 
red-refinements lead to sequence of  DGFEM solutions on the quasi-uniform meshes 
with the errors u-u__ and v-v__ and with their
empirical convergence rates in Table <ref>. The empirical convergence rate  with respect to dG norm is one as predicted from Theorem <ref>. 
Table <ref> displays the errors  and empirical rates of convergence for C^0-IP method of Section <ref> and the linear empirical rate of convergence as expected from the theory is observed.
 

		
							





 §.§ Example on L-shaped domain
 
In polar coordinates centered at the re-entering corner of the  
 L-shaped domain Ω=(-1,1)^2 ∖([0,1)×(-1,0]),
 the slightly singular functions 
u(r,θ)=v(r,θ):=(1-r^2 cos^2θ)^2 (1-r^2 sin^2θ)^2 r^1+αg_α,ω(θ) 	with the abbreviation 
	
g_α,ω(θ):=
	
    (1/α-1sin((α-1)ω)-1/α+1sin((α+1)ω))×(cos((α-1)θ)-cos((α+1)θ))
       -(1/α-1sin((α-1)θ)-1/α+1sin((α+1)θ))×(cos((α-1)ω)-cos((α+1)ω)),

are defined for the angle ω=3π/2 and the parameter
α= 0.5444837367 as the non-characteristic 
	root of sin^2(αω) = α^2sin^2(ω).
With the loads f and g according to (<ref>)
the DGFEM solutions are computed on a sequence of quasi-uniform meshes. 

Tables <ref> and  <ref>  
display the errors and the expected reduced
empirical convergence rates for the DGFE and the C^0-IP method. 















 §.§ Adaptive mesh-refinement

For the  L-shaped domain of the preceding Example <ref>

the constant load function f ≡ 1, the unknown solution 
to the (<ref>) is approximated by an adaptive mesh-refining algorithm. 

Given an initial  triangulation _0 run the steps
SOLVE, ESTIMATE, MARK and REFINE successively 
for different levels ℓ=0,1,2,….

SOLVE Compute the solution of DGFEM Ψ_ℓ:=Ψ_ (resp. C^0-IP Ψ_ℓ:=Ψ_ ) with respect to _ℓ and number of degrees of freedom given by .


ESTIMATE Compute local contribution of the error estimator from (<ref>) (resp. from (<ref>)) 

    η^2_ℓ(K):=η_K^2+∑_E∈(K)η_E^2    for all  K∈_ℓ.


MARK The Dörfler marking chooses a minimal subset ℳ_ℓ⊂_ℓ such that 

    0.3 ∑_K∈𝒯_ℓη^2_ℓ(K)≤∑_K∈ℳ_ℓη^2_ℓ(K).


REFINE Compute the closure of ℳ_ℓ and generate a new triangulation _ℓ+1 using newest vertex bisection <cit.>.





Figure <ref>(a) displays the  convergence history of 
the a posteriori  error estimator  as a function of number of degrees of freedom for uniform and adaptive mesh-refinement of DGFEM and C^0-IP method.
 
Figure <ref>(b) depicts the adaptive mesh for C^0-IP method generated by the above adaptive algorithm for level ℓ=19, and it illustrates the adaptive mesh-refinement near the re-entering corner. The suboptimal empirical convergence rate for uniform mesh-refinement is improved to an optimal empirical convergence rate 0.5 via adaptive mesh-refinement.

To show the reliability and efficiency of the estimators for DGFEM and C^0-IP, another test has been performed over L-shaped domain for the example described in Example <ref>. Figure <ref>(a) displays the convergence history of the error and the a posteriori  error estimator as a function of number of degrees of freedom for uniform and adaptive mesh-refinement of DGFEM and C^0-IP method. 
Figure <ref>(b) displays the convergence history of 
the error and the a posteriori  error estimator for uniform and 
adaptive mesh-refinement of C^0-IP method. 
The ratio between error and estimator C_rel  is plotted in 
Figure <ref>(a)-(b) and almost constant as a numerical evidence of the reliability and efficiency of the estimators for DGFEM 
and C^0-IP methods of Theorem <ref>-<ref> and Theorem <ref>.
	









§ CONCLUSIONS

This paper analyzes a discontinuous Galerkin finite element method for the approximation of regular solutions of von Kármán equations. A priori error estimate in energy norm and a posteriori error control which motivates an adaptive mesh-refinement are deduced under the minimal regularity assumption on the exact solution. The analysis suggests a novel C^0-interior penalty method and provides a priori and a posteriori error control for the energy norm.  Moreover, the analysis can be extended to hp discontinuous Galerkin finite element methods with additional jump terms for higher order derivatives of ansatz and trial functions under additional regularity assumptions on the exact solution.








§ ACKNOWLEDGEMENTS

The first and third authors acknowledge the support of National Program on Differential Equations: Theory, Computation & Applications (NPDE-TCA) and Department of Science & Technology (DST) Project No. SR/S4/MS:639/09. The second author would like to thank National Board for Higher Mathematics (NBHM) and IIT Bombay for various financial supports.

amsplain





