Nonconforming Finite Element Discretisation for Semilinear Problems with Trilinear Nonlinearity
Carsten Carstensen Gouranga Mallik Neela Nataraj




















































































































































































Nonconforming Finite Element Discretisation for 
Semilinear Problems with Trilinear Nonlinearity
    Carsten Carstensen[Department of Mathematics, Humboldt-Universität zu 
Berlin, 10099 Berlin, Germany.  Distinguished Visiting Professor, Department of Mathematics, 
Indian institute of Technology Bombay, Powai, Mumbai-400076. Email cc@math.hu-berlin.de]        
· Gouranga Mallik
[Department of Mathematics,
Indian Institute of Science, Bangalore 560012
India. Email gourangam@iisc.ac.in]
· Neela Nataraj
[Department of Mathematics, Indian Institute of Technology Bombay, Powai, Mumbai 400076, India. Email neela@math.iitb.ac.in]

    Accepted XXX. Received YYY; in original form ZZZ
======================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================

 
The Morley finite element method (FEM) is attractive for  semilinear problems 
with the biharmonic operator as a leading term in  the stream function vorticity formulation of 2D Navier-Stokes problem and in the von  Kármán equations. This paper establishes a 
best-approximation a priori error analysis and an a posteriori error analysis of discrete
solutions close to an arbitrary  regular solution on the continuous level to semilinear problems with a trilinear nonlinearity.
The analysis avoids any smallness assumptions on the data and so has to 
provide discrete stability by a perturbation analysis before the Newton-Kantorovic 
theorem can provide the existence of discrete solutions. An abstract framework 
 for the stability analysis in terms of discrete operators from the medius analysis 
 
 leads to 
 new results
 on the nonconforming  Crouzeix-Raviart FEM for 
 second-order linear non-selfadjoint and indefinite elliptic problems 
 with L^∞ coefficients. The  paper  identifies  six parameters 
 and sufficient conditions for the local a priori and  a posteriori error control of 
 conforming and nonconforming  discretisations of a class of semilinear 
 elliptic problems first in an abstract framework and then in the
 two semilinear  applications. This leads 
 to new best-approximation error estimates  
 and to 
  a posteriori error estimates in terms of explicit residual-based 
 error control for the conforming and Morley FEM.


Keywords: nonconforming, Morley finite element, elliptic, semilinear, 
stream function vorticity formulation, 2D Navier-Stokes equations, 
von Kármán equations, a posteriori,
second-order linear non-selfadjoint and indefinite elliptic, Crouzeix-Raviart



§ INTRODUCTION




 §.§ Motivation

The nonconforming finite element methods (FEMs) have recently been rehabilitated
by the medius analysis, which combines arguments from traditional a priori and a posteriori error analysis <cit.>. In particular, nonconforming finite element schemes can be equivalent <cit.> or superior to conforming finite element 
schemes <cit.>. The conforming FEMs for fourth-order problems require C^1 conformity and lead to cumbersome implementations, while the nonconforming Morley FEM 
is as simple as quadratic Lagrange finite elements; the reader may consider the 
finite element program in 
<cit.> with less than 30 lines of Matlab for a proof of its simplicity. The second-best scheme of easy implementations for  fourth-order problems is the 
C^0 interior penalty method (C0IP) <cit.> with the benefit of higher-order 
variants and the disadvantage of a (critical) stability parameter choice.  The optimal convergence rates 
are known for the  adaptive Morley FEM <cit.> in fourth-order problems, 
but open for C0IP;  cf.
<cit.>  for the state of the art in second-order applications.  
Hence the advantage of higher-order schemes is not guaranteed for C0IP and leaves 
the Morley FEM as the method of choice. 



 This relevance of the nonconforming Morley FEM  for 
fourth-order problems is not reflected in the contributions  in the literature on the attractive application 

to semilinear problems with the linear biharmonic operator as the leading term 
 (plus quadratic lower-order contributions).  
There are important  model applications of this problem in the stream-function formulation of the incompressible 2D Navier-Stokes equations <cit.> 
and in the von Kármán equations <cit.> for 
nonlinear plates in solid mechanics.  




This paper enriches the general theory of semilinear problems for trilinear low-order terms from conforming FEMs <cit.> to nonconforming FEMs with the medius analysis. This overcomes the smallness assumption (on the load f) in  <cit.> and adds 
 a posteriori error control beyond   <cit.> for a dG discretisation.
The Morley FEM allows for additional benefits and leads, for instance, to guaranteed
lower eigenvalue bounds <cit.>. 




 §.§ Discrete Stability

This paper considers the local approximation of a general regular solution u
to a nonlinear function N(u) without any
extra conditions. 
The invertible Frechét derivative DN(u) of the  nonlinear function N:X→ Y^* 
at a regular solution u is by definition 
a linear bijection between Banach spaces X and Y^*;  
this is equivalent to  an inf-sup condition on the associated 
bilinear form DN(u;∙,∙)=a+b: X× Y→ (split into two contributions
a and b  in Section 2). For a nonconforming finite element  discretisation with some finite element 
space X_h× Y_h ⊄X× Y, in the 
absence of further conditions,  the  inf-sup condition for 
a+b: X× Y→ does not
imply an inf-sup condition  for the discrete bilinear form
a_h+b_h:X_h× Y_h→.  Section 2 studies two general 
bilinear forms a and b defined on a superspace 
X×Y of X_h× Y_h and X× Y, and introduces 
four parameters in  (H1)-(H4) with a sufficient  condition for an
inf-sup condition  to hold for a_h+b_h: X_h× Y_h→
to enable a Petrov-Galerkin scheme and is the first contribution of this paper. 



There will be three applications of this abstract framework in this paper. The first
of which is on former results in  <cit.>
on a nonconforming Crouzeix-Raviart FEM for well-posed 
second-order linear self-adjoint and indefinite elliptic problems: Since the framework
applies the medius analysis tools, there are no smoothness assumptions and the
feasibility and best-approximation property  for sufficiently small mesh-sizes 
is newly established for the Crouzeix-Raviart FEM for L^∞ coefficients in this paper 
(compared to  piecewise Lipschitz continuous coefficients in  <cit.>).



 §.§ Fourth-order semilinear problems

The second and third applications of this discrete stability framework of Section 2
is on semilinear problems with a trilinear nonlinearity:  The stream function formulation of the incompressible 2D 
Navier-Stokes problem <cit.> in Section 4
and the von Kármán equations <cit.> in Section 5 with conforming and Morley FEM.  
The abstract stability result (a) overcomes the high regularity assumptions 
u ∈ H^2_0(Ω) ∩ H^3(Ω) and (b) is not restricted to small 
data  as in <cit.>.







 §.§ Overview of further results
 
The main abstract results are stated in (A)-(D) below. Here and throughout the paper, it is  assumed that the mesh needs to be sufficiently fine to well approximate the solutions to the linear problems associated with the leading elliptic differential operator. 

Throughout this subsection let N:V→ V^* be a differentiable function in a Hilbert space V with dual V^*
with one fixed  regular solution u to N(u)=0. The Hilbert space is a Sobolev space H_0^m(Ω)
associated to some  polyhedral bounded Lipschitz domain 
Ω⊂ℝ^n that is partitioned by arbitrarily fine 
shape-regular triangulations into simplices. The latter form  a  family  and given  any δ>0,
 let (δ) denote the (nonempty) 
subset of all triangulations  of maximal mesh-size smaller than or equal to δ.
For each  ∈, suppose there is a conforming or nonconforming finite element space V_h() and a 
differentiable function N_h:V_h()→ V^*_h() with additional conditions; in particular, there 
is a norm ∙_V on V+V_h() that extends the norm in V.
(Notice the simplified notation N_h≡ N_h().) 
This paper discusses conditions in  (H1)-(H6) sufficient for the subsequent consequences. 



(A). 
 There exist  ϵ,δ>0 such that, for all
∈(δ), there exists a unique discrete solution u_h∈ 
V_h() to N_h(u_h)=0 with u-u_h_V≤ϵ.

(B). 
 There exist ϵ,δ,ρ>0 such that  (A) holds and, 
for all ∈(δ) and  for any  initial iterate  
u_h^(0)∈  V_h() with u_h-u_h^(0)_V≤ρ,
 the Newton scheme converges
quadratically to  u_h .

(C).  There exist ϵ,δ,C_ qo>0 such that (A) holds and, for all 
∈(δ),

    u-u_h_V≤ 
    C_ qo(min_v_h∈  V_h() u-v_h_V+ apx())

with some approximation term apx() to be specified in the particular 
application.

A local reliable and efficient
a posteriori error control holds even for inexact solve (owing to a 
termination in an iterative solver)  in the   sense of

(D). There exist   ϵ,δ,C_ rel,C_ eff>0 such that  
any approximation   v_h∈ V_h() 
with u-v_h_V≤ϵ and ∈(δ)
satisfies

    C_ rel^-1 u-v_h_V≤ N(v_h)_V^* + min_v∈ V v_h-v _V≤ C_ eff u-v_h_V.



It is part of the abstract results in Section 2 and 3 to identify the reliability and efficiency constants in the above displayed estimate and prove that the positive constants 
ϵ, δ, ρ, C_ qo, C_ rel, and 
C_ eff are  mesh-independent. 

The abstract error control in (D) is the point of departure in the applications to the  stream function formulation of the incompressible 2D Navier-Stokes problem <cit.> in Section 4
and the von Kármán equations <cit.> in Section 5.  
This paper  establishes the first  reliable  estimate of N(v_h)_V^* and 
min_v∈ V v_h-v _V in terms of an explicit residual-based
error estimator for the conforming and Morley FEM and discusses its efficiency. 




 §.§ Outlook



This  presentation is restricted to quadratic problems in which the weak formulation involves  a trilinear form for a simple outline to cover  two important semilinear fourth-order problems. The generalisation to more general and stronger nonlinearities, however, requires appropriate  growth conditions in various norms and involves a more technical framework. The presentation matches exactly the nonconforming applications 
(Crouzeix-Raviart and Morley finite elements); other schemes
like  the discontinuous Galerkin schemes <cit.> with their discrete norms and various jump conditions could be included with more additional technicalities.

The ideas developed in this paper extend to other semilinear problems,
 to optimal control and obstacle problems <cit.>  (governed by fourth-order plates and very thin plates), fully nonlinear Monge-Ampere equations based on vanishing moment method <cit.>. 
Moreover, the  lowest-order version of   skeletal or polytopal, hybridizable  discontinuous Galerkin and   
higher-order hybrid methods   <cit.>  is a perturbation of the
nonconforming Crouzeix-Raviart finite element method for the Poisson problems. It is therefore expected that
the  Morley FEM is related to the lowest-order variant of skeletal schemes  for PDEs governed by fourth-order elliptic equations <cit.>. In this way this paper stimulates the development  of 
the a priori and a posteriori error analysis of  those schemes.




 §.§ General notation

Standard notation on Lebesgue and Sobolev spaces applies throughout the paper and ∙ abbreviates ∙_L^2(Ω) with a L^2 scalar product (∙,∙)_L^2(Ω), while the duality brackets <∙,∙>_V^*× V 
are reserved for a dual pairing in V^*× V; 
∙_∞ abbreviates the norm in L^∞(Ω); 
H^m(Ω) denotes the  Sobolev  spaces of order m with  
norm ∙_H^m(Ω); H^-1(Ω) 
(resp. H^-2(Ω)) is the dual space of 
H^1_0(Ω):={v∈ H^1(Ω): v|_∂Ω=0} 
(resp. H^2_0(Ω):={v∈ H^2(Ω): v|_∂Ω=∂ v/∂ν|_∂Ω=0}). With a regular triangulation  of the polygonal Lipschitz domain Ω⊂^n into simplices, associate 
its piecewise constant mesh-size h_∈ P_0() with 
h_|_T:=h_T:=diam(T) ≈  |T|^1/n for all T∈ and its 
maximal mesh-size  h_max:=max h_.
Here and throughout, 

    P_k():={v∈ L^2(Ω): ∀ T∈,v|_T∈ P_k(T)}

denotes the piecewise polynomials of degree at most k∈ℕ_0
and let  Π_k denote the  L^2(Ω) (resp. L^2(Ω;^n) or 
L^2(Ω;^n× n)) orthogonal projection onto P_k() 
(resp. P_k(;^m) or P_k(;^m × m)). Oscillations of 
degree k read 

    osc_k(∙,):=h_^p(I-Π_k)∙_L^2(Ω)

with its square osc_k^2(∙,):= osc_k(∙,)^2 for 
p=1 for second-order in Section 2
and p=2 for fourth-order problems in Sections 4 and 5. 
The notation A≲ B means there exists a generic 
h_-independent constant C such that A≤ CB; A≈ B abbreviates A≲ B≲ A. In the sequel, C_rel and C_eff denote generic reliability and efficiency constants. The set of all n× n real symmetric matrices is :=^n× n_sym. 



§ WELL-POSEDNESS OF THE DISCRETE PROBLEM

This section presents sufficient conditions for the stability of nonconforming discretizations of a well-posed linear problem. Subsection 2.1 introduces four parameters (H1)-(H4) and a condition on them  sufficient for a discrete inf-sup condition for the sum a+b of two bilinear forms a, b: X× Y→ extended to superspaces X⊃ X+X_h and  
Y⊃ Y+Y_h. Subsection 2.2 discusses a first application to second-order non-self adjoint and indefinite elliptic problems <cit.>. 



 §.§ Abstract discrete inf-sup condition

Let X (resp. Y) be a real Banach space with norm ∙_X (resp. ∙_Y) and suppose X and X_h (resp. Y and Y_h) are two complete linear subspaces of X (resp. Y) with inherited norms ∙_X:=(∙_X)|_X and ∙_X_h:=(∙_X)|_X_h (resp. ∙_Y:=(∙_Y)|_Y and ∙_Y_h:=(∙_Y)|_Y_h). Let a,b:X×Y→ be bounded bilinear forms and abbreviate

    a:=a|_X× Y,  a_h:=a|_X_h× Y_h and  b:=b|_X× Y,  b_h:=b|_X_h× Y_h.


Let the bilinear forms a and b be associated to the linear operators A and 
B∈ L(X;Y^*), e.g., Ax:=a(x,∙)∈ Y^* for all x∈ X. Suppose that the linear operator   
A∈ L(X; Y^*) (resp. 
A+B∈ L(X;Y^*))  associated to the  
bilinear form a 
(resp.  a+b)  is invertible and 

    0<α   :=inf_x∈X
    x_X=1sup_y∈Y
    y_Y=1a(x,y);
    
    0<β   :=inf_x∈ X
    x_X=1sup_y∈ Y
    y_Y=1(a+b)(x,y).

Suppose that three linear operators P∈ L(Y;Y_h), Q∈ L(X_h; X), 
∈ L(Y_h;Y) exist and  lead to parameters 
δ_1,δ_2,δ_3,Λ_4≥ 0 in 


  (H1)
δ_1:=sup_x_h∈ X_h
x_h_X_h=1sup_y_h∈ Y_h
y_h_Y_h=1a(A^-1(b(x_h,∙)|_Y),y_h- y_h);

  (H2)
δ_2:=sup_x_h∈ X_h
x_h_X_h=1sup_y∈Y
y_Y=1a(x_h + A^-1(b(x_h,∙)|_Y),y-Py);

  (H3)
δ_3:=sup_x_h∈ X_h
x_h_X_h=1b(x_h,(1-)∙)_ Y_h^*;

  (H4) ∃Λ_4<∞ ∀ x_h∈ X_h  (1-Q)x_h_X≤Λ_4  dist_∙_X(x_h,X).

