A saturation property for the spectral-Galerkin approximation of a Dirichlet problem in a square
Claudio Canuto Ricardo H. Nochetto Rob Stevenson Marco Verani
math.NA


Scalable Prototype Selection by Genetic Algorithms and Hashing t1
    [
    December 30, 2023
=================================================================





Both practice and analysis of adaptive p-FEMs and hp-FEMs raise
the question what increment in the current polynomial degree p
guarantees a p-independent reduction of the Galerkin error. We
answer this question for the p-FEM in the simplified context of homogeneous
Dirichlet problems for the Poisson equation in the
two dimensional unit square with polynomial data of
degree p. We show that an increment proportional to p yields a
p-robust error reduction and provide computational evidence that a
constant increment does not.  






§ MOTIVATION AND STATEMENT OF THE RESULT



High order finite element methods (FEMs) exhibit exponential
convergence for elliptic problems with piecewise analytic data, and
thus have become the methods of choice in computational science and
engineering for such problems. The seminal work of Babuška and
collaborators <cit.>
has established the mathematical foundations for
the a priori design of meshes and distribution of polynomial degrees,
and proved exponential convergence for corner and edge singularities.
In contrast, adaptive hp-FEMs hinge on a posteriori error
estimators, which help determine whether it is more convenient to
locally refine the mesh or increase the polynomial degree to improve
the resolution. Although exponential convergence is observe
experimentally, it has never been proved rigorously with the exception of
<cit.>.

Our adaptive hp-FEM of <cit.> hinges on a coarsening module due to
Binev <cit.>, which in turn guarantees instance optimality and
thus exponential convergence. As any other adaptive hp-FEM, ours
also has a module to reduce the PDE error by a fixed fraction for
piecewise polynomial data thereby avoiding data oscillation. Such
module in <cit.> relies on the a posteriori error estimator of
Melenk and Wohlmuth <cit.> for dimension d=2 and does not possess
optimal complexity. In <cit.> we turn to the equilibrated flux
residual estimator of Braess, Pillwein and Schöberl <cit.>, and
Ern and Vorahlík <cit.>,
and show that the issue of optimal complexity
reduces to studying three model problems with polynomial data in the
reference triangle for d=2. We present overwhelming computational
evidence in <cit.> supporting the fact that to reduce the Galerkin
error by a fixed factor, the polynomial degree p must be increased by
an amount proportional to p.

In this paper we take over this question again in a further simplified
setting and give a rigorous answer. We consider the Poisson equation


    -Δ u = f   in Ω,
          
      u = 0   on ∂Ω,


over the unit square Ω = (-1,1)^2 of ℝ^2 with
polynomial f of degree p.

Let us first introduce some notation. Let I = (-1,1) be the reference
element so that Ω = I^2. For p≥ 0, let ℙ_p(Ω) denote the space of polynomials of total degree ≤
p restricted to Ω,
and let _p := ℙ_p(Ω) ∩ H^1_0(Ω).
Since the latter space reduces to {0} for p<4, we will consider it only for p≥ 4.

We equip H^1_0(Ω) with the energy inner product (·,·)_H^1_0(Ω):=(∇·,∇·)_L^2(Ω)^2 and resulting norm ·_H^1_0(Ω). Analogous definitions apply to H^1_0(I).


Given any f ∈ℙ_p(Ω), let u=u(f) ∈ H^1_0(Ω) be the variational solution of (<ref>), i.e.,

    (u,v)_H^1_0(Ω) = (f,v)_L^2(Ω)    ∀ v ∈ H^1_0(Ω).

For any q ≥ 4, let u_q =u_q(f) ∈_q be the corresponding
Galerkin projection of u onto _q, i.e., 

    (u_q,v_q)_H^1_0(Ω) = (f,v_q)_L^2(Ω)    ∀ v ∈_q.


We are interested in finding sufficient conditions on q=q(p)>p so that


    u-u_q_H^1_0(Ω)≤αu-u_p_H^1_0(Ω)


for 0<α<1 independent of p. If (<ref>) holds,
we say that the error reduction is p-robust. In view of
Phytagoras equality


    u-u_p_H^1_0(Ω)^2 = u-u_q_H^1_0(Ω)^2
      + u_q-u_p_H^1_0(Ω)^2,


which is a consequence of Galerkin orthogonality u_q-u_p ⊥ u-u_q
in H^1_0(Ω), we see that (<ref>) is equivalent
to the following saturation property


    √(1-α^2) u-u_p_H^1_0(Ω)≤u_q-u_p_H^1_0(Ω).











We next
observe that (<ref>) is equivalent to the simpler
saturation property


    √(1-α^2) u_H^1_0(Ω)≤u_q_H^1_0(Ω).


To see this just define v := u - u_p ∈ H^1_0(Ω), which is the
solution of (<ref>) with polynomial forcing
f-Δ u_p ∈ℙ_p(Ω) and Galerkin solutions v_p=0 and v_q = u_q-u_p.
We aim at establishing the following rigorous result.


There exists a constant C>0 such that for all λ>1,
any mapping p↦ q = q(p) satisfying q(p) > max(λ p,p+4)
yields


    ‖ u ‖_H^1_0(Ω)≤ 
    C λ/λ-1 ‖ u_q ‖_H^1_0(Ω)    for all   p≥ 0  and  
     f ∈ℙ_p(Ω).



Since most hp-FEMs in the literature perform p-enrichment upon
adding a constant increment to p, typically 1 or 2, one may
wonder whether the preceding sufficient condition on q is also
necessary. We now investigate this question computationally upon
defining


    C_p,q,r :=max_f∈ℙ_p(Ω)u_r_H^1_0(Ω)/u_q_H^1_0(Ω)


where r≫ q is chosen computationally so that u_r is sufficiently close to u in
H^1_0(Ω); note that this is not a hidden saturation assumption
because the value of r is not predetermined but found once
the number C_p,q,r stabilizes. This calculation reduces to an eigenvalue problem,
already used in <cit.>, and leads to Figure
<ref> for q=p+k with k=2,4,6,10:



We thus realize that C_p,q,r exhibits a modest but linear growth
on q=p+k for k constant, which confirms that this choice is not
p-robust. For moderate values of p this might still be acceptable
computationally, but it could compromise computational complexity for
extreme values of p as in spectral algorithms <cit.>.

Even though the saturation property is quite delicate, it has been
often used in a posteriori error analysis of low order AFEMs until
now. It originates in the work of Bank and Weiser <cit.>, and
Bornemann, Erdmann and Kornhuber <cit.>; see Nochetto
<cit.> for related work. Dörfler and Nochetto <cit.> proved
the saturation property for p=1, q=2 and d=2 provided data
oscillation is small relative to u-u_p_H^1_0(Ω) but
showed counterexamples for piecewise constant forcing f.

We stress that (<ref>) is not asymptotic: it is valid for
any p≥0 and any f∈ℙ_p(Ω). Since u_q→ u in
H^1_0(Ω) as q→∞ it is obvious that C_p,q,∞→
1 as q→∞. It is for this reason that Theorem
<ref> has some intrinsic value in the theory of FEMs and
might have implications beyond a posteriori error analysis.

The proof of Theorem <ref> proceeds as follows. We perform a
multilevel decomposition of _q


    _q = ⊕_j=1^q _j,


where _j are polynomial subspaces of total degree j. Since this
decomposition is quasi-orthogonal in the sense that


    _j ⊥_ℓ  for all ℓ j-2, j, j+2,


we need to account for interactions between neighboring spaces
_j. We study the angle between subspaces _j and show it is
larger than π/3; this is the content of Proposition
<ref>. This in turn allows us to find the precise decay of high
frequency modes of u_q, which leads to (<ref>).

The paper is organized as follows. In section <ref> we
introduce the multilevel decomposition (<ref>) and
discuss a few properties including Proposition <ref>.
In section <ref> we analyze
the decay of high order components of u_q, whereas in section
<ref> we prove Theorem <ref>.
We conclude in section <ref> with the proof of
the rather technical Proposition <ref>.









§ MULTI-LEVEL DECOMPOSITIONS OF POLYNOMIAL SPACES



Hereafter, we recall the definition of classical polynomial bases in L^2(Ω) and in H^1_0(Ω), obtained by tensorization from
corresponding bases in L^2(I) and in H^1_0(I), where I=(-1,1) is the reference interval. The elements of theses bases enjoy certain orthogonality properties, by which a multi-level, quasi-orthogonal decomposition of H^1_0(Ω) is obtained. This will be useful in deriving the main result of this paper.

On the interval I, we consider the orthonormal Legendre basis in L^2(I)

    ϑ_k(x) =√(k+12)  L_k(x),      k ≥ 0,

(where L_k stands for the k-th Legendre orthogonal polynomial in I, 
which satisfies  deg  L_k = k and L_k(1)=1), as well as the
orthonormal Babuška-Shen (BS) basis in H^1_0(I):

    φ_k(x)   =√(k-12)∫_x^1 L_k-1(s) ds 
       =1/√(4k-2)(L_k-2(x)-L_k(x)), 
          k ≥ 2 .

The BS basis enjoys the following orthogonality properties in L^2(I)
for m≥ k:

    (φ_k,φ_m)_L^2(I) = 
    2/(2k-3)(2k+1)   if  m=k  , 
    
    - 1/(2k+1)√((2k-1)(2k+3))   if  m=k+2  , 
    
    0    otherwise.


On the square Ω=I × I, the previous bases induce, resp., the tensorized orthonormal Legendre basis in L^2(Ω):

    Θ_k(x)= ϑ_k_1(x_1)ϑ_k_2(x_2),      k ∈K̂,

where k=(k_1,k_2), x=(x_1, x_2) and  K̂=ℕ^2, and the tensorized Babuška-Shen basis in  H^1_0(Ω):

    Φ_k(x) = φ_k_1(x_1) φ_k_2(x_2),       k ∈ K ,

where K={k ∈ℕ^2  :  k_i ≥ 2  for i=1,2 }.  

The tensorized BS basis is not orthogonal in H^1_0(Ω). Indeed, from the expression


    (Φ_k,Φ_m)_H^1_0(Ω)   =
    (φ_k_1,φ_m_1)_H^1_0(I)(φ_k_2,φ_m_2)_L^2(I)
       +
     (φ_k_1,φ_m_1)_L^2(I)(φ_k_2,φ_m_2)_H^1_0(I) ,


and (<ref>) we immediately obtain 

    (Φ_k,Φ_m)_H^1_0(Ω)≠ 0  iff   
    k_1=m_1  and  k_2-m_2 ∈{-2,0,2},  or 
    
    k_2=m_2  and  k_1-m_1 ∈{-2,0,2}.

As a consequence, denoting by |k|=|k_1|+|k_2| the ℓ^1-norm in ℤ^2,  we have

    (Φ_k,Φ_m)_H^1_0(Ω) = 0     if    |k|-|m| ∉{-2,0,2}.


At last, concerning the interaction between  the Legendre basis and the BS one, we have

    (Θ_k,Φ_m)_L^2(Ω)≠ 0  iff   k_1 ∈{m_1-2,m_1} and  k_2 ∈{m_2-2,m_2},

which implies

    (Θ_k,Φ_m)_L^2(Ω) = 0     if    |k-m| > 4.



Any function v ∈ L^2(Ω) can be split uniquely into four components

    v = ∑_α∈{0,1}^2 v^α,

where v^α for α=(α_1, α_2) is even (odd, resp.) with respect to the variable x_i (i=1,2) iff α_i=0 (α_i=1, resp.). 
For instance, v^(0,1) satisfies v^(0,1)(-x_1,x_2)=v^(0,1)(x_1,x_2) and v^(0,1)(x_1,-x_2)=-v^(0,1)(x_1,x_2) for all
(x_1,x_2)∈Ω.

Components with different parity indices are always L^2(Ω)-orthogonal, and H^1_0(Ω)-orthogonal whenever v ∈ H^1_0(Ω). In particular, as a consequence of (<ref>) and (<ref>), we observe that (Φ_k,Φ_m)_H^1_0(Ω) = 0 and (Θ_k,Φ_m)_L^2(Ω) = 0 whenever k and m have at least one entry of different parity.








 §.§ Detail spaces and their projectors



For j ≥ 4, let us define the finite dimensional subspace of H^1_0(Ω)

    _j := span{Φ_k   :  |k|=j } .

Note that, thanks to (<ref>), the functions Φ_k that generate _j are mutually orthogonal in H^1_0(Ω). 
We immediately have the multi-level decompositions

    _q = ⊕_j=4^q _j    for all  q ≥ 4      and     H^1_0(Ω) = ⊕_j=4^∞_j  .

Such decompositions are `quasi-orthogonal', in the sense that by (<ref>) we have

    _j ⊥_H^1_0(Ω)_ℓ    for all  ℓ≠ j-2, j, j+2.

Furthermore, the `angle' between two non-orthogonal subspaces is uniformly bounded away from 0, as implied by the
following technical result, that will be crucial in the sequel.
We postpone its proof to section <ref>.

 
Let P_j : _j-2→_j (j ≥ 6) be the orthogonal projection with respect to the H^1_0(Ω)-inner product. 
Then,

    ‖ P_j ‖_ L(_j-2,_j) < 1/2.

Actually, there exists a constant c>0 independent of j such that 

    ‖ P_j ‖_ L(_j-2,_j)≤1/2( 1 - c/  j^2)  .






Note that the orthogonal projection P_j^*: _j→_j-2, given by the adjoint of P_j, satisfies the same estimate.







§ DECAY OF THE HIGHER-ORDER COMPONENTS OF THE GALERKIN SOLUTION



Given f ∈ℙ_p(Ω), let  u_q ∈_q be the
Galerkin solution defined in (<ref>), and let u_q =
∑_j=4^q U_j, with U_j=U_j(q) ∈_j, be its multilevel decomposition according to (<ref>). The purpose of this section is to prove that for any q sufficiently larger than p, the H^1_0(Ω)-norm of  U_q and U_q-1 decay at least proportionally to the quantity (q-p)^-1. The precise result is as follows.



For any p ≥ 0 and q > p̂ := p+4, one has 

    ‖ U_j ‖_H^1_0(Ω)≤6/q-p ‖ u_q ‖_H^1_0(Ω) ,       j=q,   q-1.



We first observe that the parity splitting (<ref>) of the forcing f induces by linearity a corresponding splitting of the Galerkin solution u_q as well as of each of its multi-level details U_j, which is nothing but the parity splitting of u_q as well as of U_j. Therefore, thanks  to the orthogonality of the components with different parity (cf. Remark <ref>), it is enough to establish (<ref>) for each component separately, and then sum-up the squares of both sides invoking Parseval's identity.

For the sake of definiteness, we will focus on the components of (even, even) type, the other types being amenable to a similar treatment. Thus, referring to (<ref>) for the notation, we consider the component u_q^(0,0) of u_q (which solves (<ref>) for the forcing f^(0,0)), as well as its details U_j^(0,0)∈_j^(0,0):= span{Φ_k   :  |k|=j  and  k_1, k_2  are even}. We aim at proving that for q ≥p̂+1 and j ∈{q-1,q}

    ‖ U_j^(0,0)‖_H^1_0(Ω)≤6/q-p ‖ u_q^(0,0)‖_H^1_0(Ω).

However, it is easily seen that U_q-1^(0,0)=0 if q is even, and similarly  U_q^(0,0)=0 if q is odd. Hence, we will prove

    ‖ U_q^(0,0)‖_H^1_0(Ω)≤6/q-p ‖ u_q^(0,0)‖_H^1_0(Ω)

under the assumption that q is even, the other situation being similar.

To avoid cumbersome notation, for the rest of the proof we will drop the superscript u^(0,0) from all entities. 
So, we will write


    u_q='∑_j=4^q U_j
        with  U_j ∈_j,


where here and in the sequel the symbol ^' indicates that the summation runs over even indices only.

From the Galerkin equations, we have for any even j ∈ [4,q]

    (u_q,W_j)_H^1_0(Ω)=(f,W_j)_L^2(Ω)    for all
       W_j∈_j   .


Since q ≥p̂+1, exploiting (<ref>) and (<ref>),
(<ref>) yields


    (U_q,W_q)_H^1_0(Ω)+(U_q-2,W_q)_H^1_0(Ω)=0    for all  W_q∈_q  ,

which can be rewritten equivalently as 

    U_q=-T_q^-1P_q U_q-2

where T_q=I and P_q is defined in Proposition <ref>.

For any even j satisfying p̂+2<j ≤ q-2, (<ref>) yields

    (U_j+2,W_j)_H^1_0(Ω)+(U_j,W_j)_H^1_0(Ω)+(U_j-2,W_j)_H^1_0(Ω)=0


for all W_j∈_j;
this is equivalent to  P_j+2^* U_j+2 + U_j + P_j U_j-2 = 0.
Assuming by induction that U_j+2=-T_j+2^-1 P_j+2 U_j with
T_j+2^-1≤ 2 (which is trivially true for j=q-2
according to (<ref>)), gives 

    (I-P_j+2^* T_j+2^-1 P_j+2) U_j=-P_j U_j-2.
 
Using Proposition <ref>, the operator

    T_j:=I-P_j+2^* T_j+2^-1 P_j+2

is invertible and satisfies

    T^-1_j≤1/1- 1/4 T^-1_j+2 .

We conclude that

    U_j=-T_j^-1P_j U_j-2.

with T^-1_j≤ 2, which proves the induction argument for all even j satisfying p̂ + 2 < j ≤ q-2.

Next, we have to bound the norm of U_j-2 for j=p̂+4 when p, hence p̂, is even, or for j=p̂+3 when p is odd. It is therefore convenient to define the even integer

    r:=p̂ = p+4    if  p  is even, 
    p̂-1 = p+3    if  p  is odd,

so that in both cases, we have to bound U_r+2_H^1_0(Ω).
To this end, let us introduce


    _r^(0,0):= '⊕_j=4^r _j,
        u̅_r:='∑_j=4^r
      U_j∈_r^(0,0);


note that u̅_r  u_r because U_j=U_j(q).
Then, in view of (<ref>), we deduce


    (u̅_r+'∑_j=r + 2^q U_j, v_r)_H^1_0(Ω)=(f,v_r)_L^2(Ω)  for
        all  v_r∈_r^(0,0) ,

which, thanks to ('∑_j=r + 2^q U_j, v_r)_H^1_0(Ω)=(P_r^*U_r+2,v_r)_H^1_0(Ω) for all v_r ∈_r^(0,0), implies that

    u̅_r+P_r^* U_r+2=u_r.








We observe that (<ref>) is also valid for j=r+2.
Since (U_r,W_r+2)_H^1_0(Ω)= (u̅_r,W_r+2)_H^1_0(Ω)
we obtain


    (U_r +4,W_r +2)_H^1_0(Ω)+(U_r +2,W_r +2)_H^1_0(Ω)+(u̅_r,W_r +2)_H^1_0(Ω)= 0

and

    T_r +2 U_r +2 = -P_r +2u̅_r

as it happened with (<ref>). This implies 

    (T_r +2-P_r +2P_r^*) U_r +2=-P_r +2 u_r

which in view of (<ref>) yields

    (I-P_r +4^* T_r +4^-1 P_r +4 -P_r +2P_r^*) U_r +2=-P_r +2 u_r.

Since P_r +4^* T_r +4^-1 P_r +4 +P_r +2P_r^*≤ 2
  1/4+1/4=3/4, thanks to Proposition <ref> and T_r +4^-1≤2, we conclude that

    U_r +2_H^1_0(Ω)≤1/1-3/41/2u_r_H^1_0(Ω)=2u_r_H^1_0(Ω)≤ 2 ‖ u_q‖_H^1_0(Ω),

where last inequality follows from the inclusion
_r^0,0)⊂_q^(0,0) and the
minimization property of the Galerkin solution.


























Collecting the above results we arrive at

    U_j = -T_j^-1P_j U_j-2     r +4 ≤   j ≤ q    (j even), 
        U_r +2_H^1_0(Ω)≤ 2 ‖ u_q ‖_H^1_0(Ω).

For q ≥p̂+4, this implies 

    U_q_H^1_0(Ω)≤'∏_j=r +4^q T_j^-1P_j  U_r +2_H^1_0(Ω)≤   2  ‖ u_q ‖_H^1_0(Ω)'∏_j=r+4^q T_j^-1P_j .


In order to bound the product on the right-hand side, let us write j=q-2m with
m=0,1, …, s and s:= 1/2(q-r)-2.
Then, by Proposition <ref>, we have T_j^-1P_j ≤1/2 T_j^-1 =: α_m. Recalling (<ref>), it holds

    α_m ≤1/2/1-1/2α_m-1= 1/2-α_m-1,     with  α_0 ≤12  .

By recurrence, it is immediate to check that α_m ≤m+1/m+2, whence

    '∏_j=r+4^q T_j^-1P_j ≤∏_m=0^s α_m ≤∏_m=0^s m+1/m+2 =1/s+2=
    2/q-r.

Since q≥ p+6 if p is even and q≥ p+5 if p is odd,
it is easily checked that

    2/q-r =2/(q-p)-4≤6/q-p   if  p  is even, 
    2/(q-p)-3≤5/q-p   if  p  is odd.

This gives the desired estimate (<ref>).







§ A SUBSPACE DECOMPOSITION IN H^1_0(Ω)  



Consider the complementary space of _q in H^1_0(Ω) given by

    _q^c :=_H^1_0(Ω)span {Φ_m   :   |m| > q }.


Therefore, H^1_0(Ω)=_q ⊕_q^c and 
any v ∈ H^1_0(Ω) can be split as


    v = v_q + z_q,     
    v_q ∈_q,    z_q ∈_q^c.


The purpose of this section is to apply once more Proposition <ref> and derive a bound on the norm of v_q and z_q in terms of the norm of v. 

We start with the following auxiliary result for any w=∑_j=4^q W_j∈_q.


For any q ≥ 4 and any w=∑_j=4^q W_j∈_q, one has

    ‖ W_q ‖_H^1_0(Ω)≤√(2)‖ w ‖_H^1_0(Ω) .
 

As in the previous section, splitting w and W_q in their orthogonal components according to the parity of the basis functions, it is enough to establish the result for each component separately. Hereafter, we detail the analysis for the `(even, even)' case, in which case we may assume q even, since otherwise W_q^(0,0)=0 and the result is trivial. 

Dropping as above the superscript u^(0,0) in functions
and subspaces, we write w=W+W_q with


    W = '∑_j=4^q-2 W_j ∈_q-2.


Keeping W_q fixed, let us first minimize the norm of w, i.e., let
us look for the minimizer W̅∈_q-2 of the quantity
Ψ(W):=‖ W+W_q‖_H^1_0(Ω)^2. Such a function
satisfies


    (W̅, Y)_H^1_0(Ω) = - (W_q, Y)_H^1_0(Ω)    for all  Y ∈_q-2

and

    Ψ(W̅) = ‖ W_q ‖_H^1_0(Ω)^2 +(W_q, W̅)_H^1_0(Ω) .

Using the orthogonality conditions (<ref>), we obtain the sequence of equations

    (W̅_4,Y_4)_H^1_0(Ω)+(W̅_6,Y_4)_H^1_0(Ω)=0    for all  Y_4∈_4 ,

and

    (W̅_j-2,Y_j)_H^1_0(Ω)+(W̅_j,Y_j)_H^1_0(Ω)+(W̅_j+2,Y_j)_H^1_0(Ω)=0    for all  Y_j∈_j

for any even j such that 4 < j <q-2, and finally

    (W̅_q-4,Y_q-2)_H^1_0(Ω)+(W̅_q-2,Y_q-2)_H^1_0(Ω)=-( W_q,Y_q-2)_H^1_0(Ω)

for all Y_q-2∈_q-2.
Setting recursively T_4=I and T_j=(I-P_j T_j-2^-1 P_j^*),  we derive W̅_j =-T_j^-1 P_j+2^* W̅_j+2 for j=4,6, …, q-4, and W̅_q-2=-T_q-2^-1 P_q^* W_q. Note that, thanks to Proposition <ref>, one can prove as in Sect. <ref> that ‖ T_j^-1‖≤ 2 for all j. Since

    (W_q, W̅)_H^1_0(Ω)   = (W_q, W̅_q-2)_H^1_0(Ω)
       =  (P_q^* W_q, W̅_q-2)_H^1_0(Ω)
      = - (P_q^* W_q, T_q-2^-1 P_q^* W_q)_H^1_0(Ω),

using once more Proposition <ref>, we deduce

    (P_q^* W_q, T_q-2^-1 P_q^* W_q)_H^1_0(Ω)   ≤‖ T_q-2^-1‖ ‖ P_q^* ‖^2 ‖ W_q ‖_H^1_0(Ω)^2 
       ≤ 2 (1/2)^2 ‖ W_q ‖_H^1_0(Ω)^2 =1/2‖ W_q‖_H^1_0(Ω)^2 .

In view of (<ref>) and the preceding estimate,
we conclude that 

    ‖ w ‖_H^1_0(Ω)^2 = Ψ(W) ≥Ψ(W̅) ≥‖ W_q ‖_H^1_0(Ω)^2 -1/2‖ W_q ‖_H^1_0(Ω)^2 = 1/2‖ W_q ‖_H^1_0(Ω)^2

for any w ∈_q^(0,0), whence the asserted estimate follows.


We now establish the main result of this section.


There exists a constant C_2>0 such that for any q ≥ 4 and any
v=v_q + z_q ∈_q ⊕_q^c, one has 

    ‖ z_q ‖_H^1_0(Ω)≤ C_2    q  ‖ v ‖_H^1_0(Ω).



Using ‖ z_q ‖_H^1_0(Ω)≤‖ v ‖_H^1_0(Ω) + ‖ v_q ‖_H^1_0(Ω), it is enough to prove the existence of a constant C_2'>0 independent of q such that for all v ∈ H^1_0(Ω) 

    ‖ v_q ‖_H^1_0(Ω)≤ C_2' q ‖ v ‖_H^1_0(Ω) .

To this end, let us focus as above on the `(even, even)' components of
v and v_q, in which case it is not restrictive to assume q even,
and drop the superscript u^(0,0) in functions and subspaces.
Let us fix any even integer r>q and assume first that v ∈_r
is written as v=v_q+V, with


    V= '∑_j=q+2^r V_j ∈'⊕_j=q+2^r _j.


By applying the same technique as above, i.e., minimizing the (squared) norm Ψ(V):=‖ v_q+V‖_H^1_0(Ω)^2, we find that

    ‖ v ‖_H^1_0(Ω)^2 = Ψ(V) ≥Ψ(V̅) = ‖ v_q ‖_H^1_0(Ω)^2 +(v_q, V̅)_H^1_0(Ω) ,

where the minimizer V̅ =  '∑_j=q+2^r V̅_j is such that V̅_q+2=-T_q+2^-1P̃_q+2 v_q, for 
T_q+2 defined recursively by (<ref>) with T_r=I, and P̃_q+2 : H^1_0(Ω) →_q+2 defined as the orthogonal projection in the H^1_0(Ω) inner product. Now,

    (v_q, V̅)_H^1_0(Ω)   = (v_q, V̅_q+2)_H^1_0(Ω) = (P̃_q+2v_q, V̅_q+2)_H^1_0(Ω)
       = - (P̃_q+2v_q, T_q+2^-1P̃_q+2 v_q)_H^1_0(Ω)

with

    | (P̃_q+2v_q, T_q+2^-1P̃_q+2 v_q)_H^1_0(Ω)| ≤‖ T_q+2^-1‖ ‖P̃_q+2 v_q ‖_H^1_0(Ω)^2 ≤ 2 ‖P̃_q+2 v_q ‖_H^1_0(Ω)^2 .
 
Writing v_q ='∑_j=4^q V_j, one has P̃_q+2 v_q = P̃_q+2 V_q =  P_q+2 V_q, whence by
Proposition <ref> with ε_j=cj^-2
and Lemma <ref> we get

    ‖P̃_q+2 v_q ‖_H^1_0(Ω) = ‖  P_q+2 V_q ‖_H^1_0(Ω)   ≤1/2(1-ε_q+2) ‖  V_q ‖_H^1_0(Ω)
       ≤1/√(2)(1-ε_q+2) ‖  v_q ‖_H^1_0(Ω) ,

which gives

    | (P̃_q+2v_q, T_q+2^-1P̃_q+2 v_q)_H^1_0(Ω)| ≤(1-ε_q+2)^2 ‖  v_q ‖_H^1_0(Ω)^2  .

Then, from (<ref>) we obtain

    ‖ v ‖_H^1_0(Ω)^2 ≥ε_q+2(2-ε_q+2) ‖  v_q ‖_H^1_0(Ω)^2 ,

which immediately yields (<ref>) for all
v ∈_r=_r^(0,0) and all r>q. 

The same result holds for all other combinations of parity indices; hence, it holds for any v ∈_r. Since polynomials vanishing on ∂Ω form a dense subset of H^1_0(Ω), we conclude that (<ref>) holds for all v ∈ H^1_0(Ω). 







§ PROOF OF THEOREM <REF>



We actually prove the equivalent condition

    ‖ u - u_q ‖_H^1_0(Ω)≲‖ u_q ‖_H^1_0(Ω),

and for that we write

    ‖ u - u_q ‖_H^1_0(Ω)  = sup_v ∈ H^1_0(Ω),  v ≠ 0 (u-u_q, v)_H^1_0(Ω)/‖  v ‖_H^1_0(Ω).

As in the previous section, let us split any v ∈ H^1_0(Ω) as v = v_q + z_q ∈_q ⊕_q^c, where _q^c is given by (<ref>). By the Galerkin orthogonality and the definition of u, we have

    (u-u_q, v)_H^1_0(Ω) = (u-u_q, z_q)_H^1_0(Ω) = (f,z_q)_L^2(Ω) - (u_q, z_q)_H^1_0(Ω) .

Recalling (<ref>) and the condition q > p̂, we have (f,z_q)_L^2(Ω) =0, hence

    (u-u_q, v)_H^1_0(Ω) =  - (u_q, z_q)_H^1_0(Ω) .

Now, recalling (<ref>), we expand the Galerkin solution u_q as u_q =∑_j=4^q U_j. Invoking (<ref>), we get

    (u_q, z_q)_H^1_0(Ω) = (U_q-1+U_q, z_q)_H^1_0(Ω).

Applying Propositions <ref> and <ref>, we get the following bound

    (u-u_q, v)_H^1_0(Ω)≤ C_1C_2 q/q-p‖ u_q ‖_H^1_0(Ω)‖ v ‖_H^1_0(Ω).

Since for q >λ p, the relation q/q-p
<λ/λ-1 holds,
and the proof is complete. 






§ PROOF OF PROPOSITION <REF>



We establish the bound ‖ P_j+2‖_ L(_j,_j+2)≤1/2(1/2 - c/j^2)
for any j ≥ 4. This will be achieved
through various steps: we bound
‖ P_j+2‖_ L(_j,_j+2) 
in section <ref> by 
the ℓ^∞-norm of a suitable matrix; in sections
<ref> and <ref> we characterize such an
ℓ^∞-norm and show that the desired bound reduces to certain
properties of a suitable function; we finally analyze such function
in section <ref>. This analysis is computer assisted.




 §.§ Bounding the operator norm by a matrix norm



Recalling the definition (<ref>) of the subspaces _j as well as Remark <ref>, we can split each _j into its two nontrivial orthogonal components according to parity; precisely, if j is even we have _j= _j^(0,0)⊕_j^(1,1), whereas if j is odd we have _j= _j^(1,0)⊕_j^(0,1). Furthermore, again by orthogonality it holds P_j+2∈ L(_j^ α, _j+2^ α) for any α∈{0,1}^2; hence, our target result can be achieved by considering each parity component separately. Hereafter, we will analyze the case j even and α =(0,0), i.e., we will bound ‖ P_j+2‖_ L(_j^(0,0), _j+2^(0,0)); the other three cases can be treated similarly.

Let us set d_j := dim _j^(0,0), note that d_j = j/2-1
because j=|h|=h_1+h_2 with h_1,h_2≥ 2 even, and let us
introduce the normalized basis functions Φ̂_h := Φ_h /
‖Φ_h ‖_H^1_0(Ω). For the sake of definiteness, let
us order the basis functions in each _j^(0,0) by increasing the
first index h_1. Any v ∈_j^(0,0) is represented as


    v='∑_|h|=jv̂_h Φ̂_h
      with  𝐯:=(v̂_h) ∈ℝ^d_j;

the summation symbol means that only indices h=(h_1,h_2) with even
components are considered. Similarly, any w ∈_j+2^(0,0) is
represented as


    w='∑_|k|=j+2ŵ_k Φ̂_k
      with  𝐰:=(ŵ_k) ∈ℝ^d_j+2.


Therefore, if w=P_j+2 v, then 𝐰= 𝐀^T_j 𝐯,
where 𝐀_j  ∈ℝ^d_j× d_j+2 is the
matrix whose entries are


    a_hk:=(Φ̂_h,Φ̂_k)_H^1_0(Ω)  for   |h|=j,  |k|=j+2,


and  h_1,  h_2,   k_1,  k_2 even (recall that the Φ_h's that span _j^(0,0) form an orthogonal basis for this space).
Note that 𝐀_j is a sub-block of the (even, even) block
𝐀^0,0 of the stiffness matrix 𝐀 for
the normalized Babuška-Shen basis in H^1_0(Ω). More
precisely, denoting by 𝐈_j ∈ℝ^d_j ×
    d_j the identity matrix of order d_j, we have
𝐀^0,0 =  tridiag (𝐀^T_j-2 , 𝐈_j ,  𝐀_j) (with j ≥ 4).

Now, one immediately has 

    ‖ P_j+2‖_ L(_j^(0,0),_j+2^(0,0)) = ‖𝐀^T_j ‖_2 = ‖𝐀_j ‖_2 = ‖ P_j+2^* ‖_ L(_j+2^(0,0),_j^(0,0))

(where ‖·‖_p denotes the p-norm of a matrix), which together with the inequality ‖𝐀_j ‖_2 = 
ρ(𝐀_j 𝐀^T_j)^1/2≤‖𝐀_j 𝐀^T_j ‖_∞^1/2, yields the bound 

    ‖ P_j+2‖_ L(_j^(0,0),_j+2^(0,0))≤‖𝐀_j 𝐀^T_j ‖_∞^1/2 .





 §.§ A first expression for the matrix entries



In order to compute the norm on the right-hand side of
  (<ref>), we observe that 𝐀_j is a
bi-diagonal matrix by condition (<ref>). In fact, for any
index h with |h|=j the only indexes h',h” with |h'|=|h”|=j+2 
that give rise to entries a_h,h' and a_h,h” different from 0 are
h'=h+(0,2) and h”=h+(2,0).
The explicit value of these entries is computable via the following formulas (in which all inner-products and norms are those of H^1_0(Ω)):

    a_h,h'= (Φ̂_h,Φ̂_h') = (Φ_h,Φ_h')/‖Φ_h ‖ ‖Φ_h'‖        
    a_h,h”= (Φ̂_h,Φ̂_h”) = (Φ_h,Φ_h”)/‖Φ_h ‖ ‖Φ_h”‖
 
with

    (Φ_h,Φ_h')    =  -1/(2h_2+1)√((2h_2-1)(2h_2+3)), 
    
    (Φ_h,Φ_h”)    =-1/(2h_1+1)√((2h_1-1)(2h_1+3)),

and

    ‖Φ_h‖^2=2/(2h_1-3)(2h_1+1)+2/(2h_2-3)(2h_2+1),

whence

    ‖Φ_h'‖^2=2/(2h_1-3)(2h_1+1)+2/(2h_2+1)(2h_2+5),


    ‖Φ_h”‖^2=2/(2h_1+1)(2h_1+5)+2/(2h_2-3)(2h_2+1).


Since 2≤ h_1 ≤ j-2 (h_1 even), it is convenient to set n :=j/2 and h_1 :=2i, with 1 ≤ i ≤ n-1; consequently,
h_2:=j-h_1=2(n-i). Substituting these expressions in the previous formulas, we obtain

    (Φ_h,Φ_h')=  -1/(4(n-i)+1)√((4(n-i)-1)(4(n-i)+3)) =: a_i


    (Φ_h,Φ_h”)=-1/(4i+1)√((4i-1)(4i+3)) =: b_i,

and

    ‖Φ_h‖^2=2/(4i-3)(4i+1)+2/(4(n-i)-3)(4(n-i)+1) =: ϕ_i,


    ‖Φ_h'‖^2=2/(4i-3)(4i+1)+2/(4(n-i)+1)(4(n-i)+5) =: ψ_i,


    ‖Φ_h”‖^2=2/(4i+1)(4i+5)+2/(4(n-i)-3)(4(n-i)+1) =: η_i.

Note that b_i=a_n-i and ϕ_i=ϕ_n-i, η_i=ψ_n-i. Hence, for 1 ≤ i ≤ n-1,

    (𝐀_j)_ii= a_i/√(ϕ_i ψ_i) =: δ_i,     
    (𝐀_j)_i,i+1= b_i/√(ϕ_i η_i) = a_n-i/√(ϕ_n-iψ_n-i) = δ_n-i,

i.e., 𝐀_j =  bidiag ( δ_i ,  δ_n-i
  ). Consequently, the nonzero entries of the matrix 𝐀_j𝐀^T_j ∈ℝ^d_j× d_j are

    (𝐀_j𝐀^T_j)_i,i-1= δ_iδ_n-i+1,  
    (𝐀_j𝐀^T_j)_ii= δ_i^2 + δ_n-i^2,  
    (𝐀_j𝐀^T_j)_i,i+1= δ_i+1δ_n-i,

i.e., 𝐀_j𝐀^T_j=  tridiag (δ_iδ_n-i+1, δ_i^2 + δ_n-i^2 ,  δ_i+1δ_n-i ). Let us denote by 
s_i^(j) the sum
of the entries in the i-th row of the matrix 𝐀_j𝐀^T_j, which are all non-negative. Setting for convenience
δ_n=0, we thus have

    s_i^(j)=
    δ_iδ_n-i+1 + δ_i^2 + δ_n-i^2 + δ_i+1δ_n-i,     1 ≤ i ≤ n-1.

It is easily seen that s_i^(j)= s_n-i^(j) for 1 ≤ i ≤n/2. Since

    ‖𝐀_j 𝐀^T_j ‖_∞ = max_1 ≤ i ≤ n-1 s_i^(j) ,

in view of (<ref>) we are left with the problem of proving the existence of a constant C>0 such that

    max_1 ≤ i ≤ n-1 s_i^(j)≤1/4  - C/j^2     for all    j ≥ 4;

indeed, thanks to √(1/4-x)≤1/2 -x for x≤1/4, we obtain Proposition <ref> with c=C.

A direct computation shows that s_1^(j) and  s_n-1^(j)
satisfy the bound in (<ref>) for a suitable C,
because s_1^(j)= s_n-1^(j) < 1/4 for all j ≥ 4 and s_1^(j)→3/28 as j →∞. Thus, in the sequel we focus on the rows indexed from 2 to n-2, for j ≥ 8 (i.e., n ≥ 4).




 §.§ A second expression for the matrix entries



We now apply a change of variables. Observing that all quantities
a_i, ϕ_i, ψ_i, η_i, δ_i defined above
depend upon 4i or 4(n-i) for 2≤ i≤ n-2, we first set
I :=4i and N :=4n ≥ 16.
To introduce the new variables (t,r), we first go back to the original
range 1≤ i ≤ n-1, i.e. 4 ≤ I ≤ N-4, and parametrized
I as follows


    I=4(1-t)+(N-4)t=4+Rt,
        
    0 ≤ t ≤ 1,


with 
R:=N-8 ≥ 8. Similarly, we write


    N-I = 4+Rτ,
        τ=τ(t) := 1-t.


At last, we introduce the second parameter r:=1/R≤1/8. With these notation at hand, we easily obtain the following expressions for  a_i, ϕ_i and ψ_i:

    a_i^2 = 1/R^41/(τ +3r)(τ+5r)^2(τ+7r) =: 1/R^4 A(t, r),


    ϕ_i=1/R^2( 2/(t+r)(t+5r)+2/(τ+r)(τ+5r)) =: 1/R^2 B(t,r),


    ψ_i=1/R^2( 2/(t+r)(t+5r)+2/(τ+5r)(τ+9r)) =: 1/R^2 C(t,r).

Hence, we arrive at

    δ_i^2 = a_i^2/ϕ_i ψ_i = A(t,r)/B(t,r)C(t,r) =: D(t,r).

Straightforward computations show that

    δ_i+1^2 =D(t+4r,r),
        δ_n-i^2 =D(τ,r),
        δ_n-i+1^2 =D(τ+4r,r).

We conclude that the sum of the entries in the i-th row of 𝐀_j𝐀^T_j, given by (<ref>), can be expressed as follows:



    s_i^(j)   = √(D(t,r)D(τ+4r,r)) + D(t,r)
    
        + √(D(t+4r,r)D(τ,r)) + D(τ,r) =: S(t,r)


for  2 ≤ i ≤ n-2, which is equivalent to  4r ≤ t ≤ 1-4r.




 §.§ Bounding the matrix norm



Since the function S(t,r) is symmetric with respect to t=1/2 for any r, we may restrict it to the triangle 0 ≤ t ≤1/2, 
0 ≤ r ≤1/4 t. Fig. <ref> displays two
plots of the function 1/4 - S(t,r), and suggests
clearly that S(t,r) < 1/4 whenever r>0, with a quadratic behavior in r at the origin.
However, establishing such results rigorously is somehow complicated
by the fact that S(t,r) is singular at (t,r)=(0,0), where it becomes multi-valued.

To remove this singularity, we apply the Duffy transform (t,a) ↦ (t,r)=(t,at), which maps the rectangle 0 ≤ t ≤1/2, 0 ≤ a ≤1/4 onto the triangle 0 ≤ t ≤1/2, 0 ≤ r ≤1/4 t.
Correspondingly, we are led to consider the function σ(t,a) := S(t,at), which turns out to be smooth everywhere in this rectangle;
a plot of the function 1/4 - σ(t,a) is depicted in Fig. <ref>.







With the help of a symbolic manipulator, we obtain, for all t ∈ [0,12],

    σ(t,0)=1/4,     ∂σ/∂ a(t,0) = 0,     ∂^2 σ/∂ a^2(t,0) = -G(t)/τ^2,

where

    G(t) := 3t^10 + 9t^8τ^2 - 8t^7τ^3 + 16t^6τ^4 + 24t^5τ^5 + 16t^4τ^6 - 8t^3τ^7 + 9t^2τ^8 + 3τ^10.

We note that the polynomial G(t) is strictly decreasing in
[0,1/2] between G(0)=3 and G(1/2)=1. We thus easily see that
∂^2 σ/∂ a^2(t,0) ≤ -1 for 0 ≤ t ≤1/2; hence, by continuity we get the existence of two constants C_*>0 and  a_* ∈ (0,14] such that ∂^2 σ/∂ a^2(t,a)  ≤  -C_* for 0 ≤ t ≤1/2 and 0 ≤ a ≤ a_*. With these constants at hand, by Taylor's expansion with Lagrange's reminder, we are entitled to write

    σ(t,a)   =  1/4 + ∂^2 σ/∂ a^2(t,a̅) a^2  ≤ 1/4 - C_* a^2 
        for   0≤ t ≤12,  0 < a ≤ a_*,

with some a̅ = a̅(t,a) ∈ (0, a). 



By computing the symbolic expression of the function ∂^2
  σ/∂ a^2(t,a) and by examining its level sets
(via a numerical procedure), one finds that a_* safely satisfies a_* > 1/10 (see Figure <ref>). 
Therefore, going back to our function S(t,r)=σ(t,r/t), we deduce that

    S(t,r) ≤1/4 - C_*/t^2 r^2 ≤1/4 - 4C_* r^2     for   0 < r ≤110t,    t ≤12.

Recalling (<ref>) and using the expressions t=4r(i-1) and r=1/2 1/j-4 > 1/2j, we immediately obtain

    s_i^(j)≤1/4 - C_*/j^2    for   4 ≤ i ≤n2.


Note that we require the restriction i≥ 4 to satisfy the constraint


    r ≤1/10 t = 2/5 r(i-1)
          ⇒    
      i ≥7/2.


Therefore, we are left with the task of establishing a similar bound for
s_2^(j) and s_3^(j) by different means.
It is easily checked that  for j →∞ it holds s_2^(j)→65/308 < 1/4 and s_3^(j)→55/220 < 1/4, while both s_2^(j) and s_2^(j) are <1/4 for all j ≥ 8. This implies the desired bound for a suitable constant C_**>0.   

The proof of (<ref>) is thus complete,
whence Proposition <ref> is established.




§ ACKNOWLEDGEMENTS

The authors wish to thank Valeria Simoncini for helpful suggestions concerning Section 6. The first and fourth authors are members of the INdAM research group GNCS, which granted partial support to this research. The second author has been artially supported by the NSF Grant DMS -1411808, the Institut Henri Poincaré (Paris) and the Hausdorff Institute (Bonn).


99

  
Babuska-1991
I. Babuška, A. Craig, J. Mandel, and J. Pitkäranta.
Efficient preconditioning for the p-version finite element method in
  two dimensions.
  SIAM J. Numer. Anal., 28(3):624Ð–661, 1991.

BW85
R.E. Bank and A. Weiser.
Some a posteriori error estimators for elliptic partial 
differential equations
Math. Comp., 44:285–301, 1985.

BEK:96
F. A. Bornemann, B. Erdmann, and R. Kornhuber,
A posteriori error estimates for elliptic problems in
	two and three space dimensions,
SIAM J. Numer. Anal., 33:1188–1204, 1996.
        
69.1
M. Bürg and W. Dörfler.
Convergence of an adaptive hp finite element strategy in higher
  space-dimensions.
Appl. Numer. Math., 61(11):1132–1146, 2011.

21
P. Binev, W. Dahmen, and R. DeVore.
Adaptive finite element methods with convergence rates.
Numer. Math., 97(2):219 – 268, 2004.

Binev-2015
P. Binev.
Tree approximation for hp-adaptivity.
IMI Preprint Series 2015:07, 2015.

BPS-2009
D. Braess, V. Pillwein, and J. Schöberl.
Equilibrated residual error estimates are p-robust.
Comput. Methods Appl. Mech. Engrg., 198(13-14):1189–1197,
  2009.

CHQZ06
C. Canuto, M. Y.Hussaini, A. Quarteroni, and T. A. Zang.
Spectral Methods. Fundamentals in Single Domains.
Springer 2016.

CNSV16 
C. Canuto, R.H. Nochetto, R. Stevenson, and M. Verani.
Convergence and optimality of hp-AFEM.
Numer. Math. 135 (2017), 1073-1119

CNSV17 
C. Canuto, R.H. Nochetto, R. Stevenson, and M. Verani.
On p-robust saturation for hp-AFEM.
Comput. & Math. with Appl.  73 (2017),  2004–2022.

64.145
L. Demkowicz, J. Gopalakrishnan, and J. Schöberl.
Polynomial extension operators. Part III.
Math. Comp., 81(279):1289–1326, 2012.

DN02
W. Dörfler and R.H. Nochetto, 
Small data oscillation implies the saturation assumption.
Numer. Math. 91:1–12, 2002.

ErnVohralik-2015
A. Ern and M. Vohralík.
Polynomial-degree-robust a posteriori estimates in a unified setting
  for conforming, nonconforming, discontinuous Galerkin, and mixed
  discretizations.
SIAM J. Numer. Anal., 53(2):1058–1081, 2015.

ErnVohralik-2016
A. Ern and M. Vohralík.
Stable broken H1 and H(div) polynomial extensions for
  polynomial-degree-robust potential and flux reconstruction in three space
  dimensions.
hal 01422204, Inria Paris-Rocquencourt, 2016.

GuoBabuska-1
B. Guo and I. Babuška.
The hp version of finite element method, Part 1: The basic
  approximation results.
Comp. Mech, 1:21–41, 1986.

GuoBabuska-2
B. Guo and I. Babuška.
The hp version of finite element method, Part 2: General results
  and application.
Comp. Mech, 1:203–220, 1986.

MW01
J. M. Melenk and B. I. Wohlmuth.
On residual-based a posteriori error estimation in hp-FEM.
Adv. Comput. Math., 15(1-4):311–331, 2002.

Noch93
R.H. Nochetto,
Removing the saturation assumption in a posteriori error
analysis,
Istit. Lombardo Accad. Sci. Lett. Rend. A,
127:67-82, 1993.