Abbreviate the bound 
b_X× Y^* 
 of the bilinear form b|_X× Y^*
 simply by  b
 and  set a:=A_L(X;Y^*)  as well as 
 A^-1:=A^-1_L(Y^*;X) — whenever there is no risk of confusion 
 (e.g. with the L^2 norm ∙ of a Lebesgue function). 
If  (H4) holds with 0≤Λ_4<∞, set

    β:=β/Λ_4β+a(1+Λ_4(1+A^-1b))>0.

In the applications discussed in this paper, δ_1+δ_2+δ_3 
from (H1)-(H3) will be smaller than αβ so that the subsequent result provides a discrete inf-sup condition with β_h>0.


Under the aforementioned notation, (<ref>)-(<ref>) 
and (H1)-(H4) imply 

    αβ-(δ_1+δ_2+δ_3)≤β_h:=inf_x_h∈ X_h
     x_h_X_h=1sup_y_h∈ Y_h
     y_h_Y_h=1 
    (a_h+b_h)(x_h,y_h).




Given any x_h∈ X_h with x_h_X_h=1, define 
    x:=Qx_h, ξ=A^-1(b(x_h,∙)|_Y)∈ X, and η=A^-1(b(x,∙)|_Y)∈ X .
 

The inf-sup condition (<ref>) and Aη=Bx lead to

    βx_X≤Ax+Bx_Y^*=A(x+η)_Y^*≤ax+η_X.

This and triangle inequalities imply

    β/a x_X≤x+η_X≤x-x_h_X+x_h+ξ_X+η-ξ_X.

The definition of ξ and η, the boundedness of the operator A^-1 and 
of the bilinear form b|_X× Y show

    ξ-η_X=A^-1(b(x-x_h,∙)|_Y)_X≤A^-1bx-x_h_X.

The combination of (<ref>)-(<ref>) reads

    β/a x_X≤x_h+ξ_X+(1+A^-1b)x-x_h_X.

Since (H4) implies 

    x-x_h_X≤Λ_4 x_h+ξ_X,

the estimate (<ref>) results in 

    x_X≤a/β(1+Λ_4(1+A^-1b))x_h+ξ_X.

The triangle inequality and (<ref>)-(<ref>) lead to

    1=x_h_X_h   ≤x-x_h_X+x_X
       ≤(Λ_4+a/β(1+Λ_4(1+A^-1b)))x_h+ξ_X.

With the definition of β in (<ref>), this reads

    β≤x_h+ξ_X.


For given x_h+ξ∈X and for any 0<ϵ<α, the inf-sup condition (<ref>) implies the existence of some y∈Y with y_Y=1 and

    (α-ϵ)x_h+ξ_X≤a(x_h+ξ,y)=a(x_h+ξ,y-Py)
    +a(x_h+ξ,Py).

Since a(ξ, y_h)=b(x_h, y_h) for y_h:=Py, the latter term is equal to

    a(x_h+ξ,y_h)=a_h(x_h,y_h)+b_h(x_h,y_h)+a(ξ,y_h- y_h)+b(x_h, y_h-y_h).

Let γ_h:=a_h(x_h,y_h)+b_h(x_h,y_h), then (H1)-(H3) and (<ref>) lead to

    a(x_h+ξ,y)≤γ_h+δ_1+δ_2+δ_3.

The combination of (<ref>)-(<ref>) and 
ϵ↘0 in the end result in

    αβ-(δ_1+δ_2+δ_3)≤γ_h≤a_h(x_h,∙)+b_h(x_h,∙)_Y_h^*.

The last estimate holds for an arbitrary x_h with x_h_X_h=1 and so proves the discrete inf-sup condition  αβ-(δ_1+δ_2+δ_3)≤β_h. 


It is well known that  a  positive β_h>0 in  (<ref>) implies the 
best-approximation for the Petrov-Galerkin scheme 
<cit.>
in the following sense.

  
Suppose (X,a) is a  Hilbert space  and 
u∈ X, u_h∈ X_h, and  F∈Y^* satisfy 
(a+b)(u,∙)=F:=F|_Y∈ Y^* and 
(a_h+b_h)(u_h,∙)=F_h:=F|_Y_h∈ Y_h^*. Then  

    β_h   u- u_h _X≤ M min_x_h∈ X_h u - x_h _X
    +
    sup_y_h∈ Y_h
     y_h_Y_h=1(F_h(y_h)- (a+b)( u ,y_h))

with the bound M:=  a+b_X× Y_h ≤a+b of the bilinear form 
(a+b)|_X× Y_h.


The proof of the quasi-optimal convergence 
for a stable discretisation is nowadays standard  in all finite element textbooks 
in the context of the Strang-Fix lemmas.  





 §.§ Second-order linear non-selfadjoint 
  and indefinite elliptic problems


This subsection applies (H1)-(H4) to second-order linear self-adjoint and indefinite elliptic problems and establishes a priori estimates for conforming and nonconforming FEMs under more general conditions on the smoothness of the coefficients of the elliptic operator and for Ω⊂ℝ^n  vis-à-vis <cit.>.



  §.§.§ Mathematical model

The strong form of a second-order problem with L^∞ coefficients A,  b and γ  
reads: Given f∈ L^2(Ω) seek u∈ V:=H^1_0(Ω) such that 

    ℒu :=  -∇· ( A∇ u+u 𝐛)+ γ  u=f.

The coefficients A∈ L^∞(Ω;), 𝐛∈ L^∞(Ω;^n),γ∈ L^∞(Ω) 
satisfy 0<λ≤λ_1( A(x))≤⋯≤λ_n( A(x))
≤λ<∞ for the eigenvalues λ_j( A(x)) 
of the SPD A(x)  for a.e. x∈Ω. 

For u,v∈ V, the expression

    a(u,v):=∫_Ω( A∇ u)·∇ v

defines a scalar product on V (and V is endowed with this scalar product in the sequel)
equivalent to the standard scalar product in the sense  that 
the H^1-seminorm |∙|_H^1(Ω):=∇∙ in V satisfies

    λ^1/2| ∙|_H^1(Ω)≤∙_a
    :=a(∙,∙)^1/2≤λ^1/2| ∙|_H^1(Ω).


Given the bilinear form b:V× V→ with

    b(u,v)   :=∫_Ω(u 𝐛·∇ v+γ uv)  for   all   u,v ∈ V

and the linear form F∈ L^2(Ω)^*⊂ H^-1(Ω)=: V^* with F(v):=∫_Ωfv for all v∈ V, the weak formulation of (<ref>) seeks the solution u∈ V to

    (a+b)(u,v):=a(u,v)+b(u,v)=F(v)   for all  v∈ V.

In the absence of further conditions on the smoothness of the coefficients, any 
higher regularity of the weak solution u∈ H^1_0(Ω) (<ref>)
in the form  u∈ H^s(Ω) for any s>1 
 is not guaranteed even for  f∈ C^∞(Ω)
<cit.>.



  §.§.§ Triangulations

Throughout this paper,  is a set of shape-regular triangulations of the polyhedral 
bounded Lipschitz domain Ω⊂^n into simplices.
Given an initial  triangulation _0 of 
Ω, let
 the newest-vertex bisection define local  mesh-refining that leads to a set of  
shape-regular triangulations ∈. 

Shape-regularity means that there exists  a universal constant 
κ>0  such that the maximal diameter diam(B) 
of a ball B⊂ K satisfies 
κ  h_K≤ diam(B)≤ diam(K)=:h_K for any K∈∈.

Given ∈, let h_∈ P_0() be piecewise constant with h_|_K=h_K
= diam (K) for K∈ and let h_max :=h_max():=max h_;  recall  
(δ):={∈ :  h_max()≤δ} for any  δ >0.

The set of all sides of the shape-regular triangulation  of 
Ω into simplices is denoted by .
The set of all internal vertices (resp. boundary vertices) and  interior sides (resp.  boundary sides) 
of  are denoted by (Ω) (resp.  (∂Ω)) and 
(Ω) (resp. (∂Ω)).



  §.§.§ Conforming FEM

Let P_1() denote the piecewise affine functions in L^∞(Ω) with respect to the triangulation  so that the associated P_1 conforming finite element function spaces without and with (homogeneous) boundary conditions read

    S^1():=P_1()∩ C(Ω̅) and  S^1_0():={v_C∈ S^1():v_C=0  on ∂Ω}.

The interior nodes (Ω) label the nodal basis functions φ_z with patch ω_z={φ_z>0}= int(  supp φ_z) around z∈(Ω).

Given some finite-dimensional finite element space V_h with 
S_0^1()⊆ V_h⊂ V≡ H^1_0(Ω), the discrete formulation 
of (<ref>) seeks u_h∈ V_h with 

    a(u_h,v_h)+b(u_h,v_h)=F(v_h) for all v_h∈ V_h .


The arguments of <cit.>  are rephrased in the following lemma
(proven in the appendix) 
that allows the application of Theorem <ref>  in the subsequent 
theorem. 



For any ϵ>0 there exists some δ>0 such that 







the solution z∈ V≡ H^1_0(Ω) to a(z,∙)=g∈⊂ H^-1(Ω) 
for  g∈  satisfies, for all ∈(δ), that 

    min_z_C∈ S^1_0() z-z_C_a
    + min_Q_0∈ P_0(;^n) A∇ z - Q_0  ≤ϵg.





Adopt the aforementioned assumptions on a and b  in 
(<ref>) -(<ref>) and suppose that (<ref>) 
is well-posed in the sense that it allows for a unique solution u 
for all right-hand sides f∈ L^2(Ω). Then 

    0<β:=inf_x∈ V
     x_h_a=1sup_y∈ V 
     y_a=1(a+b)(x,y)

and for any positive  β_0<β, there exist δ>0 such that

    β_0≤β_h:=inf_x_h∈ V_h
     x_h_a=1sup_y_h∈ V_h 
     y_h_a=1(a+b)(x_h,y_h)

holds for all S_0^1()⊂ V_h:=X_h=Y_h⊂ V with respect to ∈(δ). Moreover,  the solution u to (<ref>) and u_h to  (<ref>) satisfy 

    u-u_h _a≤a+b/β_0min_v_h∈ V_h u-v_h_a.





The invertibility of a linear operator from one Banach space into the dual 
of another is equivalent to an inf-sup condition 
<cit.>; in particular, 
the well-posedness of the theorem implies β>0. The remaining assertions 
follow from Theorem <ref> with a=a, b=b, 
S_0^1()⊂ V_h=X_h=Y_h ⊂ 
X=Y=V=H^1_0(Ω) endowed with the norm ∙_a. 
Then  α=α=1=a and β is the constant in  (<ref>).

To conclude the discrete inf-sup condition, it is sufficient to verify that the parameters involved in (H1)-(H4) can be chosen such that the discrete inf-sup constant in Theorem <ref> is positive.  Moreover, the discrete inf-sup constants of a+b are equal to those of the dual problem a+ b^* with b^*(u,v)=b(v,u). Therefore,
Theorem <ref> is applied to a and b^* (rather than a and b).

 Let Q and  be the identity, while P∈ L(V;V_h) denotes
  the Galerkin projection onto V_h with respect to a, i.e. a(v-Pv,∙)=0 in V_h for all v∈ V. Then the parameters in  (H1), (H3), and (H4) are δ_1=δ_3=Λ_4=0. 
The choice of the parameter δ_2 in (H2) concerns v∈ V and 
u_h∈ V_h with v_a=1= u_h_a and the solution 
z:=A^-1(b^*(u_h,∙))∈ V to a(z,∙)=b(∙, u_h). 
Notice that g:= b·∇ u_h+γ u_h∈ satisfies
 
    b(φ, u_h)=∫_Ω( b·∇ u_h+γ u_h)φ= gφ   for   all  φ∈ V

and (with the Friedrichs constant C_F for ∙≤λ^-1/2C_F∙_a in V)

    g≤ u_h_a   A^-1/2 b_∞+γ_∞u_h≤   ( b_∞+C_Fγ_∞)λ^-1/2
    =: C
    .

The Galerkin orthogonality with  P, the definition of z,  and a Cauchy 
inequality with v_a=1 in the end show

    a(u_h+A^-1(b^*(u_h,∙)),v- P v)=a(z,v- P v)=a(z- P z,v) ≤z-Pz_a.

Given any ϵ>0, Lemma <ref> leads to  δ >0 such that
for all ∈(δ) there exists some  z_C ∈ S_0^1()  with 

    z-Pz_a≤ z-z_C_a ≤ϵg≤ϵ C

with  (<ref>)  in the last step. The combination of the previous inequalities 
proves  (H2) with  δ_2 := ϵ C.

Theorem <ref> applies with β=β and
β_h≥β- ϵ C. This proves the assertion  on β_h≥β_0
for sufficiently small ϵ and δ. 

The quasi-optimal convergence (<ref>) follows from  
Corollary <ref> without the second term in the conforming
discretisation.



The proof requires that the discrete space V_h solely satisfies S_0^1()⊂ V_h⊂ H^1_0(Ω) and so allows for conforming hp finite element  spaces. 
The condition ∈(δ) allows for local mesh-refining as long as
max h_0 is sufficiently small.





  §.§.§ Nonconforming FEM

This subsection establishes the first best-approximation-type a priori 
error estimate for the  lowest-order nonconforming
FEM in any space dimension ≥ 2 under the 
assumptions on the coefficients of Subsubsection <ref>
as an application of Theorem <ref>.
This  generalises  <cit.> from piecewise Lipschitz
continuous to L^∞ coefficients. 

The nonconforming Crouzeix-Raviart (CR) finite element spaces read

    CR^1(𝒯): ={v ∈ P_1(𝒯): ∀ E ∈ℰ,  v  is continuous at mid(E) },
           CR^1_0(𝒯):={v ∈ CR^1(𝒯):v(mid (E))
    =0   for all E ∈ℰ (∂Ω) }.

Here mid (E) denotes the  mid operator for a simplex obtained by taking the arithmetic mean of all vertices.
The CR finite element spaces give rise to the  bilinear forms 
 a_pw,b_pw:CR_0^1()× CR_0^1()→ defined, for all u_CR,v_CR∈ CR_0^1(), by

    a_pw(u_CR,v_CR):=∑_T∈∫_T(𝐀∇ u_CR)·∇ v_CR,
       b_pw(u_CR,v_CR):=∑_T∈∫_T(u_CR𝐛·∇ v_CR+γ u_CRv_CR).

The nonconforming FEM seeks the discrete solution u_CR∈ CR_0^1() to

    a_pw(u_CR,v_CR)+b_pw(u_CR,v_CR)=F(v_CR) v_CR∈ CR_0^1().

Notice that ∇_pw∙ with the piecewise action 
∇_pw of the gradient ∇ is a norm on CR_0^1() and so is 
∙_pw:= A^1/2∇_pw∙.
The subsequent theorem implies the unique solvability and boundedness of discrete solutions for sufficiently fine meshes. 


Suppose that  is a bijection and so ^-1 is bounded and 
(<ref>) 
holds with β=^-1>0. Then 
there exist positive δ and β_0  such that any ∈(δ) satisfies

    β_0≤β_h:=inf_w_ CR∈ CR_0^1()
     w_ CR_ pw=1sup_ v_ CR∈ CR_0^1()
     v_ CR_ pw=1(a_ pw+b_ pw)(w_ CR,v_ CR).




Let H^1():={v∈ L^2(Ω) | ∀ T∈,  v|_T∈ H^1(T)} and endow the vector space

    V:=X:=Y:={v∈ H^1() | ∀ E∈, ∫_E[v]_E=0}⊃ V+CR_0^1()

with the norm
∙_pw:= A^1/2∇_pw∙. 
Here and throughout the paper, the jump of v∈V 
across any interior face 
E=∂ K_+∩∂ K_-∈(Ω) shared by 
two simplices K_+ and K_-  reads

    [v]_E:=
    v|_K_+-v|_K_-  on  E=∂ K_+∩∂ K_-

(then  ω_E:=int( K_+∪ K_-)),  
while [v]_E:=v|_E along any boundary face 
E∈(∂Ω) according to the homogeneous boundary condition on 
∂Ω (and then 
ω_E:=int(K) for K∈ with E∈(K)).

The boundedness of a+b follows from a piecewise Friedrichs inequality

    v̂≤ C_pwF( ∑ _E∈  |  ω_E|^-1   |∫_E [v̂]ds |^2  
    + ∇_pwv̂^2 )^1/2

known for all v̂∈V with the volume |ω_E| of the side-patch
|ω_E|≈ h_E^n. For  v̂∈V and   E∈, 
the integral  ∫_E [v̂]ds=0 vanishes; hence the piecewise Friedrichs 
inequality reduces to v̂≤ C_pwF∇_pwv̂.
This enables a proof that (V,a) is a Hilbert space that 
b is a  bounded bilinear form with respect to those norms. Consequently,
α=1=a and 
 (<ref>)  holds with some β=^-1>0.

Define the nonconforming 
interpolation operator  I_CR∈ L(V;CR_0^1())   by 

    I_CRv    := ∑_E ∈ℰ(_E v ds) ψ_E 
      for all v ∈V

with the side-oriented basis functions ψ_E of CR_0^1() with  
ψ_E( mid(F))=δ_EF, the Kronecker symbol,   for all sides
E,F ∈ℰ.
 For any v_CR∈ CR_0^1(), the conforming companion operator Q:=:= J ∈ L(CR_0^1();V) with  J  v_CR∈ P_4()∩ C^0(Ω̅) from <cit.> satisfies (a) that w:=v_CR-J  v_CR⊥ P_1() 
 is L^2 orthogonal to the space P_1() of piecewise first-order polynomials, 
 (b)  the integral mean property of the gradient

    Π_0(∇_pw(v_CR-J v_CR))=0,

and  (c) the approximation and stability property  (with a universal constant 
Λ_CR) 

    h_^-1(v_CR-J v_CR)
    +∇_pw(  v_CR-J v_CR) ≤Λ_CRmin_v∈ H^1_0(Ω)∇_pw ( v_CR-v) .

(The proofs in <cit.> are in 2D, but can be generalised to any dimension). Note that J is a right inverse to I_CR in the sense that I_CRJ  v_CR=v_CR holds for all v_CR∈ CR_0^1(). The inequality (<ref>) implies (H4)
with Λ_4=(λ/λ)^1/2Λ_CR.

The bilinear forms  
a≡ a_pw, b≡ b_pw:V×V→
read, for all u,v∈V, as 

    a(u,v):=∑_T∈∫_T(𝐀∇u)·∇v  and  b(u,v):=∑_T∈∫_T(u𝐛·∇v +γuv).

As in the stability proof of the conforming FEM, Theorem <ref> applies
to a and  b^* (rather than to a and  b).

The proof of (H1) concerns  u_CR,v_CR∈ CR_0^1() with 
u_CR_pw=1= v_CR _pw and  the solution 
z=A^-1(b(∙, u_CR) |_V)∈ V to 
a(z,∙)= b(∙, u_CR) in V. The right-hand side is the L^2
scalar product of the test function in V with
g:= b·∇_pw u_CR+γ u_CR∈ bounded with 
the discrete Friedrichs inequality ∙≤ C_dF∇_pw∙
in CR_0^1() <cit.> by
 g≤ ( b_∞ +C_dFγ_∞)λ^-1/2=: C_0.
 Since ∇_pw w⊥  P_0(;^n)  in L^2(Ω;^n) for 
 w:= v_CR-Jv_CR, 
 Lemma <ref> applies for any ϵ>0 
 and leads to δ>0 so that, for ∈(δ) 
 with the L^2 projection Π_0,   

    a(A^-1(b^*(u_CR,∙)|_V),v_CR-J v_CR)=
    a_pw(z,w)  =∫_Ω ((1-Π_0) A∇ z )·∇_pww 
         ≤ (1-Π_0) A∇ z  ∇_pww ≤ϵg Λ_CRλ^-1/2≤ C_1ϵ =:δ_1

with (<ref>) for v=0 and ∇_pw v_CR≤λ^-1/2 v_CR_pw =λ^-1/2 
in the end for C_1:=  C_0  Λ_CRλ^-1/2. This 
concludes the proof of  (H1). 

The proof of (H2) concerns  u_CR∈ CR_0^1(), 
 v∈V with u_CR_pw=1
 =v_pw, 
 and   the solution  z∈ V to a(z,∙)= b(∙, u_CR) in V as before.
The operator P:V→ CR^1_0(), however, is not I_CR because 
the oscillating coefficients A prevent the immediate cancellation property for
a(u_CR,v-Pv )=0. The latter  is a consequence of 
the best-approximation P in the Hilbert space V onto its linear and closed subspace  CR_0^1();   so let   Pv∈ CR_0^1() be the unique
minimiser in 

    v -Pv_pw=min_v_CR∈ CR_0^1()v - v_CR_pw≤v_pw =1.

Lemma <ref> applies for any ϵ>0 and leads to δ>0 so that, for 
each ∈(δ), there exists some z_C∈ S^1_0()⊂ CR_0^1() 
with z-z_C_pw≤ϵ C_0. This,  
a(u_CR+z_C,v-Pv )=0,  and 
v- P v_pw≤ 1 in the end 
 provide  (H2) (with b^* replacing b):

    a(u_CR+ A^-1(b^*(u_CR,∙)|_V), 
    v-Pv) 
       =∫_Ω( A∇_pw(z-z_C)·∇_pw(v-Pv)
       ≤ z-z_C_pwv- P v_pw≤    C_0 ϵ  =:δ_2

The proof of (H3) concerns u_CR,v_CR∈ CR_0^1() with u_CR_pw=1= v_CR_pw and w:= v_CR- J  v_CR. 
This and  (<ref>) (with v=0) 
prove

    b^*(u_CR,v_CR- J  v_CR)
    =∫_Ω( b·∇_pwu_CR+γ u_CR) w
       = ∫_Ω g w≤  h_maxgΛ_CR∇_pw v_CR≤ C_0Λ_CRλ^-1/2δ.

Without loss of generality, assume δ≤ϵ. Then 
 (H3) follows with  
 δ_3:= C_3ϵ for 
 C_3:=C_0Λ_CRλ^-1/2.
 (It is remarkable that  in the last inequalities,  the extra property 
 w:=v_CR-J  v_CR⊥ P_1()  leads to the bound 
 λ^-1/2Λ_CRosc_1(g,),
 but  that can easily be exploited solely for piecewise smooth 
 or at least piecewise continuous b and γ).
  
 Since  (H1)-(H4) hold for a and  b^*,  
 Theorem <ref> proves 
 β_h≥β - (C_1+C_2+C_3)ϵ
 with  positive β<β defined in   (<ref>).
Any positive  ϵ<  β/(C_1+C_2+C_3) concludes the proof;
in fact, any  constant β_0 with 0<β_0< β can be realised
in (<ref>)
by small δ>0.


The following best-approximation-type error estimate generalises a result
in  <cit.>.


Let u∈ H^1_0(Ω) solve (<ref>) and set p:= A∇ u+u b∈ H( div,Ω). There exists δ>0 such that for all ∈(δ), the discrete problem (<ref>) has a unique solution u_ CR∈ CR_0^1() and u,   u_ CR,  p and its piecewise integral mean Π_0 p satisfy

    u-u_ CR_ pw   ≲ u-I_ CRu_ pw+p-Π_0p+ osc_1(f-γ u,).


 

Given  e_CR:=I_CRu-u_CR, the discrete inf-sup condition of Theorem <ref> implies the existence of v_CR∈ CR_0^1() with v_CR_pw≤ 1/β_0 and

    e_CR_pw
    = a_pw(e_CR,v_CR)+b_pw(e_CR,v_CR).

Recall from (a)-(b) in the proof of Theorem <ref> that 
v:= J  v_CR satisfies I_CRv=v_CR and Π_1v=Π_1v_CR. 
Since a(u,v)=-b(u,v)+F(v) and u_CR solves (<ref>),  w:=v-v_CR 
satisfies  

    a_pw(e_CR,v_CR)   =a_pw(u,v_CR)-a_pw(u_CR,v_CR)
       =F(w) -a_pw(u,w)-b(u,v)+b_pw(u_CR,v_CR).

This leads in (<ref>) to 

    e_CR   = F(w) -a_pw(u,w)-b_pw(u,w) -b_pw(u-I_CRu, v_CR)
       =∫_Ω(f-γ u) w - ∫_Ω p·∇_pww -b_pw(u-I_CRu, v_CR).

Since ∇_pww   P_0(;^n) in L^2(Ω;^n) and w   P_1(;^n), an upper bound for the first two terms on the right-hand side is

    (I-Π_1)(f-γ u)w - (p-Π_0 p)·∇_pw w
       ≤(p-Π_0 p+ osc_1(f-γ u,))
    λ^-1/2 w_pw
       ≤Λ_CRλ^-1/2β_0^-1(p-Π_0 p+ osc_1(f-γ u,)).

This and a triangle inequality conclude the proof.






















































  §.§.§ A modified CR-FEM for general right-hand sides
 
The nonconforming scheme of the previous subsection allows for a right-hand side in L^2(Ω), 
while conforming variants directly apply to f ∈ H^-1(Ω). This subsection briefly discusses the modification for  f ∈ H^-1(Ω) and a surprising 
analog to Theorem <ref> without  oscillation terms. 

Given  f ∈ H^-1(Ω),  the modified CR-FEM 
seeks  u_CR∈ CR_0^1() such that 

    a_pw(u_CR,v_CR)+b_pw(u_CR,v_CR)   =<f,Jv_CR>_H^-1(Ω) × H^1_0(Ω)  for all  v_CR∈ CR_0^1()

with the duality brackets <∙ ,∙ >_H^-1(Ω) × H^1_0(Ω) on the right-hand side of (<ref>) acting on 
 f ∈ H^-1(Ω) and the test function Jv_CR∈ H^1_0(Ω). The bilinear form b(∙,∙) is replaced in (<ref>) by a modification 
 b_pw(∙,∙) defined, for 
 u_CR,v_CR∈ CR_0^1(),   by 

    b_pw(u_CR,v_CR):=∑_T∈∫_T(u_CR𝐛·∇_pw  v_CR+γ u_CRJv_CR).

The difference to  b_pw (u_CR,v_CR) from  (<ref>) is in the final 
application of Jv_CR rather than v_CR with the conforming companion
operator J from the proof of Theorem <ref>.


Let u∈ H^1_0(Ω) solve (<ref>) with the right-hand side f≡ f∈ H^-1(Ω) and 
 set p:= A∇ u+u b∈ H( div,Ω). There exists δ>0 such that for all ∈(δ), the discrete problem (<ref>) has a unique solution u_CR∈ CR_0^1() 
 and u,  u_ CR,  p and its piecewise integral mean Π_0 p satisfy

    u-u_ CR_ pw≲ u-I_ CRu_ pw+p-Π_0p.




The stability of the modified bilinear form a_pw(∙, ∙) +b_pw(∙, ∙) follows from the methodology of this section. An immediate 
proof follows from the stability (<ref>) and a perturbation argument: For any v_CR, w_CR∈ CR_0^1() with
v_CR_pw=1=  w_CR_pw, 

    | b_pw(v_CR,   w_CR ) - b_pw(v_CR,  w_CR) |
     ≤γ_∞ v_CR    w_CR - J w_CR≤γ_∞ C_pwFΛ_CR    h_max

with a piecewise  Friedrichs inequality and  (<ref>) in the last step. 
The combination with   (<ref>) and a triangle inequality 
prove stability of the modified scheme

    β_0/2≤β_h:=inf_w_CR∈ CR_0^1()
     w_CR_pw=1sup_ v_CR∈ CR_0^1()
     v_CR_pw=1(a_pw+b_pw)(w_CR,v_CR)

for any ∈𝕋(   β_0/(2γ_∞ C_pwFΛ_CR)). 
The proof of the a priori error estimate follows the arguments of the proof of Theorem <ref>.  Given  e_CR:=I_CRu-u_CR, the stability of the modified scheme leads to some v_CR∈ CR_0^1() with norm
v_CR_pw≤ 2/β_0 and 

    e_CR_pw
    = a_pw( I_CRu- u_CR ,v_CR)+b_pw(I_CRu- u_CR,v_CR).

Let  v:=Jv_CR and w:= v-v_CR. Then   (<ref>)  and  (<ref>)  imply  

    e_CR_pw
     =  a_pw( I_CRu-u ,v_CR)  
     -  (p, ∇_pw w)_L^2(Ω) 
    - (u- I_CRu ,   γ v +   b·∇_pw v_CR  )_L^2(Ω) .

The point is that all terms with v disappear and no oscillation terms remain. In fact, all other terms
are controlled by u- I_CRu  _pw or by p-Π_0 p as in the proof
of  Theorem <ref>; further details are omitted.





§ A CLASS OF 
SEMILINEAR PROBLEMS WITH
 TRILINEAR NONLINEARITY
 
This section is devoted to an abstract framework for an a priori and a posteriori analysis to solve a class of semilinear problems that includes the applications in Section 4 and 5.



 §.§ A priori error control

Suppose  X and Y are real Banach spaces and let the 
quadratic function N:X→ Y^* be of the form

    N(x):= x+Γ(x,x,∙)

with a leading linear operator A∈ L(X;Y^*) and  F∈ Y^* for the affine operator x:=Ax-F for all x∈ X and a bounded trilinear form Γ: X× X× Y→.

To approximate a regular u solution to N(u)=0, the discrete version involves some discrete spaces X_h and Y_h plus  a discrete function  F_h∈ Y_h^*, 
_hx_h:=A_hx_h-F_h, and a 
bounded trilinear form Γ_h:X_h× X_h× Y_h→ with N_h(x_h)=_h x_h+Γ_h(x_h,x_h,∙). The discrete problem seeks 
u_h∈ X_h such that

    _hu_h+Γ_h(u_h,u_h,∙)=0   in  Y_h^*.

The a priori error analysis is based on the Newton-Kantorovich theorem and 
adapts  the abstract discrete inf-sup results of Subsection <ref>.
Some further straightforward notation is required for this. Suppose that there exists some   invertible bounded linear operator 
operator A (i.e. Av=a(v,∙) in Y for all v∈X) on extended Banach spaces X and Y and 
suppose that  there exists   a bounded  extension 

    Γ:X×X×Y→  with  Γ 
    :=Γ_X×X×Y:=
    sup_x∈X
    x_X=1sup_ξ∈X
    ξ_X=1sup_y∈Y
    y_Y=1Γ(x,ξ,y)<∞

of Γ=Γ|_X× X× Y with 
Γ_h=Γ|_X_h× X_h× Y_h.
Given the regular solution u∈ X to N(u)=0 in (<ref>), let the 
bilinear form b:X×Y→  
be the linearisation of  Γ at the solution u , i.e.,

    b(∙,∙):= Γ(u,∙, ∙)+Γ(∙,u,∙),

and be bounded by b:=b_X×Y≤ 
2 u_X  Γ. 
Adopt the notation  (<ref>) for the bilinear forms a,a_h,b, and b_h
as respective restrictions of a and b and suppose 
F∈Y^* exists with F:=F|_Y and 
F_h:=F|_Y_h.

Recall 
that the bounded linear operator  A is invertible and the so the associated 
bilinear form a is bounded and satisfies (<ref>) 
with some positive α.

Recall  that u is a regular solution to N(u)=0 in the sense that 
N(u)=0 and DN(u)∈ L(X;Y^*) with DN(u)=(a+b)(∙,∙) satisfies the 
 inf-sup condition (<ref>).
 
Suppose all the aforementioned bilinear forms satisfy (H1)-(H4) with some operators  P∈ L(Y;Y_h), Q∈ L(X_h; X), and ∈ L(Y_h;Y).
In addition to  (H1)-(H4)  suppose 
that δ_5,δ_6≥ 0 satisfy 


(H5)  δ_5:=(F-Au)(1-)∙_Y_h^*;


(H6)  ∃ x_h∈ X_h such that δ_6:=u-x_h_X.


The non-negative  parameters δ_1,δ_2,δ_3,δ_5,δ_6 
and α,  β, b all
depend  on the fixed regular solution u to N(u)=0
and this dependence is suppressed in the notation for simplicity. 

Under the present assumptions and with the additional 
smallness assumption 4 δΓ < β_0
(in the notation of (<ref>)-(<ref>))
the properties 
(A)-(B) hold for the  fixed discretisation at hand 
in the following sense. Suppose that Γ>0 for 
otherwise N is a linear equation with a unique solution
and the results of Section 2 apply. 


Given a regular solution  u∈ X to N(u)=0, assume the existence of 
extended  bilinear forms a and b with 
 (<ref>)-(<ref>) and α>0 (resp. β>0 in  (<ref>) 
 and β>0 in  (<ref>)). 
Suppose that 
(H1)-(H6) hold with parameters δ_1,…,δ_6≥ 0 and 
that  x_h∈ X_h  satisfies  (H6). Suppose that  

    β_0   :=αβ-(δ_1+δ_2+δ_3+2Γδ_6) >0  and 
    δ   :=β_0^-1(δ_5+ aδ_6+δ_6 (x_h_X_h+  u_X)Γ +δ_3/2 )≥ 0

satisfy   4 δΓ < β_0.
Then  ϵ:=δ_6 +δ + r_- with m:=2 Γ/β_0 >0,
h:= δ m≥ 0, 
 
    r_-:= (1-√(1-2h ))/m  - δ≥ 0 ,   and  ρ:=(1+√(1-2h))/ m>0

 satisfy  (i) there exists a solution u_h∈ X_h to N_h(u_h)=0
 with u-u_h _X≤ϵ and (ii)  given any v_h∈ X_h
 with v_h-u_h_X_h≤ρ, the Newton scheme 
 with initial iterate  v_h converges R-quadratically
 to the discrete solution  u_h in (i). 
 If even 4 ϵΓ≤β_0, then (iii)
 there is at most one solution u_h∈ X_h to N_h(u_h)=0
 with u-u_h _X≤ϵ.


The proof is based on the Newton-Kantorovich convergence theorem found, e.g., in <cit.> for X= Y=^n and in <cit.> for Banach spaces. The notation is adopted to the setting of 
Theorem <ref>.

 
Assume the Frechét-derivative DN_h(x_h) of N_h at some x_h∈ X_h satisfies

    D N_h(x_h)^-1_L( Y_h^*; X_h)≤ 1/β_0
          and  D N_h(x_h)^-1N_h(x_h)_X_h≤δ.

Suppose that D N_h is Lipschitz continuous with Lipschitz constant
2 Γ 
and that  4 δΓ≤β_0.
Then there exists a unique root u_h∈ B(x_1,r_-) to  N_h in the ball around the first iterate x_1 := x_h - D N_h(x_h)^-1N_h(x_h) 
  and this is the only root in B(x_h,ρ)
with r_-, ρ from (<ref>). If even 4 δΓ < β_0,
then 
the Newton scheme with   initial iterate x_h leads to a sequence in   B(x_h,ρ)  
that  converges R-quadratically to u_h. 

Proof of Theorem <ref>.
Suppose that δ≥ 0 and Γ>0  so 
that  r_-≥ 0 in (<ref>)  well defined.
The bounded trilinear form 
Γ_h=Γ|_ X_h×  X_h×  Y_h 
leads to the Frechét-derivative DN_h( x_h)∈ L(X_h;Y_h^*) with 

    DN_h( x_h;ξ_h,η_h)=a_h(ξ_h,η_h)+Γ_h(x_h, ξ_h,η_h)
    +Γ_h( ξ_h,x_h,η_h)  for all  x_h, ξ_h∈ X_h, η_h∈ Y_h.

The definitions of a and b and their extensions and discrete versions with 
 (H1)-(H4) allow  in  Theorem <ref> for a positive 
 inf-sup constant 
 β_1:=αβ-(δ_1+δ_2+δ_3)  in (<ref>) for the bilinear form 
 
    DN(u)|_X_h× Y_h= 
     a_h+Γ(u,∙,∙)+Γ(∙,u,∙)
     = a_h+b_h
 
 for the extended nonlinear form N(x)=A(x)
 -F+ Γ(x,x,∙) ∈Y^*
 for x∈X and its derivative  DN(u) at u.
 This discrete inf-sup condition (<ref>)
and  a triangle inequality with x_h from (H6)  lead to an inf-sup constant

    0< β_0:=
    
    β_1 - 2Γδ_6
    ≤β_h:=inf_ξ_h∈ X_h
    ξ_X_h=1sup_η_h∈ Y_h
    η_h_Y_h=1DN_h(x_h;ξ_h,η_h)

for the bilinear form 
DN_h(x_h;∙,∙)=a_h+Γ_h(x_h,∙,∙)
+Γ_h(∙,x_h,∙). 
The discrete inf-sup constant is a singular value and equal to the norm of the inverse operator; 
1/β_0 is an upper bound of the operator norm of the 
discrete inverse. This proves  the first estimate of (<ref>). 
It also proves  in the second estimate of
(<ref>) that

    DN_h(x_h)^-1N_h(x_h)_X_h≤β_0^-1N_h(x_h)_Y_h^*

and it remains to estimate N_h(x_h) in the norm of Y_h^*.
Given any  y_h∈ Y_h with y_h_Y_h=1 and   y:= y_h∈ Y,  an exact Taylor expansion with  N(u;y)=0 shows 

    N_h(x_h;y_h)= N_h(x_h;y_h)-N(u;y)
       = F(y-y_h) +a_h( x_h,y_h)-a(u,y)
    +Γ_h( x_h, x_h,y_h)-Γ(u,u,y)
       =F(y-y_h) -a( u,y-y_h)
    +a(x_h-u,y_h)+Γ_h( x_h, x_h,y_h)
    -Γ(u,u,y).

In  abbreviated duality brackets, the  first two terms in (<ref>)
are equal to

    F(y-y_h) -a( u,y-y_h)
    =⟨F-Au,(-I)y_h⟩≤δ_5

with (H5). The definition of δ_6 in (H6) proves

    a(x_h-u,y_h)≤aδ_6.

Up to the factor 2, the last two terms in (<ref>) are equal to

    2Γ_h( x_h, x_h,y_h)-2Γ(u,u,y)
    =Γ( x_h-u, x_h,y_h)+Γ( x_h,x_h-u,y_h)
         +Γ(u,x_h-u,y)+Γ(x_h-u,u,y)-b(x_h,y-y_h).
       ≤2δ_6(x_h_X_h+ u_X)Γ +δ_3.

The combination of the preceding three displayed estimates with 
(<ref>) implies 

    β_0^-1N_h(x_h)_Y_h^*≤δ

with δ≥ 0 from  (<ref>). The 
combination of (<ref>) and  (<ref>)
shows the second inequality  in (<ref>). The smallness assumption 
reads h<1/2 and is stated  explicitly in the theorem; hence 
the Newton-Kantorovich Theorem <ref> applies. 

Let us interrupt the proof for a brief discussion of the extreme but possible  case 
δ=0 with the implications  δ_6=δ_5=δ_3=0 and 
x_h=u in (H6). The proof of   (<ref>) remains valid 
in this case and then 
N_h(x_h)=0 guarantees  that u=x_h is the discrete solution u_h.
In this very particular situation, the Newton scheme converges  and leads to
the constant sequence x_h=x_1=x_2=... with the  limit x_h=u_h. 
Theorem <ref> applies with r_-=0=ϵ and provides 
(i)-(iii).

Therefore, throughout the remainder of this proof suppose that δ>0 and so 
ρ,ϵ, r_- >0 in  Theorem <ref> show   the existence of a discrete solution  u_h to N_h(u_h)=0 in B(x_1,r_-)
 and this is the only discrete solution in B(x_h,ρ). This and 
 triangle inequalities lead to

    u-u_h_X≤u-x_h_X+
    x_1-x_h_X_h
    +x_1-u_h_X_h≤δ_6+δ+r_-= ϵ

for  the Newton correction x_1-x_h is estimated in the second inequality of 
(<ref>). This proves the existence of a discrete solution 
u_h in  X_h∩B(u,ϵ)  as asserted in (i). 

Theorem <ref> implies (ii) and it remains to prove the uniqueness of discrete solutions in 
B(u,ϵ) under the additional assumption that 
 4 ϵΓ≤β_0, i.e.,  2mϵ≤ 1.
Recall  that the limit u_h∈B(x_1,r_-)  in (i)-(ii) is the only discrete
solution in B(x_h,ρ).
Suppose there exists a second solution 
u_h∈ X_h∩B(u,ϵ) to  N_h(u_h)=0.
The uniqueness in B(x_h,ρ) and a triangle inequality imply that 

    ρ<   x_h- u_h_X≤  u- u_h_X
    +  u- x_h_X≤ϵ+δ_6≤ 2ϵ≤ 1/m

with the smallness assumption on ϵ in the end. But this leads to a
contradiction with 
the  definition of ρ in  (<ref>) and so  concludes the proof of (iii).



In the applications, if h_max is chosen sufficiently small, the parameters δ_1,δ_2,δ_3, δ_5, and  δ_6 are also small.
 In particular, δ from (<ref>) is small and 
 so is ϵ. This ensures 
 4 δΓ≤  4 ϵΓ< β_0
so that  Theorem <ref> applies. 




The convergence speed in the Newton-Kantorovich theorem is known 
to be h=δ m and this parameter is uniformly smaller than one in the applications.
Hence the number of iterations in the Newton scheme does not increase 
as the mesh-size decreases. 





 §.§ Best-approximation

This subsection discusses the best-approximation result (C) for regular solutions of semilinear problems with trilinear nonlinearity under the assumption
(H1)-(H6) with parameters δ_1,…,δ_6 and α (resp. β) from (<ref>)  (resp.  (<ref>)).

The extra term  N̂ (u)_ Y_h^* in the best-approximation result in 
Theorem <ref>  will be discussed afterwards and leads to 
some best-  and data-approximation term.

 
If u is a regular solution to N(u)=0 and 
 δ and  ϵ:=δ_6 +δ + r_- from  (<ref>)-(<ref>)
 satisfy 2mϵ≤ 1, then there exists C_ qo>0 
 such that  the  unique discrete solution 
 u_h∈ X_h∩B(u,ϵ )  satisfies the best-approximation property 

    u-u_h_X≤ C_ qo(min_v_h∈ X_hu-v_h_X
    +N(u)_Y_h^*).




Given the best-approximation u_h^* to u in X_h with respect to the norm in 
X, set e_h:=u_h^*-u_h∈ X_h and apply the discrete inf-sup 
condition  (<ref>) to the bilinear form  
DN(u)|_X_h× Y_h with  the constant 
β_1:=αβ-(δ_1+δ_2+δ_3)
from the proof of Theorem <ref>. 
This leads to  y_h∈ Y_h with y_h_Y_h≤ 1/β_1 and 

    e_h_X_h= DN(u;e_h,y_h).

Since the quadratic Taylor expression of N at u for N_h(u_h;y_h)=0 is exact, e:=u-u_h∈X satisfies

    0=N(u;y_h)-DN(u;e,y_h) - D^2N(u;e,e,y_h).

The sum of (<ref>) and (<ref>),   
D^2N(u;e,e,y_h)=2 Γ(e,e,y_h), and y_h_Y_h≤ 1/β_1
prove

    β_1 e_h_X_h≤N(u)_Y_h^*
    +DN(u)u-u_h^*_X
    +Γe_X^2.

This,  a triangle inequality, and  
min_x_h∈ X_hu-x_h_X  = u-u_h^* _X show

    (β_1-Γe_X)e_X≤N(u)_Y_h^*+(β_1+ DN(u))min_x_h∈ X_hu-x_h_X.

Recall 4ϵΓ≤β_0≤β_1 and 
e_X≤ϵ from Theorem <ref>, so that
3β_1/4  ≤β_1-Γe_X leads 
in (<ref>) to 
C_qo=3/4 max{  1/ β_1,  1 + DN(u)/β_1   }
and  apx():=N(u)_Y_h^* 
in the asserted best-approximation. 
This concludes the proof.


Two examples for the term apx():=N(u)_Y_h^* conclude this subsection. 


If Y_h⊂ Y, then apx()=N(u)_Y_h^*≤N(u)_Y^*=0. 
Hence, Theorem <ref> implies the quasi-optimality of the conforming FEM.



For the second-order linear non-selfadjoint and indefinite elliptic problems of Subsection 2.2,  
Γ=0 and β_0=β_1 etc. is feasible in 
Theorem <ref> and the best-approximation estimate holds. The 
approximation  term apx() is the norm of 
the functional F-(a_pw+b_pw) (u,∙ ) in V_h^*. This is exactly 
the extra  term in Corollary <ref> that leads to the additional
two terms in Theorem <ref>.




 §.§ A posteriori error control

The regular solution u to N(u)=0 is approximated by some v_h∈ X_h sufficiently close to u such that the Theorem <ref> below asserts reliability (<ref>) and efficiency 
(<ref>)-(<ref>).


Any v_h∈ X_h with u-Qv_h_X<β/Γ satisfies

    u-v_h_X≤N(Qv_h)_Y^*/β-Γu-Qv_h_X+Qv_h-v_h_X,
       Qv_h-v_h_X≤Λ_4u-v_h_X,
       N(Qv_h)_Y^*≤ (1+Λ_4)(DN(u)+
    β
    )u-v_h_X.




Abbreviate ξ=Qv_h and e:=u-ξ. Recall that the bilinear form a+b is associated 
to the derivative DN(u;∙,∙)∈ L(X;Y^*) 
with an inf-sup constant β>0. Hence for any 
0<ϵ<β there exists some y∈ Y with y_Y=1 and 

    (β-ϵ)e_X≤ DN(u;e,y).

Since N(u)=0 and N is quadratic, the finite Taylor series

    N(ξ,y)=- DN(u;e,y)+ D^2N(u;e,e,y)

is exact. This, D^2N(u;e,e,y)=2Γ(e,e,y), and (<ref>) imply

    (β-ϵ)e_X   ≤ -N(ξ,y)+Γ(e,e,y)
    ≤N(ξ)_Y^*+Γe_X^2.

With ϵ↘ 0 and β-Γe_X>0, this leads to

    e_X≤N(ξ)_Y^*/β-Γe_X.

A triangle inequality 
u-v_h_X≤e_X+Qv_v-v_h_X 
concludes the proof of (<ref>).

Recall that (H4) implies (<ref>). This and a triangle inequality show

    e_X≤u-v_h_X(1+Λ_4).

The identity (<ref>) results in

    N(ξ)_Y^*≤DN(u;e)_Y^*+Γ(e,e,∙)_Y^*≤(DN(u)+Γe_X)e_X.

The combination of the previous two displayed estimates proves (<ref>).


The discrete function v_h can be estimated in the sense of (D) from the introduction.


In addition to  the assumptions of Theorem <ref> suppose that 
u- v_h_X≤ϵ≤κβ /(Γ (1+Λ_4)) holds
for some positive κ<1 and  v_h∈ X_h. Then 
C_ rel,1:=1/(β(1-κ)) and 
 C_ rel,2:= 1+  LC_ rel,1  for 
 L:=a+2Λ( u _X+ϵ(1+Λ_4))
satisfy reliability in the sense that 

    u- v_h_X≤ C_ rel,1N(v_h)_Y^*
    +C_ rel,2Qv_h-v_h_X

and efficiency with  (<ref>) and with C_ eff,1:= ((1+Λ_4) (DN(u)+β)  + LΛ_4) in 


    N(v_h)_Y^*≤ C_ eff,1  u- v_h_X .
 



Recall the abbreviations  ξ=Qv_h and e:=u-ξ.  A triangle inequality and (H4)  
show that e_X ≤ (1+Λ_4)  u- v_h_X≤ϵ(1+Λ_4)≤κβ/ Γ.
This and Theorem <ref> imply 

    u-v_h_X≤N(Qv_h)_Y^*/β(1-κ)
     +Qv_h-v_h_X.

The derivative DN is globally Lipschitz continuous with a Lipschitz constant
2Λ, the function 
N is  Lipschitz continuous 
 in the closed  ball   B(u,  ϵ(1+Λ_4))  in X 
 with a Lipschitz constant  L. 
Since  v_h, Qv_h∈B(u,  ϵ(1+Λ_4)),

    N(Qv_h)_Y^*≤N(v_h)_Y^* + L Qv_h-v_h_X.

The combination of the previous displayed estimates proves the asserted reliability. The
efficiency employs the Lipschitz continuity as well and then utilises 
(<ref>)-(<ref>)
to verify 

    N(v_h)_Y^*≤N(Qv_h)_Y^*+  L Qv_h-v_h_X≤ C_ eff,1u-v_h_X.

This  concludes the proof. 




§ STREAM FUNCTION 
VORTICITY FORMULATION OF THE
 INCOMPRESSIBLE 2D NAVIER-STOKES PROBLEM

This section is devoted to the stream function vorticity formulation of 2D Navier-Stokes 
equations with right-hand side f∈ in a 
polygonal bounded Lipschitz domain Ω⊂ℝ^2: 
There exists <cit.> at least one 
distributional solution u∈ V:= to

    Δ^2 u+∂/∂ x_1((-Δ u)∂ u/∂ x_2)-∂/∂ x_2((-Δ u)∂ u/∂ x_1)=f in Ω.

The analysis of extreme viscosities lies beyond the scope of this paper and the viscosity 
(the factor in front of the bi-Laplacian in (<ref>))
is set  one throughout this paper.




 §.§ Continuous problem

The weak formulation to (<ref>) seeks u∈ V such that 

    a(u,v)+Γ(u,u, v)=F( v)  v∈ V.

The associated bilinear form a: V× V→ and the trilinear form Γ: V× V× V→ read

    a(η, χ):=ΔηΔχ,   Γ(η,χ,ϕ):=Δη(∂χ/∂ x_2∂ϕ/∂ x_1-∂χ/∂ x_1∂ϕ/∂ x_2),

and F∈ V^* is given by F( ϕ):= f ϕ for all η,χ,ϕ∈ V. 
The Hilbert space V≡ H^2_0(Ω) with the scalar product 
 a(∙,∙) is endowed with the H^2 seminorm 
∙:=|∙|_H^2(Ω)
and ∙_V^* denotes the dual norm.
The bilinear form a(∙,∙) is equivalent to the  scalar product in V
and the trilinear form 
Γ(∙,∙,∙) is bounded (owing to the continuous embedding 
V⊂ H^2(Ω)↪ W^1,4(Ω))  with

    ⟨ N(u), v⟩=N(u; v):=a(u, v)-F(v)+Γ(u,u, v)   for all  u,v∈ V.

The 2D Navier-Stokes equations in the weak stream function vorticity 
formulation (<ref>) seeks u∈ V with N(u)=0.

The regularity results for the biharmonic operator Δ^2  in <cit.> ensure that z∈ V with a(z,∙)∈ H^-1(Ω)⊂ V^* belongs to 
H^2+s(Ω) for some elliptic regularity index  s ∈  (1/2,1]  and 
z_H^2+s(Ω)≤ C a(z,∙)_H^-1(Ω). 
The regularity results for the 
Navier-Stokes problem in <cit.> ensure that 
any weak solution u∈ V to N(u)=0 satisfies u∈ H^2+s(Ω).
This makes the  continuous embeddings  
H^2+s(Ω) ↪ W^1,∞(Ω) (for s>0) and 
H^2+s(Ω)↪ W^2,4(Ω) (for s>1/2) available
throughout this (and the subsequent) section. The embeddings 
and inequalities  imply 
for u∈ H^2+s(Ω) and for θ∈ V, ϕ∈ H^1(Ω) that

    Γ(u,θ,ϕ)
    
    
    
    
       ≲u_H^2+s(Ω)θ_H^2(Ω)ϕ_H^1(Ω).

Consequently,  the derivative b(∙,∙) := DN(u;∙,∙):= 
Γ(u,∙,∙)+Γ(∙,u,∙) at the solution u 
is a bounded bilinear form in H^2(Ω)× H^1(Ω) and will be key in the subsequent analysis.



 §.§ Conforming FEM

Let V_C be a conforming finite element space contained in
C^1(Ω)∩ V; for example, the spaces associated with Bogner-Fox-Schmit, HCT, or Argyris elements <cit.> and a regular triangulation 
of Ω into triangles. The conforming finite element formulation  seeks u_C∈ V_C with

    N_h(u_C;v_C):=N(u_C;v_C):=a(u_C, v_C)-F( v_C)+Γ(u_C,u_C, v_C)=0  v_C∈ V_C.



If u is a regular solution to N(u)=0, then there exist positive 
ϵ, δ, and ρ such that 
 (A)-(C) hold with apx()≡ 0 
for all ∈(δ). 



Set X=Y=V,  X_h=Y_h=V_C,  
a(∙,∙):=a(∙,∙), and 
b(∙,∙):=b(∙,∙)
:=Γ(u,∙,∙)+Γ(∙,u,∙). 
For  and Q chosen as identity, 
the parameters in the hypotheses (H1) and (H3)-(H5) are  
δ_1=δ_3=Λ_4=δ_5=0.

For the proof of (H2), suppose 
 θ_h≡ x_h∈ V_C⊂ V with θ_h =1 and recall
 from the end of the previous subsection that 
 b(θ_h, ∙)∈ H^-1(Ω).  Hence  the solution 
z∈ V to the biharmonic problem

    a(z,ϕ)=b(θ_h,ϕ)ϕ∈ X

satisfies z∈ H^2+s(Ω) and 
z_H^2+s(Ω)≤ C θ_h= C. 
(Note that z is called A^-1(b(x_h,∙)|_Y) in
Subsection <ref>).
This regularity and the Galerkin projection P with the 
Galerkin orthogonality and the approximation 
property z-Pz≲ h_max^s <cit.> lead 
for any y∈ Y≡ V with y =1 to 

    a(x_h + z, y-Py ) =a(z,y-Py )
    = a(z-Pz,y) ≲ h_max^s.

This proves   (H2) with  δ_2≲ h_ max^s.
The choice x_h=Pu implies (H6) with δ_6≲ h_ max^s 
(from the higher regularity of u and 
u-Pu≲ h_max^s). Consequently,   
for sufficiently small maximal mesh-size h_max, 
Theorem <ref> provides 
the discrete inf-sup condition 1≲β_h and 
Theorem <ref> applies.
Since V_C is a conforming finite element space, 
Theorem <ref> holds with 
apx():= N(u)_Y_h^*≡ 0. 
This concludes the proof.


The explicit residual-based a posteriori error estimator for the stream function vorticity formulation of 2D Navier-Stokes equations requires some notation for the differential operators: For any scalar function v, vector field Φ=(ϕ_1,ϕ_2)^T, 
and tensor σ with the 4 entries σ_11, σ_12,
σ_21, and σ_22 in form of a 2× 2 matrix,  





    ∇ v=[ ∂ v/∂ x_1; ∂ v/∂ x_2 ],  Curl  v=[ -∂ v/∂ x_2;  ∂ v/∂ x_1 ],  curl[ ϕ_1; ϕ_2 ]
    =∂ϕ_2/∂ x_1-∂ϕ_1/∂ x_2,
       
    DΦ=
    [ ∂ϕ_1/∂ x_1 ∂ϕ_1/∂ x_2; ∂ϕ_2/∂ x_1 ∂ϕ_2/∂ x_2 ],
 


    [ ϕ_1; ϕ_2 ]=∂ϕ_1/∂ x_1+∂ϕ_2/∂ x_2
    ,  Curl[ ϕ_1; ϕ_2 ]
    =[ -∂ϕ_1/∂ x_2  ∂ϕ_1/∂ x_1; -∂ϕ_2/∂ x_2  ∂ϕ_2/∂ x_1 ]
    ,  and σ
    =[ ∂σ_11/∂ x_1 +∂σ_12/∂ x_2; ∂σ_21/∂ x_1 +∂σ_22/∂ x_2 ].

For any K∈ and E∈(Ω), define the volume and edge  error estimators by 

    η_K^2   :=h_K^4Δ^2 u_C
    - curl(Δ u_C∇ u_C) -f^2_L^2(K), 
    η_E^2   :=h_E^3 [ (D^2   u_C)]_E ·ν_E^2_L^2(E)  
    +h_E[Δ u_C ]_E^2_L^2(E)
         +h_E^3[Δ u_C∇ u_C]_E·τ_E^2_L^2(E)

with   the unit tangential (resp. normal) vector  τ_E 
(resp. ν_E) along  the edge E∈. 
Recall osc_m(∙,):=h_^2(I-Π_m)∙_L^2(Ω) for 
m∈ℕ_0 in all
fourth-order applications.

  
If  u∈ V is a regular solution to N(u)=0 and m∈_0, then there exist positive 
ϵ,δ, C_ rel, and C_ eff such that, for any ∈(δ),
the unique discrete solution  u_C∈  V_C to (<ref>) with 
u-u_C<ϵ satisfies 

    C_ rel^-2 u-u_C^2   ≤∑_K∈η_K^2
    +∑_E∈ (Ω)η_E^2≤ C_ eff^2( u-u_C^2
    + osc_m^2(f)).



The proof utilizes a quasiinterpolation operator.

 
For any ∈ there exists an interpolation operator Π_h: H^2_0(Ω)→ V_C such that,
 for 0≤ k≤ m≤ 2 and  φ∈ H^2_0(Ω),  

    φ-Π_hφ_H^k(K)≲
    h_K^m-k|φ|_H^q(ω_K)

holds for any in the triangle K∈ 
and  the  interior ω_K of the union  ω_K    of the triangles in  
sharing a vertex with K. 



This follows from <cit.> once the required scaling properties of the degrees
of freedom are clarified. The Argyris or the HCT finite element schemes involve some
normal derivative and  do not form an affine finite element family, but 
an almost affine finite element element family <cit.>. It is by now understood that
this guarantees the appropriate scaling properties. This is  explicitly  calculated in  <cit.> for the
HCT finite elements and also follows  for the Argyris finite elements, as employed 
e.g. in <cit.>. Since the result is frequently accepted
<cit.>,  further details are omitted. 


Proof of Theorem <ref>.
Continue the notation of the proof of Theorem <ref>  with
X=Y=V, X_h=Y_h=V_C, Q=1,  etc. and recall that, for sufficiently small δ,
Theorem <ref>  guarantees u-u_C < β/Γ.
Hence Corollary <ref> implies (for v_h≡ u_C) 

    u-u_C≤  C_ rel,1N(u_C)_V^*.

With Π_h from Lemma <ref>,  some 
appropriate  ϕ∈ V with ϕ=1  satisfies 

    N(u_C)_V^*= N(u_C;ϕ)= N(u_C;ϕ-Π_hϕ).
 
Two successive integrations by parts result in

    a(u_C,ϕ-Π_hϕ)= (Δ^2 u_C) (ϕ-Π_hϕ)  +[Δ u_C]_E ∇(ϕ-Π_hϕ)·ν_E 
                       - 
     (ϕ-Π_hϕ)[(D^2 u_C)]_E ·ν_E.

An integration by parts in the nonlinear term 
Γ(u_C,u_C,ϕ-Π_hϕ) leads to

    Γ(u_C,u_C,ϕ-Π_hϕ)=∑_K∈∫_K Δ u_C∇ u_C· Curl(ϕ-Π_hϕ)
       =∑_K∈∫_K (ϕ-Π_hϕ) curl(-Δ u_C∇ u_C)
    +∑_E∈∫_E(ϕ-Π_hϕ) [Δ u_C∇ u_C]_E·τ_E.

Those identities show that  (<ref>) is equal to a sum 
over edges of jump contributions plus a sum
over triangle of volume contributions; the latter is  

    ( Δ^2 u_C- curl (Δ u_C∇ u_C)-f   ) (ϕ-Π_hϕ)  ≲∑_K∈η_K h_K^-2 ||   ϕ- Π_h ϕ_L^2(K)

and controlled  with standard manipulations 
based on Lemma <ref> (with k=0 and m=2) and the finite
overlap of the patches (ω_K:K∈). The jump contributions include 
some trace inequality as well and are otherwise standard as in linear problems
that involve the bi-Laplacian. For instance, the nonlinear jump contribution 
for each edge E reads

    ∫_E(ϕ-Π_hϕ) [Δ u_C∇ u_C]_E·τ_E
    =∫_E(ϕ-Π_hϕ) [Δ u_C]_E  ∇ u_C·τ_E

in case of  an interior edge E shared by the two triangles T_+ and T_- that form the patch
ω_E and  vanishes in case of a boundary edge E⊂∂Ω
(with ϕ=Π_hϕ=0 on ∂Ω). 
The continuity of ∇ u_C leads to the previous equality. This term is
controlled by the residual  
h_E^3/2  [Δ u_C]_E  ∇ u_C·τ_E_L^2(E) times 

    h_E^-3/2ϕ-Π_hϕ_L^2(E)≲   h_E^-2ϕ-Π_hϕ_L^2(T_±) 
    +h_E^-1ϕ-Π_hϕ_H^1(T_±)≲ | ϕ |_H^2(ω_T_±)

with a trace inequality on one of the two triangles T_± in the first and 
 Lemma <ref> (for k=0,1)  in the second estimate. 
 The remaining terms are
controlled in a similar way.


Some words are in order about the term 
h_E^3/2  [Δ u_C]_E  ∇ u_C·τ_E_L^2(E), in which  an inverse inequality along the interior edge E=∂ T_+∩∂ T_- 
(shared by T_±∈) of the polynomial
∇ u_C·τ_E (unique as a trace from T_±) shows 
∇ u_C·τ_E _L^∞(E)≲ h_E^-1 u_C _L^∞(E).
This and the global continuous embedding  
H^2(Ω)↪ L^∞(Ω)
leads to

    h_E^3/2  [Δ u_C]_E  ∇ u_C·τ_E_L^2(E)≲ h_E^1/2  [Δ u_C]_E _L^2(E) u_C .

Since u_C ≲ 1, the nonlinear edge contribution 
 is controlled  by another contribution
h_E^1/2  [Δ u_C]_E _L^2(E) to η_E;
in other words,
this nonlinear edge contribution can be omitted. 

The overall strategy in the efficiency proof 
follows the  bubble-function technique due to Verfürth 
 <cit.>. The emphasis in this paper is on the nonlinear contributions 
and on the interaction of the various nonlinear terms with the volume estimator. 
We will give two examples only to illustrate some details and start with 
the cubic bubble-function  b_K ∈ W^1,∞_0(K) 
(the product of all three barycentric coordinates times 27)
of the triangle K∈ with 0≤ b_K≤max b_K=1. 
Let f_K:=Π_m f∈ P_m(K) be the
 L^2(K)  orthogonal polynomial projection of f∈ L^2(K) for  degree 
 m∈ℕ_0 
so that f-f_K_L^2(K)=h_K^-2 osc_m(f,K).

Since g:=  Δ^2 u_C- curl (Δ u_C∇ u_C)-f_K is a polynomial
of degree at most max{k-4,(k-2)(k-1)-1,m} (recall that k is the 
degree of the  finite element functions), an inverse estimate 
reads    g_K^2≲∫_Kρ_K g  for the test function
ρ_K:=b_K^2g∈ H^2_0(K)⊂ V. 
The above integrations by parts (<ref>)-(<ref>)  with the test function 
ϕ-Π_hϕ replaced by ρ_K are restricted to K for 
the support of b_K and ∇ b_K  is K. 
This leads to the first equality in 

    ∫_K g ρ_K   =a(u_C, ρ_K)+  Γ(u_C,u_C,ρ_K) -∫_Kρ_Kf_K
       =a(u_C-u, ρ_K)+  Γ(u_C,u_C,ρ_K) - Γ(u,u,ρ_K)
    +∫_Kρ_K(f-f_K)

and  (<ref>)  leads to the second.  
Except for the last term (that leads to oscillations in the end), 
elementary algebra,  Γ(u,u,ρ_K)- Γ(u_C,u_C,ρ_K) 
=Γ(u-u_C,u,ρ_K)+Γ(u_C,u-u_C,ρ_K),
Cauchy, and Hölder inequalities bound the 
above terms upto a constant  by

    u-u_C_H^2(K)( (1 +  | u |_W^1,∞(Ω)) ρ_K_H^2(K) 
    +| u_C|_H^2(Ω)   | ρ_K|_W^1,∞(K)).

The  inverse estimates  
ρ_K_H^2(K) + | ρ_K|_W^1,∞(K)≲ h_K^-2ρ_K_L^2(K)≤ h_K^-2 g _L^2(K)
lead in the preceding estimates (after division by  h_K^-2 g _L^2(K))
to

    h_K^2    g_L^2(K)≲ u-u_C_H^2(K)+ osc_m(f,K).

This and a  triangle inequality  prove efficiency 
η_K≲ u-u_C_H^2(K)+ osc_m (f,K) of the volume contribution.

The patch  ω_E of an  interior edge  E∈ is the interior of the union of the two neighbouring 
triangles in  sharing the edge E and may be a non-convex quadrilateral.  Observe that the
shape-regularity in  implies the  shape-regularity of 
the largest rhombus  R   contained in the patch ω_E that has E as one diagonal. 
Let b_R∈ H^1_0(R)⊂ H^1_0(ω_E) be the (piecewise quadratic) 
edge-bubble function of E in R (with 0≤ b_R≤max b_R=1) 
and let Φ_E∈ P_1(R) be the affine function that vanishes along E
and satisfies ∇Φ_E =  h_E^-1ν_E.
Then b_E:= Φ_E b_R^3∈ H^2_0(R)⊂ H^2_0(ω_E)  satisfies 
∇ b_E ·ν_E= h_E^-1 b_R^3 along E and
| b_E|_L^∞(ω_E)≲ 1 as in  <cit.>. 
Extend [Δ u_C]_E constantly in the normal direction to E and set 
ϱ_E:=h_E^2 [Δ u_C]_E  b_E ∈ H^2_0(R)⊂ H^2_0(ω_E).  An inverse 
estimate in the beginning,
∇ϱ_E·ν_E=h_E b_R^3[Δ u_C]_E on E, 
and piecewise integrations by parts lead to 

    h_E [Δ u_C]_E_L^2(E)^2  ≲  h_Eb_R^3/2[Δ u_C]_E_L^2(E)^2
    =∫_E ∇ϱ_E·ν_E [Δ u_C]_E
            =∫_ω_E(Δ u_CΔϱ_E -ϱ_EΔ^2_pw u_C).

The test-function  ϱ_E  in (<ref>) shows  that,
the right-hand side reads

    a(u_C-u,ϱ_E)+Γ(u_C,u_C,ϱ_E)-Γ(u,u,ϱ_E)
     + ∫_ω_E(f-Δ_pw ^2 u_C+ curl_pw (Δ u_C∇ u_C))ϱ_E.

A Cauchy inequality in the first, the arguments for (<ref>) in the second term, and the  bound  
(η_T_++η_T_-)h_E^-2ϱ_E _L^2(ω_E) for the third term lead to 

    h_E [Δ u_C]_E_L^2(E)^2  ≲(  u-u_C_H^2(ω_E)  +η_T_++η_T_-)  
     (h_E^-2ϱ_E _L^2(ω_E)+ |ϱ_E |_H^2(ω_E)).

The function ϱ_E is polynomial in each of the two open triangles in R∖ (E∪∂ R)
and allows for inverse estimates. Since  |b_E|≲ 1 a.e., this proves that the last factor is
controlled by

    h_E^-2ϱ_E_L^2(ω_E)≲[Δ u_C]_E_L^2(ω_E)≲  h_E^1/2[Δ u_C]_E_L^2(E)

for the constant extension of  [Δ u_C]_E in the direction of ν_E in the last step.
The combination of the previous two displayed inequalities with the above efficiency of the 
volume contribution concludes the proof of 

    h_E^1/2[Δ u_C]_E_L^2(E)≲u-u_C_H^2(ω_E)+η_T_++η_T_-≲u-u_C_H^2(ω_E)+ osc_m (f,{T_+,T_-}).

The efficiency of h_E^3[ (D^2   u_C)]_E ·ν_E^2_L^2(E) is also established through an adoption of the corresponding arguments in  <cit.>. Hence the
straightforward details are omitted. 




 §.§ Morley FEM

The nonconforming Morley element space V_h:=() 
associated with the  triangulation  of the polygonal domain 
Ω⊂ℝ^2 into triangles reads

    ():={ v_M∈ P_2()        v_M  is continuous at (Ω)  and vanishes at (∂Ω), 
    
       ∫_E[∂ v_M/∂ν]_E=0
      for all   E ∈ (Ω), 
    
       ∫_E∂ v_M/∂ν=0
      for all E∈ (∂Ω)
    }.

The discrete formulation 
seeks u_M∈() such that 

    N_h(u_M;v_M):=a_pw( u_M, v_M)-F( v_M)
    +Γ_pw( u_M, u_M, v_M)=0  v_M∈().

Here and throughout this section, V:=V+ () is endowed with the mesh-dependent 
norm 
φ_pw
:=√(a_pw(φ,φ))
for φ∈V and, for all η,χ,ϕ∈(), 

    a_pw(η,χ)
       :=∑_K∈∫_K D^2 η:D^2χ,
    Γ_pw(η,χ,ϕ)
       := ∑_T∈∫_T Δη(∂χ/∂ x_2∂ϕ/∂ x_1-∂χ/∂ x_1∂ϕ/∂ x_2).


The a priori error estimate  means  best-approximation upto first-order terms and so refines 
<cit.> for the Morley FEM and generalises it for any regular solution.  


If u∈ H^2_0(Ω) is a regular solution to N(u)=0, then there exist positive 
ϵ, δ, and ρ such that  (A)-(C) hold for all ∈(δ) with

    apx() ≲u-I_ M u_ pw
    +   h_Δ u ∇ u 
    +  osc_0(f,)≲ h_ max^s.



The proof requires the following four lemmas. 

 
For any v∈ V+ (), the Morley interpolation 
I_M(v)∈()  defined by

    (I_M v)(z)=v(z)  for any  z∈(Ω)  and ∫_E∂ I_M v/∂ν_E=∫_E∂ v/∂ν_E for any  E∈


satisfies (a) D^2_ pw I_M =Π_0 D^2 and (b)


    h_K^-2(1-I_M)v _L^2(K)+h_K^-1∇(1-I_M) v_L^2(K)
    +D^2 I_Mv_L^2(K)≲D^2v_L^2(K).



Let H^2((ω_K)) denote the piecewise H^2 functions on 
the neighbourhood  ω_K, piecewise with respect to  the triangulation 
(ω_K) of all triangles T with zero distance to K∈. Let
|∙|_H^2((ω_K)) be the corresponding seminorm
as the local contributions of ∙_pw associated with ω_K.

 
 
There exists an enrichment operator E_M:()→ V such that φ_M∈() satisfies

    (a)  ∑_m=0^2 h_K^2m|φ_M-E_Mφ_M|_H^m(K)^2 ≲  h_K^4|φ_M|_H^2((ω_K))^2 K∈; 
       (b)   h_^-2(φ_M-E_Mφ_M)_^2≲∑_E∈ h_E [D^2φ_M]_Eτ_E_L^2(E)^2
       ≲φ_M-E_Mφ_M_ pw^2≤Λmin_φ∈ V D_h^2(φ_M-φ)_^2; 
       (c)    I_ME_Mφ_M=φ_M,  and  φ_M-E_Mφ_M⊥ P_0() in L^2(Ω).



The Sobolev embeddings for conforming functions depend on the domain 
Ω, while their discrete counterparts for nonconforming functions require
particular attention.



For any 1≤ p<∞, there exists a constant C=C(Ω,p,∢) 
(which depends on p, Ω, and the shape regularity of ) with

    v_L^∞(Ω)+v_W^1,p()≤ C 
    v_ pw  for all v∈ H^2_0(Ω)+  ().
 



The main observation is that the enrichment operator E_M from 
Lemma <ref> maps into the HCT finite element space  plus squared 
bubble-functions  <cit.>; 
so  v_M-E_Mv_M is a piecewise polynomial of degree at most 6 for
any v_M∈()
(with respect to some refinement of , where each triangle T is divided into three 
sub-triangles by connecting each vertex with its center of inertia). This leads to 
inverse estimates such as 

    | v_M-E_Mv_M |_W^1,∞(T)≲ h_T^-1 | v_M-E_Mv_M |_H^1(T)≲ h_T^-2 v_M-E_Mv_M _L^2(T).

Lemma <ref>.b shows for v∈ H^1_0(Ω) that the right-hand side is controlled by
 
    h_^-2 (v_M-E_Mv_M) _L^2(Ω)≲ v_M-E_Mv_M _pw≤Λmin_φ∈ V v_M-φ_pw≤Λ v+v_M _pw.

Since T∈ is arbitrary, this proves

    | v_M-E_Mv_M |_W^1,∞(Ω,)≲ v+v_M _pw.

Since  v_M-E_Mv_M is Lipschitz continuous with  Lipschitz constant 
| v_M-E_Mv_M |_W^1,∞() and
 vanishes at the vertices of T∈,  

    || v_M-E_Mv_M ||_L^∞(Ω)≤ h_max | v_M-E_Mv_M |_W^1,∞(Ω,)
 
holds for the maximal mesh-size 
h_max≤diam(Ω).
This and (<ref>)    imply  (with C_1≈ 1)

    v_M-E_Mv_M _L^∞(Ω)≤  C_1  v+v_M_pw.
 
The boundedness of  the  continuous 2D Sobolev embedding 
H^2(Ω)↪  L^∞ (Ω) leads to  ∙_L^∞(Ω)≤ C_2 ∙ in H^2_0(Ω).
Consequently, with  a triangle inequality in the beginning, 

    v+v_M _L^∞(Ω)   ≤v_M-E_M v_M _L^∞(Ω)+
      v+ E_Mv_M _L^∞(Ω)
       ≤  C_1 v+v_M_pw+  C_2   v+ E_M v_M .

The  triangle inequality and Lemma <ref>.b (again with φ=-v) show

    v+E_M v_M≤ v+v_M _pw+  v_M -E_Mv_M _pw≲ v+v_M _pw.

The combination of (<ref>) with the previously displayed estimate shows 
the first assertion 
v+v_M _L^∞(Ω)≲ v+v_M _pw.
The proof of the second assertion is similar with 
(<ref>)-(<ref>).  
The boundedness of  the  continuous 2D Sobolev embedding 
H^2(Ω)↪ W^1,p(Ω) leads to 
| ∙  |_W^1,p(Ω)≤ C(p,Ω) ∙ in H^2_0(Ω).
Consequently,

    | v+v_M |_W^1,p(Ω,)   ≤  | v+E_M v_M  |_W^1,p(Ω)
    +| v_M -E_Mv_M|_W^1,p(Ω,)
       ≤ C(p,Ω)    v+E_M v_M + |Ω|^1/p
    | v_M -E_Mv_M|_W^1,∞(Ω,)

with the area |Ω|≈ 1≈ C(p,Ω). 
Recall  (<ref>)  and (<ref>) in the end to control the 
previous upper bound in terms of
v+E_M v_M +   v+ v_M _pw≲  v+ v_M _pw. This concludes the proof of the second assertion
 | v+v_M |_W^1,p(Ω,)≲  v+ v_M _pw.



The bound for  a_pw is immediate from 
(<ref>) for the norm  ∙_pw
in V≡ V+ ().  The bound   
Γ_pw =√(2) C(Ω,4,∢)^2
in

    |  Γ_pw(η,χ,ϕ)|≤Γ_pwη_pwχ_pwϕ_pw  for all η, χ, ϕ∈V≡ V+ ()

follows  from (<ref>) 
with Hölder inequalities and Lemma <ref>.


 For 1/2<s≤ 1 there 
exists a positive constant C such that any η∈ H^2+s(Ω)
and  φ_M∈() satisfy 
a_pw(η,φ_M-E_Mφ_M)≤ Ch_max^sη_H^2+s(Ω)φ_M_ pw. 

Proof of Theorem <ref>. 
Set X=Y=V, X_h=Y_h=V_h, X=V+V_h, a(∙,∙):=a_pw(∙, ∙), b(∙,∙):=Γ_pw( u,∙,∙)+Γ_pw(∙,u,∙) and P=I_M, Q=𝒞= E_M. The regularity u∈ H^2+s(Ω) of Subsection <ref> with s>1/2 allows for  the 
bounded global Sobolev embeddings   
H^2+s(Ω) ↪ W^2,4(Ω)↪ W^1,∞(Ω).
This and   Lemma <ref>
lead for θ∈X and ϕ∈ H^1(Ω) to

    |Γ_pw(u,θ,ϕ)|+|Γ_pw(θ,u,ϕ)|
       ≲(  u_W^2,4(Ω)θ_W^1,4(Ω,)
    +  u_W^1,∞(Ω)θ_pw)ϕ_H^1(Ω)
       ≲u_H^2+s(Ω) θ_pwϕ_H^1(Ω).

For θ_M∈() with θ_M_pw=1, the aforementioned estimates imply that b(θ_M,∙)∈ H^-1(Ω) 
and so the solution z∈ V to the biharmonic problem

    a(z,ϕ)=b(θ_M,ϕ)ϕ∈ V

satisfies z∈ H^2+s(Ω) and z _H^2+s(Ω)≲ 1
<cit.>.
The regularity  z∈ H^2+s(Ω) and Lemma <ref> 
(resp. Lemma <ref>) 
imply  (H1)  (resp.  (H2)) 
with δ_1≲ h_ max^s (resp. δ_2≲ h_ max^s).
The estimate (<ref>) and Lemma <ref> verify (H3) with 
δ_3≲ h_ max.  Lemma <ref>.b leads to (H4)
with Λ_4=Λ. For any y_M∈ M() with 
y_M_pw=1, 
Lemma <ref> guarantees 

    a_pw(u, y_M-E_My_M)≲ h_max^s  u_H^2+s(Ω)≈ h_max^s,

while  Lemma <ref> shows 

    F(y_M-E_My_M)≲ h_max^2 f≲  h_max^s.
 
This implies  (H5) with δ_5≲  h_max^s. 
Choose x_h=I_M u so that (H6) holds with δ_6≲ h_ max^s.

In conclusion,  for sufficiently small mesh-size h_max, 
the discrete inf-sup inequality of Theorem <ref> holds  with β_h≥β_0>0. 
Moreover, Theorems <ref> and <ref> apply and prove
(A)-(C). To compute apx()=N(u)_M()^*,  let 
ϕ_M ∈ M() satisfy  ϕ_M_pw=1
and apx() = N(u; ϕ_M). Since N(u, E_Mϕ_M)=0, 

the difference 
ψ
:=ϕ_M-E_M ϕ_M∈V satisfies 

    apx()     = N(u;ψ ) 
    = a_ pw (u- I_M u, ψ)-F((1-Π_0)ψ)+
    Γ_ pw(u,u,ψ)

with  Lemma <ref>.c for a_ pw (I_M u, ψ)=0 and
Π_0 ψ=0 a.e. in the last step.
This, the finite overlap of (ω_K:K∈) in Lemma <ref>.a and <ref>.b 
 for ψ_pw≲ 1 
lead to (<ref>).
 



 §.§ A posteriori error estimate

For any K∈ and E∈, define the volume and edge error estimators  by 

    η_K^2:=h_K^4 curl(-Δ u_M∇ u_M)-f_L^2(K)^2 and 
       η_E^2 :=h_E[D^2 u_M]_Eτ_E_L^2(E)^2
    +h_E^3[Δ u_M∇ u_M]_E·τ_E_L^2(E)^2
    +h_E^3{Δ u_M∇ u_M}_E·τ_E_L^2(E)^2.

Here and throughout this section, the average of ϕ∈X across the interior edge E=∂ K_+∩∂ K_-∈(Ω) shared by two triangles K_±∈  reads
{ϕ}_E:=(ϕ|_K_++ϕ|_K_-)/2,
while {ϕ}_E:=ϕ|_E along any boundary edge E∈(∂Ω).


If u∈ V is a regular solution to (<ref>), then there exist positive δ,ϵ, and   C_ rel 
such that, for any ∈(δ), the discrete solution u_M∈() to (<ref>) 
with u-u_M_ pw≤ϵ  satisfies

    C_ rel^-2 u-u_M_ pw^2   ≤∑_K∈η_K^2+∑_E∈η_E^2.



 

Let u_M be the solution to (<ref>) close to u and  
apply Theorem <ref> with  
X=Y=V, X_h=Y_h=V_h,  v_h= u_M, and Q:=E_M from Lemma <ref>. 
Suppose that  ϵ,δ satisfy Theorem <ref> and,
if necessary, are chosen smaller
such that, for any ∈(δ), exactly one discrete solution
u_ M∈ X_M  to (<ref>) satisfies  u-u_ M_ pw≤ϵ≤β/( 2(1+Λ )Γ
 ).  Lemma <ref>.b implies
u_M-E_M u_M  _ pw≤Λ   u -u_M _ pw≤Λϵ.
This and  triangle inequalities show

    E_M u _M +  u_M_pw   ≤ u_M-E_M u_M _pw + 2   u_M_pw≤  2  u  +(2+ Λ)ϵ=:M;
    
     u- E_Mu_M   ≤  u- u_M_pw+   u_M- E_M u_M_pw≤ (1+Λ)ϵ≤β/( 2Γ
    ).

Consequently,  the abstract residual (<ref>) in Theorem <ref> implies

    u- u_M_pw≤ 2β^-1N(E_M u_M)_ V^*
    +  u_M- E_M u_M_pw.

There exists some  ϕ∈ V with ϕ=1 and

    N(E_Mu_M)_V^*= N(E_M u_M;ϕ)= a(E_M u_M,ϕ)-F(ϕ)
    +Γ(E_M u_M,E_M u_M,ϕ)
       =N( u_M;ϕ)+ a_pw(E_M u_M- u_M,ϕ)+Γ(E_M u_M,E_M u_M,ϕ)
    -Γ_pw( u_M, u_M,ϕ)

with the definition of N and of   N. This,
the  bound  of  a_pw,    elementary arguments  with the trilinear form
and its bound Γ_ pw from Remark <ref>,   
and M prove

    N(E_M u_M)_ V^*≤N( u_M;ϕ)+(1+ M Γ_ pw)
     u_M- E_M  u_M_ pw.
 

Since  u_M  solves  (<ref>),  
N_h( u_M;ϕ)=N_h( u_M;χ) holds for χ:=ϕ-I_M ϕ with 
the  Morley interpolation I_M ϕ of ϕ. 

Since Lemma <ref>.a  implies a_pw( u_M,ϕ-I_Mϕ)=0, 
an integration by parts in the nonlinear term Γ_pw(∙,∙,∙) leads to

    N( u_M;ϕ)   =∑_K∈∫_K Δ u_M∇ u_M· curlχ
    -F(χ)
       =∑_K∈∫_K ( curl(-Δ u_M∇ u_M)-f)χ
    +∑_E∈∫_E [Δ u_M∇ u_M·τ_E]_E{χ}_E
         +∑_E∈∫_E{Δ u_M∇ u_M}_E·τ_E
      [χ]_E.

This and standard arguments with Cauchy and trace inequalities
plus  Lemma <ref>.b with ϕ=1 eventually lead to
some constant C_A≈ 1 with

    C_A ^-2N( u_M;ϕ)^2 ≤∑_K∈η_K^2+∑_E∈η_E^2.


Piecewise inverse estimates  
u_M- E_M  u_M_pw≲  h_^-2 ( u_M- E_M  u_M )_L^2(Ω)

and  Lemma <ref>.b with the tangential jump residuals 
lead to some constant C_B≈ 1 with 

    C_B ^-2 u_M- E_M  u_M_pw^2 ≤∑_E∈ h_E [D^2 u_M]_Eτ_E_L^2(E)^2.

This is bounded by ∑_E∈η_E^2.

The combination of  (<ref>)-(<ref>)
concludes the proof with 
 C_ rel=   2 β^-1 C_A+  (1+ 2β^-1(1+M Γ _pw)) C_B.



The efficiency of the estimator remains  as an open question  owing to the average term 
{Δ u_M∇ u_M·τ_E}_E_L^2(E) in η_E 
(for the remaining contributions
are efficient). 
The sum of all those contributions associated to those terms,  however, converge
 (at least) with linear rate in  that
 
    S:= (∑_E∈h_E^3{Δ u_M∇ u_M·τ_E}_E_L^2(E)^2
     )^1/2≲  h_max u _H^2+s(Ω) u_M_pw
     =O(h_max).

 Before a sketch of the proof concludes this remark, it should be stressed that 
 (<ref>)  can be 
 a higher-order term:  Consider a uniform mesh in a singular situation with 
 re-entering corners (with  an exact solution of reduced regularity  
 u∉ H^3(Ω) <cit.>) with a  suboptimal convergence rate
s<1. Then   S in  (<ref>) is of higher-order. 

The proof of (<ref>) starts with a  triangle inequality 

    S^2 
     ≤∑_T∈∑_E∈(T)
     h_E^3 (Δ u_M∇ u_M) |_T_L^2(E)^2.

The discrete trace inequality  (i.e. a trace inequality followed by an inverse inequality) 
for each summand shows
 
    h_E^3 (Δ u_M∇ u_M) |_T_L^2(E)^2
     ≲ h_E^2Δ u_M∇ u_M_L^2(T)^2.
 
Recall  the piecewise constant mesh-size  h_∈ P_0(),
h_|_T:=h_T:=diam(T) for T∈, with maximum
h_max:=max h_≤δ.
The shape regularity of  shows

    S≲  h_Δ_pw u_M∇_pw u_M_L^2(Ω)≤  h_Δ u ∇_pw u_M_L^2(Ω) +
       h_Δ_pw (u-u_M)∇_pw u_M_L^2(Ω)

with a triangle inequality in the last step.
Recall that u∈ H^2+s(Ω) for s>1/2 enables 
the bounded embedding H^s(Ω)↪ L^2p(Ω)
for any p with 1<p<1/(1-s).
This and a Hölder inequality with 1/p+1/p'=1 leads to

    Δ u ∇_pw u_M_L^2(Ω)≤Δ u _L^2p(Ω)∇_pw u_M_L^2p'(Ω).

Lemma <ref> shows that the last term is controlled 
by  u_M_pw. Consequently,

    h_Δ u ∇_pw u_M_L^2(Ω)≲ h_max u _H^2+s(Ω) u_M_pw.

The analysis of the second term  starts with 0<r≤ s≤ 1 and 
the elementary observation 

    h_Δ_pw (u-u_M)∇_pw u_M_L^2(Ω)≤√(2) h_max^1-r u-u_M_pw 
    | h_^r    u_M |_W^1,∞(Ω,).

The  asserted convergence rate follows  with
 u-u_M_pw≲  h_max^s  u _H^2+s(Ω). The maximum of the remaining term 
| h_^r    u_M |_W^1,∞(Ω,)=h_T^r |   u_M |_W^1,∞(T)
is attained for (at least) one T∈. An inverse inequality and 
Lemma <ref>  in the end show

    h_T^r |   u_M |_W^1,∞(T)≲ | u_M |_W^1,2/r(T)≤ 
    | u_M |_W^1,2/r(Ω,)≲ u_M_pw.

Consequently, 
h_Δ_pw (u-u_M)∇_pw u_M_L^2(Ω)≲  h_max^1+s-r u_M_pw.
The combination of the previous estimates proves  (<ref>). 



The lack of local efficiency  is part of a  more general structural 
difficulty. Whenever volume terms require a piecewise integration by parts with Morley  
finite element test functions, there arise average 
terms like {ϕ}_E in Theorem <ref>, which are not residuals.
This prevents an efficiency analysis in this section as well as in <cit.> or 
<cit.>. It is left as an open problem for future research and may cause 
a modification of the discrete scheme. In the vibration of a biharmonic plate or in the von Kármán
equations of the subsequent section, this difficulty does not arise.




§ VON KÁRMÁN EQUATIONS

Given a load function f∈ L^2(Ω), the von Kármán equations model the deflection of a very thin elastic plate with vertical displacement u∈ and the Airy stress function v∈ such that

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

With  the co-factor matrix  (D^2 v )  of D^2 v, the von Kármán brackets read

    [u,v]:=∂^2 u/∂ x_1^2∂^2 v/∂ x_2^2
    +∂^2 u/∂ x_2^2∂^2 v /∂ x_1^2
    -2∂^2 u/∂ x_1∂ x_2∂^2 v/∂ x_1∂ x_2=(D^2 u):D^2 v.




 §.§ Continuous problem

The  weak formulation  of  the von Kármán equations (<ref>)
seeks  u,v∈ V:=H^2_0(Ω) with


    a(u,φ_1)+ γ(u,v,φ_1)+γ(v,u,φ_1)   =
    (f,φ_1)_L^2(Ω)φ_1∈ V
    
    a(v,φ_2)-γ(u,u,φ_2)      =0            φ_2 ∈ V.

 
Here and throughout this section  abbreviate, for all η,χ,φ∈ V, 

    a(η,χ):= D^2 η:D^2χ  and  γ(η,χ,φ):=- [η,χ]φ.

The abstract theory of Sections <ref>-<ref> applies for the real
Hilbert space  X:=V× V  with its dual X^* to  the operator N:  X→  X^*  defined by

    N(Ψ;Φ):=⟨  N(Ψ), Φ⟩:= 
    A(Ψ,Φ)-F(Φ)+Γ(Ψ,Ψ,Φ)

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

    A(Θ,Φ)   :=a(θ_1,φ_1)+a(θ_2,φ_2),
    
    F(Φ)   :=(f,φ_1)_,
    Γ(Ξ,Θ,Φ)   :=γ(ξ_1,θ_2,φ_1)+γ(ξ_2,θ_1,φ_1)-γ(ξ_1,θ_1,φ_2).

Note that   A(∙,∙) is a scalar product in X  and the trilinear 
form Γ(∙,∙,∙) is bounded <cit.>.

It is known <cit.> that there exist
a solution Ψ∈ X with N(Ψ)=0. Any solution has the regularity 
Ψ∈ H^2+s(Ω):=(H^2+α(Ω))^2 for  1/2<s≤ 1 depending on the polygonal
bounded Lipschitz domain Ω <cit.>. This allows for the
boundedness 

    Γ(Ψ,Θ,Φ)≤ C Ψ_H^2+s(Ω)ΘΦ_H^1(Ω)  for any Θ∈ X and Φ∈ H^1_0(Ω;^2).




 §.§ Conforming FEM

With the notation of Section <ref> on V_C⊂ H^2_0(Ω), the conforming 
finite element formulation seeks Ψ_C=(u_C,v_C)∈ X_h:=V_C× V_C such that

    N(Ψ_C;Φ_C)=0   for all  Φ_C ∈  X_h.



If Ψ∈ X is a regular solution to N(Ψ)=0, then  
there exist positive  ϵ, δ, and ρ such that  (A)-(C) hold with
apx()≡ 0  for all ∈(δ). 



The proof  is analogous to that of Theorem <ref> 
and hence omitted. The a priori error analysis is derived
in  <cit.> with a fixed point iteration (of linear convergence).


For any K∈ and E∈, define  the volume and edge error
estimators by 

    η_K^2   :=h_K^4Δ^2 u_C-[u_C,v_C]-f^2_L^2(K)
    +h_K^4 Δ^2 v_C+1/2[u_C,u_C]^2_L^2(K),
    η_E^2   :=h_E^3[ (D^2   u_C)]_E·ν_E^2_L^2(E)
    +h_E^3[ (D^2   v_C)]_E·ν_E^2_L^2(E)
             +h_E[D^2 u_Cν_E]_E·ν_E^2_L^2(E)
    +h_E[D^2 v_Cν_E]_E·ν_E^2_L^2(E) .



If Ψ∈ X  is a regular solution  to N(Ψ)=0,  then there exist positive 
δ, ϵ, C_ rel, and C_ eff 
such that, for all ∈(δ), the unique discrete solution  
Ψ_C=(u_C,v_C)∈  X_h to (<ref>) with Ψ-Ψ_C   <ϵ  satisfies 

    C_ rel^-2Ψ-Ψ_C^2   ≤∑_K∈η_K^2
    +∑_E∈η_E^2≤ C_ eff^2(Ψ-Ψ_C^2+ osc_0^2(f)).


 

For Y=X,  Y_h=X_h, we proceed as in the proof of Theorem <ref>
and, for sufficiently small δ, derive  (H1)- (H6) and  
u-u_C < β/Γ from 
Theorem <ref>.  
Hence Corollary <ref> implies for v_h≡Ψ_C= (u_C, v_C) that

    Ψ-Ψ_C≲N(Ψ_C)_X^*=N(Ψ_C;Φ)

for some  Φ∈ X with Φ=1 and its approximation 
Π_hΦ∈ X_h
(Π_h from Lemma <ref> applies componentwise). 
Abbreviate (χ_1,χ_2):=χ :=Φ-Π_hΦ and deduce from  
 (<ref>) that 
N(Ψ_C)_X^*=N(Ψ_C;(χ_1,χ_2)).
Successive integrations by parts show

    A(Ψ_C,χ )= (Δ^2 u_C) χ_1   + (Δ^2 v_C) χ_2  
          +[D^2 u_C]_Eν_E·∇χ_1
    + [D^2 v_C]_Eν_E·∇χ_2 
         - χ_1 [(D^2 u_C)]_E·ν_E
    - χ_2 [(D^2 v_C)]_E·ν_E .

This and the definition of Γ(∙,∙,∙) lead to the residual

    A(Ψ_C,Φ-Π_hΦ)-F(Φ-Π_hΦ)+Γ(Ψ_C,Ψ_C,Φ-Π_hΦ)
         = (Δ^2 u_C-[u_C,v_C]-f)χ_1+(Δ^2 v_C+[u_C,u_C])χ_2
           -([(D^2 u_C)]_E·ν_E)χ_1
    +[D^2 u_C]_Eν_E·∇χ_1
            -([(D^2 v_C)]_E·ν_E)χ_2
    + [D^2 v_C]_Eν_E·∇χ_2.

The two edge terms in the above expression that involve ∇χ_j for j=1,2  can be rewritten as    

    [D^2 u_C]_Eν_E·∇χ_1+[D^2 v_C]_Eν_E·∇χ_2
         =[D^2 u_Cν_E]_E·ν_E ∂χ_1/∂ν+[D^2 v_Cν_E]_E·ν_E ∂χ_2/∂ν
            +[D^2 u_Cν_E]_E·τ_E ∂χ_1/∂τ+[D^2 v_Cν_E]_E·τ_E ∂χ_2/∂τ.

The last two terms involve tangential derivatives and so vanish for u_C and v_C belong to .

Standard arguments analogous to <cit.> with a Cauchy  inequality, an inverse inequality, and Lemma <ref>  conclude the proof of the reliability. 

The proof of the efficiency of the volume term η_K is immediately adopted from  that of 
<cit.>. The arguments in the proof of efficiency for the edge terms 
h_E[D^2 u_Cν_E]_E·ν_E_L^2(E)  
and h_E[D^2 v_Cν_E]_E·ν_E_L^2(E) are the same 
as for the (linear) biharmonic  equation and can be adopted from 
<cit.> or <cit.>. 
Further  details are omitted. 




 §.§ Morley FEM

The Morley FEM
seeks Ψ_M∈  X_M:=() ×()⊂X:=X+X_M 
(endowed with the norm ∙_pw) such that

    N_h(Ψ_M;Φ_M):=A_pw(Ψ_M,Φ_M)
    +Γ_pw(Ψ_M,Ψ_M,Φ_M)-F(Φ_M)=0 Φ_M ∈  X_M.

Here and throughout this subsection,  for all Ξ=(ξ_1,ξ_2), 
Θ=(θ_1,θ_2), Φ=(φ_1,φ_2)∈X,

    A_pw(Θ,Φ):=a_pw(θ_1,φ_1)+a_pw(θ_2,φ_2), 
       F(Φ):= fφ_1, 
       Γ_pw(Ξ,Θ,Φ):=b_pw(ξ_1,θ_2,φ_1)+b_pw(ξ_2,θ_1,φ_1)-b_pw(ξ_1,θ_1,φ_2),
 

and, for all η,χ,  φ∈V:=H^2_0(Ω)+ (),

    a_pw(η,χ):= D^2 η:D^2χ and  
     b_pw(η,χ,φ):=- [η,χ]φ.
 
(The boundedness of  a_pw is immediate and  that of Γ_pw 

follows from Lemma <ref>.) 


If Ψ∈ X is a regular solution to N(Ψ)=0, then  
there exist positive  ϵ, δ, and ρ such that  (A)-(C) hold
for any ∈(δ) with 

    apx()≲Ψ-I_MΨ_ pw +  osc_0(f+ [u,v],)
    +  osc_0([u,u],)≲  h_ max^s.



 
Set Y=X, Y_h=X_M, X= X+X_M, a(∙,∙):=A_pw(∙,∙), b(∙,∙):=2Γ_pw( Ψ,∙,∙) and P=I_M,  Q=𝒞=E_M. 
Given  Ψ∈ H^2+s(Ω),  
Θ, Φ∈X, piecewise inequalities and 
the bounded global  Sobolev imbedding  H^2+s(Ω)↪ W^2,4(Ω) (for s>1/2)
show

    Γ_pw(Ψ,θ,Φ)
    
    
    ≲Ψ_H^2+s(Ω)Θ_pwΦ_L^4(Ω).

For Θ_M ∈ X_M with Θ_M_pw=1, 
the linear functional Γ(Ψ,Θ_M,∙)∈ H^-1(Ω)
leads to a  unique solution Z∈  X to the biharmonic problem 
A( Z,Φ)=Γ(Ψ,Θ_M,Φ) for all Φ∈ X 
with Z∈ H^2+s(Ω) <cit.>. For  φ_M∈(), the inverse estimate

    φ_M-E_Mφ_M_L^4(K)≤ C h_K^-1/2φ_M-E_Mφ_M_L^2(K) for all   K∈,

the bound for Γ_pw, and 
Lemma <ref>.a imply δ_3≲ h_ max^3/2.  The 
remaining conditions for the parameters in the (H1)-(H2) and (H4)-(H6) 
are verified as in the proof of Theorem <ref>.
For some Φ_M∈ X_M with Φ_M_pw=1,  
apx()=N(Ψ)_X_h^*=N(Ψ; Φ_M).
This, N(Ψ; E_MΦ_M)=0,  (<ref>),  and 
Lemmas <ref>-<ref> 
lead for (χ_1,χ_2):=χ:=Φ_M-E_MΦ_M to

    apx()     = N(Ψ; Φ_M-E_M Φ_M) 
    = A_pw(Ψ, χ )-F(χ)
    +Γ_pw(Ψ,Ψ,χ) 
        = A_pw(Ψ-I_M Ψ,χ)
    - (f+ [u,v], χ_1)_L^2(Ω)+1/2  ([u,u], χ_2)_L^2(Ω)
       ≲Ψ-I_M Ψ_ pw + osc_0(f+ [u,v],)
    + osc_0([u,u],)≲  h_ max^s

with arguments as in  the final part of the proof  of Theorem <ref>.
Hence Theorems <ref> and <ref> 
apply and prove (A)-(C).


For any K∈ and E∈, define the volume and edge 
error estimators  by 

    η_K^2   := h_K^4[u_M,v_M]+f_L^2(K)^2+
    h_K^4 [u_M,u_M]_L^2(K)^2,
    η_E^2    := h_E[D^2 u_M]_Eτ_E_L^2(E)^2
    +h_E[D^2 v_M]_Eτ_E_L^2(E)^2.
 


If Ψ=(u,v)∈  X is a regular solution to N(Ψ)=0, then there exist 
δ,ϵ, C_ rel, and C_ eff such that,
for any ∈(δ), the discrete solution
Ψ_M=(u_M,v_M)∈ X_M:=()×() to (<ref>)
with Ψ-Ψ_ M_ pw≤ϵ  satisfies

    C_ rel^-2Ψ-Ψ_ M_ pw^2
    	   ≤∑_K∈η_K^2+∑_E∈η_E^2≤ C_ eff^2 
    	(Ψ-Ψ_M_ pw^2+ osc_0^2(f)).





Let Ψ_M be the solution to (<ref>) close to Ψ and  
apply Theorem <ref> with  Y=X, Y_h=X_M,  v_h=Ψ_M, and Q=E_M.
Suppose that  ϵ,δ satisfy Theorem <ref> and,
if necessary, are chosen smaller such that, for any ∈(δ), exactly one discrete solution
Ψ_M∈ X_M  to (<ref>) satisfies 
 Ψ-Ψ_M_ pw≤ϵ≤β/( 2(1+Λ )Γ
 ). Lemma <ref>.b implies
 Ψ_ M- E_MΨ_ M_ pw≤ΛΨ-Ψ_M_ pw≤Λϵ.
This and  triangle inequalities show

    E_MΨ_M + Ψ_M_pw   ≤Ψ_M-E_MΨ_M _pw + 2  Ψ_M_pw≤  2 Ψ +(2+ Λ)ϵ=:M;   
    Ψ- E_MΨ_M   ≤Ψ- Ψ_M_pw+ Ψ_M- E_MΨ_M_pw≤ (1+Λ)ϵ≤β/( 2Γ
     ).

Consequently,  the abstract residual (<ref>) in Theorem <ref> implies

    Ψ-Ψ_M_pw≤  2β^-1N(E_MΨ_M)_ X^*+Ψ_M-E_MΨ_M_pw.

There exists Φ∈  X with Φ=1 and 

    N(E_MΨ_M)_ X^* =  N(E_MΨ_M;Φ)=A_ pw(E_MΨ_M,Φ)-F(Φ)+Γ(E_MΨ_M,E_MΨ_M,Φ)

with the definition of N. This and the definition of N(Ψ_M;Φ) lead to

    N(E_MΨ_M)_ X^*   =N(Ψ_M;Φ)+A_pw(E_MΨ_M-Ψ_M,Φ)
    
         +Γ(E_MΨ_M,E_MΨ_M,Φ)-Γ_pw(Ψ_M,Ψ_M,Φ)
       ≤N(Ψ_M;Φ) +
    (1+MΓ_ pw )Ψ_M-E_MΨ_M_pw

with the  bound  of  A_pw, elementary arguments  with the trilinear form and 
its bound  Γ_ pw 
(deduced from Lemma <ref>
as in Remark <ref>), 
and M  in the second step.
Since Ψ_M solves  (<ref>), 
N(Ψ_M;Φ)=N(Ψ_M; χ)
holds for χ:=(χ_1,χ_2):=Φ-I_MΦ  with the  Morley interpolation I_M Φ of Φ.  
Since Lemma <ref>.a  implies A_pw(Ψ_M,χ)=0, the  
definitions of Γ_ pw(∙,∙,∙)  and F(∙) lead to

    N(Ψ_M; Φ)   =
    Γ_pw(Ψ_M,Ψ_M,χ)  - F(χ) 
       = 1/2( [u_M,u_M], χ_2 )_L^2(Ω)   - (f+[u_M,v_M], χ_1)_L^2(Ω)
       ≤ ( ∑_K∈η_K^2 )^1/2 h_^-2(Φ-I_MΦ)_L^2(Ω)≤ C_A ( ∑_K∈η_K^2 )^1/2

with weighted Cauchy inequalities in the second last step and the constant C_A≈ 1 from Lemma <ref>.b  with Φ=1  in the end. 

The combination with (<ref>) reads

    Ψ-Ψ_M_pw≤  2 β^-1  C_A( ∑_K∈η_K^2 )^1/2
    +(1+  2β^-1 (1+MΓ _pw )    )  
    Ψ_M-E_MΨ_M_pw.

The last term  is  controlled as in (<ref>) and  this concludes the proof of the 
reliability estimate with 
C_ rel= max{ 2β^-1 C_A ,(1+ 2β^-1(1+M Γ _pw))C_B}.

The proof of the efficiency of the volume term η_K 
 is immediately adopted from  that of  <cit.>.
The arguments in the proof of efficiency for the edge term η_E are the same 
as for the (linear) biharmonic  equation and can be adopted from  <cit.>.
Further details are omitted. 



The adaptation of the nonconforming scheme that allows for a right-hand side f ∈ H^-1(Ω) is possible by arguments in Section <ref>.




§ ACKNOWLEDGEMENTS


The authors thank for the comments by the anonymous referee that led to Subsection <ref>.
The research of the first author has been supported by the Deutsche Forschungsgemeinschaft in the Priority Program 1748 under the project "foundation and application of generalized mixed FEM towards nonlinear problems
in solid mechanics" (CA 151/22-2).  The research of the second author is supported by the NBHM Grant 0204/58/2018/R&D-II/14721. The finalization of this paper has been supported by   DST SERB MATRICS grant of the second author MTR/2017/000199 and SPARC project 
(id 235) entitled the mathematics and computation of plates.





§ APPENDIX



















































 §.§ Proof of Lemma <ref>

The first part of the assertion is included in <cit.> 
and so is merely outlined for convenient reading of the second. 
The Rellich compact embedding theorem 
H^1_0(Ω)c↪ L^2(Ω) 
leads to c↪ H^-1(Ω) in the sequel. 
Hence S:={ g∈  | g=1} is pre-compact in V^*. The operator 
A∈ L(V;V^*), associated to the scalar product a via Av=a(v,∙) for all 
v∈ V (note A in contrast to the coefficients A); A is 
invertible and A^-1∈ L(V^*;V) maps S onto W:=A^-1(S) pre-compact in H^1_0(Ω). The open balls B(z,ϵ/6) in V 
around  z∈W  with radius ϵ/6 with respect to the norm ∙_a
form an open cover of the compact set W and so have a finite sub-cover for z_1,…, z_J∈W, 

    W⊂∪_j=1^J  B(z_j,ϵ/6)⊂ V.

Since (Ω) is dense in V, there exists ζ_j∈(Ω) with 
z_j-ζ_j_a <ϵ/6. The smoothness of ζ_j proves 

    ζ_j-I_Cζ_j_a  ≤ C h_max≤ C δ

for any triangulation ∈(δ) and the nodal interpolation I_C in S_0^1(); the constant C depends on max_j=1, …,  JD^2ζ_j, 
the shape-regularity parameter κ, and on λ. 
For any g∈∖{0} with 
z=A^-1(g)/g∈ W from (<ref>), there exists at least one index j∈{1,…,J} with 
z∈ B(z_j,ϵ/6). This, the choice of ζ_j, and   (<ref>) 
with δ:=ϵ/(6C)   prove

    z-I_Cζ_j_a ≤ z-z_j_a+
     z_j-ζ_j_a+ζ_j-I_Cζ_j_a<ϵ/3+Cδ <ϵ/2.

A rescaling of this leads to 
A^-1(g)-g I_Cζ_j_a =g z-I_Cζ_j_a≤ϵg/2.
This proves that the first term in the asserted inequality  is bounded by the right-hand side. 
The analysis of the second term considers the pre-compact subset 
A∇ W={ A∇ z: Tz=g∈ S} of L^2(Ω;^n), where T: L^2(Ω) ⟶ H^1_0(Ω) is the solution map with z=Tg ∈ H^1_0(Ω).  Since 
the open balls B(Q,ϵ/6)  around Q∈ A∇ W 
in the L^2 norm form an open cover
of the compact closure A∇ W in L^2(Ω;^n), 
there exists Q_1,…, Q_K in  A∇ W with

    A∇ W ⊂∪_k=1^K  B(Q_k,ϵ/6)⊂ L^2(Ω;^n).

Since (Ω;^n) is dense in L^2(Ω;^n), 
there exists Φ_k∈(Ω;^n) with 
Q_k-Φ_k <ϵ/6. The smoothness of Φ_k and a Poincaré inequality 
(on simplices with constant h_T/π) prove

    Φ_k-Π_0 Φ_k≤ ||∇Φ_k||   h_max/π≤ C δ

for any triangulation ∈(δ) with the L^2 projection Π_0 onto 
P_0(;^n). The constant C=max{ ||∇Φ_1||,…, ||∇Φ_K||} depends on
the smoothness of the  functions Φ_1 ,…, Φ_K. 

For any g∈∖{0} with 
z=A^-1(g)/g∈ W from (<ref>), there exists at least one index k∈{1,…,K} with 
A∇ z∈ B(Q_k,ϵ/6). This, the choice of Φ_k, and   (<ref>) 
with δ:=ϵ/(6C)   prove

    (1-Π_0)  A∇ z  ≤ A∇ z - Π_0 Φ_k ≤ A∇ z -Q_k+Q_k -Φ_k+(1-Π_0)Φ_k<ϵ/2.

A rescaling of this proves
(1-Π_0)  A∇ z  ≤ϵg/2 for all Az=g∈ L^2(Ω)
(with arbitrary norm g≥ 0).
This concludes the proof. 

amsplain


