A new mixed functional-probabilistic approach for finite element accuracy
Joël Chaskalovic Franck Assous

Radiation-reaction electromagnetic fields in metasurfaces, a complete description of their optical properties
    Michele Merano
    December 30, 2023
=============================================================================================================





The aim of this paper is to provide a new perspective on finite element accuracy. Starting from a geometrical reading of the Bramble-Hilbert lemma, we recall the two probabilistic laws we got in previous works that estimate the relative accuracy, considered as a random variable, between two finite elements P_k and P_m, (k < m). Then, we analyze the asymptotic relation between these two probabilistic laws when the difference m-k goes to infinity. New insights which qualified the relative accuracy in the case of high order finite elements are correspondingly obtained. 


keywords: Error estimates, Finite elements, Bramble-Hilbert lemma, Probability.





§ INTRODUCTION


Error estimates play a crucial role in advancement of finite element methods. The development and the actual use of a given numerical method is often, if not always, guided by its performance in terms of reliability and accuracy. For this reason, it is still an active subject of research for a wide range of applied mathematicians.
Since the seminal papers of Strang&Fix <cit.>, Ciarlet&Raviart <cit.>,  Babuska<cit.>, Bramble&Hilbert <cit.>, with co-workers, a plethora of work has been published to elaborate and improve the error estimates in various configurations. The goal is to find bounds for the error u-u_h, between the exact solution u of a partial differential equation and its finite element approximation u_h.
Usually, the main property of error estimates which is considered concerns the rate of convergence of a given finite element. In general, these estimates tell us that the finite element error u-u_h, for a given chosen norm, is O(h^k), where h is the mesh size, namely the largest diameter of the elements in a given mesh, and k a positive integer. As the constant involved in O(h^k) is in most cases unknown, it very seldom considered in the analysis.
The aim of these estimates is generally to give a measure of the efficiency of the considered finite element method, and tell us how fast the error decreases as we decrease the mesh size h. However, in these estimates, this constant depends, among others, on the basis functions of the concerned finite element method, and on a semi-norm of the exact solution u, (see for instance <cit.>).
Furthermore, quantitative uncertainties are commonly produced in the mesh generation so that quantitative uncertainties also exist in the approximate solution u_h. For this reason, we have considered the approximation error as a random variable <cit.>, and we aimed to evaluate the probability of the difference between two approximation errors u- u_h^(k) and u- u_h^(m), for a suited functional norm, corresponding to Lagrange finite elements P_k and P_m, (k<m). It is the reason why we introduced in <cit.> a probabilistic framework to compare the relative accuracy between these two finite elements.
This paper is mainly devoted to the asymptotic relation between the two probabilistic laws we derived in <cit.>. Amongst other, it will highlight the relative accuracy between high order finite elements. To the best of our knowledge, this is the first time that such mixed functional and probabilistic approaches are combined to provide new perspectives on finite element accuracy. 
The paper is structured as follows: Section <ref> summarizes the results of <cit.> necessary for understanding the rest of our analysis: the geometrical interpretation of error estimates and the probabilistic laws we got for the relative finite element accuracy.  In Section <ref> properties of P_k basis polynomials are derived whereas Section <ref> is devoted to the asymptotic relationship between the two probability laws of Section <ref>. Concluding remarks follow.



§ GEOMETRICAL INTERPRETATION OF ERROR ESTIMATES AND RELATED PROBABILISTIC LAWS


We consider an open bounded and non empty subset Ω of ^n, and Γ its boundary assumed to be C^1-piecewise. Let u be the solution to the second order elliptic variational formulation:


    (VP){[                   u ∈   V; a(u,v) = l(v),   ∀ v ∈ V, ].


where V is a given Hilbert space endowed with a norm ._V, a(·,·) is a bilinear, continuous and V-elliptic form defined on V × V, and l(·) a linear continuous form defined on V.
Classically, variational problem (VP) has, one and only one, solution u ∈ V (see for example <cit.>). In this paper, we will restrict ourselves to the case where V is the usual Sobolev space of distributions H^1(Ω). 
Let us also consider an approximation u_h of u, solution  to the approximate variational formulation:


    (VP)_h{[                        u_h∈   V_h; a(u_h,v_h) = l(v_h),  ∀ v_h∈ V_h, ].


where V_h is a given finite-dimensional subset of V. 
To state a corollary of Bramble-Hilbert's lemma and a corresponding error estimate, we follow <cit.> or <cit.>, and we assume that Ω is exactly covered by a mesh 𝒯_h composed by N_K n-simplexes K_μ, (1 ≤μ≤ N_K), which respects classical rules of regular discretization, (see for example <cit.> for the bidimensional case and <cit.> in ^n). Moreover, we denote by P_k(K_μ) the space of polynomial functions defined on a given n-simplex K_μ of degree less than or equal to k, (k ≥ 1). 
Then, we remind below the result of <cit.> from which our study is developed: 




Suppose that there exists an integer k ≥ 1 such that the approximation u_h of V_h is a continuous piecewise function composed by polynomials which belong to P_k(K_μ), (1≤μ≤  N_K). 
Then, u_h converges to u in H^1(Ω):


    lim_h→ 0u_h-u_1,Ω=0.


Moreover, if the exact solution u belongs to H^k+1(Ω), we have the following error estimate:


    u_h-u_1,Ω≤𝒞_k h^k   |u|_k+1,Ω ,


where 𝒞_k is a positive constant independent of h, ._1,Ω the classical norm in H^1(Ω) and |.|_k+1,Ω denotes the semi-norm in H^k+1(Ω).




Consider now two families of Lagrange finite elements P_k and P_m corresponding to two values (k,m)∈^*2, (k < m), the corresponding inequalities given by (<ref>), assuming that the solution u to (VP) belongs to H^m+1(Ω), are:




    u^(k)_h-u_1,Ω   ≤   𝒞_k h^k  |u|_k+1,Ω, 
    
    u^(m)_h-u_1,Ω   ≤   𝒞_m h^m  |u|_m+1,Ω ,


where u^(k)_h and u^(m)_h respectively denotes the P_k and P_m Lagrange finite element approximations of u.


Now, if one considers a given mesh for the finite element P_m that contains the mesh processed for the P_k approximation, then for the particular class of problems where (VP) is equivalent to a minimization formulation (MP), (see for example <cit.>), one can show that the approximation error for P_m is always lower than the one for P_k, and P_m is more accurate than P_k, for all values of the mesh size h.
In this paper we consider a more general where, for a given mesh size h, two independent meshes for P_k and P_m are built by a mesh generator. So, usually, to compare the relative accuracy between these two finite elements, one asymptotically considers inequalities (<ref>) and (<ref>) to conclude that, when h goes to zero, P_m is more accurate that P_k, as h^m goes faster to zero than h^k. 
However, for each application, h has a fixed value and this way of comparison is no longer valid. For this reason, our viewpoint will be to determine the relative accuracy between P_k and P_m, (k<m), for a value of h corresponding to two independent meshes.
To this end, let us set:


    C_k = 𝒞_k |u|_k+1,Ω C_m = 𝒞_m |u|_m+1,Ω.


Therefore, instead of (<ref>) and (<ref>), we consider in the sequel the two following inequalities:

    u^(k)_h-u_1,Ω   ≤    C_k h^k, 
    
    u^(m)_h-u_1,Ω   ≤    C_m h^m.


So, we proposed in <cit.> a geometrical interpretation of (<ref>)-(<ref>) which enabled us to consider the values of u^(k)_h-u_1,Ω and u^(m)_h-u_1,Ω as two random variables, respectively denotes X^(k) and X^(m), whose values belong to [0, C_i h^i], (i=k  i=m), due to (<ref>)-(<ref>).
Then, we derived two probabilistic laws of the event:


    {X^(m)(h) ≤ X^(k)(h)}≡{u^(m)_h-u_1,Ω≤u^(k)_h-u_1,Ω},


which corresponds to the relative accuracy between the two finite elements P_k and P_m, for a given value of the mesh size h.
More precisely, let us introduce the two random events A and B as follows:




    A    ≡   {u^(m)_h-u_1,Ω≤u^(k)_h-u_1,Ω}, 
    
    B    ≡   {u^(k)_h-u_1,Ω∈ [C_m h^m,C_k h^k]},    h<h^*,


where h^* is defined by:


    h^* ≡( C_k/C_m)^1/m-k.


Then we showed in <cit.> the two following results:


Let A and B be the two events defined by (<ref>) and (<ref>) and let us assume they are independent. Then, the probability law of the event {X^(m)(h) ≤ X^(k)(h)} is given by:





    Prob{ X^(m)(h) ≤ X^(k)(h)} =  |
    [            1 0 < h < h^*,;            0      h> h^*. ].




Furthermore, if we replace the hypothesis of independency between A and B by considering the two random variables X^(i), (i=k  i=m), independent and uniformly distributed [0, C_i h^i], we also proved in <cit.> the following theorem:


Let u be the solution to the second order variational elliptic problem (VP) defined in (<ref>) and u^(i)_h, (i=k  i=m, k<m), the two corresponding Lagrange finite element P_i approximations, solution to the approximated formulation (VP)_h defined by (<ref>).


We assume the two corresponding random variables X^(i)(h), (i=k  i=m), are independent and uniformly distributed on [0, C_i h^i], where C_i are defined by (<ref>). 


Then, the probability of the event { X^(m)(h) ≤ X^(k)(h)} is given by:


    Prob{ X^(m)(h) ≤ X^(k)(h)} =  |
    [ 1 - 1/2(h/h^*)^m-k       0 < h ≤ h^*,;     1/2(h^*/h)^m-k           h ≥ h^*. ].




The global shapes of the probabilistic laws (<ref>) and (<ref>) are plotted in Figure <ref> and new features of the relative finite elements accuracy are described in <cit.>. Amongst other, these laws clearly indicate that there exist cases, (if h>h* then Prob{ X^(m)(h) ≤ X^(k)(h)}≤ 0.5), where P_m finite elements probably must be overqualified and a significant reduction of implementation time and execution cost can be obtained without loss of accuracy by implementing P_k finite element.


Another interesting property concerns the "sigmoid" probability law (<ref>) and its relationship with the two steps law (<ref>). To prove this relationship, as a first step, we need some new features of the Lagrange finite element P_k. This is the purpose of the next section.



§ P_K CANONICAL BASIS ESTIMATES


In this section we follow the definitions and properties of the P_k finite element in ^n described by P. A. Raviart and J. M. Thomas in <cit.>. 
Let us consider K ⊂^n a n-simplex which belongs to a regular mesh 𝒯_h. Since a complete polynomial of order k which belongs to P_k(K) contains


    N ≡(
      [ n+k;   n ])
    = (n+k)!/n! k!


terms, each n-simplex of the mesh 𝒯_h must be associated to N independent degrees of freedom to assure the unisolvence of the finite element. 
It is convenient to carry out all analysis of n-simplexes in terms of the so-called barycentric coordinates λ_1,…,λ_n+1 which satisfy ∑_i=1^n+1λ_i=1.
A regularly spaced set of points M_i_1,…,i_n+1 may be defined in a n-simplex by the barycentric coordinates values, namely:


    M_i_1,…,i_n+1 = (i_1/k,…,i_n+1/k),  0 ≤ i_1,…,i_n+1≤ k,


satisfying:


    i_1 + … +i_n+1=k.


One can verify that the number of points defined by (<ref>)-(<ref>) is equal to N, the dimension of P_k(K) in (<ref>). 
Therefore, we introduce the canonical basis of functions p_i_1,…,i_n+1 of the variables (λ_1, …, λ_n+1) which belongs to P_k(K) defined by:


    p_i_1, …, i_n+1(λ_1, …, λ_n+1) ≡∏_j=1^n+1P_i_j(λ_j),

where the auxiliary polynomial P_i_j(λ_j) is given by:


    P_i_j(λ_j) ≡ |
    [ ∏_c_j=1^i_j(k λ_j - c_j +1/c_j),                         i_j ≥ 1,;                               1,                         i_j = 0. ].


P_i_j is clearly a polynomial of order i_j in λ_j, and therefore, due to condition (<ref>), p_i_1, …, i_n+1 given by (<ref>) is a polynomial of order k. 
In the sequel, notice that we will also use a simple index numbering to substitute the multi-index numbering. It will be the case for the N points M_i_1,…,i_n+1 simply denoted (M_i)_i=1,N, as well as for the N canonical functions p_i_1,…,i_n+1 denoted (p_i)_i=1,N, and so on. 
Therefore, the main property of the canonical basis of functions p_i ≡ p_i_1,…,i_n+1 defined in (<ref>) is that for a given set of N values φ_i ≡φ_i_1, …, i_n+1 known at the N points M_i ≡ M_i_1, …, i_n+1, the polynomial Q in P_k(K) given by:


    ∀ M ∈ K: Q(M)     =     Q(λ_1, …, λ_n+1) 
        =    ∑_i_1 + … + i_n+1=kφ_i_1, …, i_n+1 p_i_1, …, i_n+1 (λ_1, …, λ_n+1) 
        =    ∑_i=1^Nφ_i p_i(λ_1, …, λ_n+1),


is the unique one in P_k(K) such that Q(M_i)= φ_i. 

The following result concerns the features of the canonical basis (p_i)_i=1,N defined by (<ref>), where N is given by (<ref>), regarding the semi-norm |.|_m,p,K in W^m,p(K) in the particular cases p=2 and m=0 or m=1.
First of all, we remark that the structure of the elementary polynomials P_ij defined by (<ref>) looks like the numerator of the famous Lagrange polynomials. Then, we will establish the first estimate.


Let [a,b], (a<b), be a given interval and N_p a given non-zero integer. We consider a set of N_p+1 uniform distributed points x_j, (j=0,…,N_p), in [a,b] defined by:

    ∀ j=0,…,N_p: x_j=a+jh, h=b-a/N_p.

Let also Π be the function defined on [a,b] by: Π(x)= ∏_j=0^N_p(x-x_j). 
Then,

    ∀ x ∈ [a,b]: |Π(x)| ≤ (N_p+1)! h^N_p+1.




 / 
First, remark that for x=b, Π(b) vanishes, so that (<ref>) is satisfied. Hence, let x be a fixed value in the interval [a,b[. It exists a unique i∈{0,…,N_p-1} such that x∈[x_i,x_i+1[. Therefore, we write the function Π(x) as follows:


    Π(x) = (x-x_0)…(x-x_i)(x-x_i+1)…(x-x_N_p).


Furthermore, we have the following inequalities:


    ∀ j=0,…,i : |x-x_j| ≤ (i+1-j)h, 
          ∀ j=1,…,N_p-i : |x-x_ i+j| ≤ jh ≤ (i+j+1)h.


Therefore, Π(x) written in (<ref>) can be controlled by the help of (<ref>) and (<ref>) by:


    |Π(x)| ≤[(i+1)h×(ih)… (2h)× (h)]×[(i+2)h×(i+3)h…(N_p+1)h],


and after reorganizing the right side of (<ref>), we get (<ref>).



The following lemma gives us the first point-to-point estimates for the polynomials p_i defined by (<ref>).



Let p_i, (i=1,…,N), be the basis functions of the space of polynomials P_k(K) which are defined by (<ref>)-(<ref>). 
Then, ∀ i=1,…,N,  ∀ l=1,…,n+1:


    |p_i(λ_1, …, λ_n+1)|≤ k^n+1, |∂ p_i/∂λ_l(λ_1, …, λ_n+1)|≤ k^n+2.




 /

Let us introduce the integer n_i, (0 ≤ n_i≤ n+1), which corresponds to the number of polynomials P_i_j(λ_j) such that:




    ∀ j=1,…,n_i,  (n_i ≥ 1),     :     P_i_j(λ_j)=P_1(λ_j)=kλ_j, (i_j=1), 
    ∀ j=n_i+1,…,n+1,  (n_i ≤ n)     :     P_i_j(λ_j)=kλ_j(kλ_j - 1)…(kλ_j - i_j+1)/i_j!, (i_j>1).


When n_i=0, then p_i has the following structure:


    p_i(λ_1,…,λ_n+1) = ∏_j=1^n+1(kλ_j(kλ_j - 1)…(kλ_j - i_j+1)/i_j!), (i_j >1, ∀ j=1,…,n+1),


and when n_i=n+1, then p_i corresponds to:


    p_i(λ_1,…,λ_n+1) = ∏_j=1^n+1(kλ_i_j).


∙ Let us begin by fixing a given value of n_i, (1 ≤ n_i ≤ n).
Concerning the control of the polynomials p_i, we split it into two groups of elementary polynomials P_i_j as follows:


    p_i(λ_1,…,λ_n+1) = ∏_j=1^n_i(kλ_i_j).∏_j=n_i+1^n+1(kλ_j(kλ_j - 1)…(kλ_j - i_j+1)/i_j!).


Now, on the first hand, the barycentric functions λ_j, (j=1…,n+1) satisfy: 
    ∀ M∈ K: 0≤λ_j(M)≤ 1.


On the other hand, by applying lemma <ref> by setting x=kλ_j and h=1, we have the following estimate:


    ∀ j=n_i+1,…,n+1: |P_i_j(λ_j)| = |kλ_j(kλ_j - 1)…(kλ_j - i_j+1)/i_j!| ≤  1,


and finally,

    |p_i(λ_1,…,λ_n+1)| ≤|∏_j=1^n_i(kλ_i_j)| ≤ k^n_i≤ k^n,


as n_i ≤ n.
Let us now consider the partial derivative ∂ p_i/∂λ_l, for a given pair of non zero integers (i,l). 
By (<ref>) we can write the concerned partial derivative as:


    ∂ p_i/∂λ_l =P_i_1…∂ P_i_l/∂λ_l… P_i_n+1.


Thus, two cases have to be considered. The first one corresponds to the case when P_i_l is a single monomial (i_l=1): P_i_l(λ_l)= kλ_l.
Therefore, (<ref>) gives:


    ∂ p_i/∂λ_l =k P_i_1… P_i_l-1P_i_l+1… P_i_n+1,


and similarly to (<ref>)-(<ref>), we get the following estimate:


    |∂ p_i/∂λ_l| ≤ k.k^n_i-1≤ k^n.


Let us now consider the case when i_l>1. It means that the polynomial P_i_l has the structure of (<ref>), composed at least by two monomials. Then, its partial derivative with respect to λ_l is equal to:


    ∂ P_i_l/∂λ_l =1/i_l ![/k(kλ_l-1)…(kλ_l-i_l+1) + … + kλ_l(kλ_l-1)…(kλ_l-i_l+2)k].


So, by using the same arguments we implemented to upper bound the function Π defined in lemma <ref>, on can increase each term of the right hand side of (<ref>) to finally obtain:


    |∂ P_i_l/∂λ_l| ≤k/i_l ![/  i_l! i_l ]≤ k^2,


as ∀ l=1 to n+1: i_l ≤ k.
Finally we get the estimate for the partial derivative of p_i with respect to λ_l:


    |∂ p_i/∂λ_l| ≤ k^2   k^n_i≤ k^n+2.


Let us now consider the two cases when n_i=0 or n_i=n+1.
∙ If n_i=0 then, due to (<ref>) we have:


    |p_i(λ_1,…,λ_n+1)| ≤  1 ≤ k^n,


thanks to (<ref>).
In the same way, we have the following inequalities:


    |∂ p_i/∂λ_l| ≤ k^2 ≤ k^n+2,


where we used (<ref>)-(<ref>) which correspond to the present case.
∙ If n_i=n+1 then p_i is given by (<ref>) and we have:


    |p_i(λ_1,…,λ_n+1)| ≤ k^n+1|∂ p_i/∂λ_l| ≤ k^n+1,


due to the basic barycentric functions features.
Therefore, from all the above upper bounds, we get (<ref>).



We can now get two estimates of the canonical basis (p_i)_i=1,N with respect to the semi-norms |p_i|_0,2,K and |p_i|_1,2,K.



Let (p_i)_i=1,N be the canonical basis defined in (<ref>). Then, if k >n/2, we have:


    ∑_i=1^N|p_i|_0,2,K = 𝒪(/k^n+1 (k+n)^n) ∑_i=1^N|p_i|_1,2,K = 𝒪(/k^n+2 (k+n)^n)/ρ_K,


where 𝒪 denotes Landau's notation and ρ_K the diameter of the largest inscribed sphere within K.




 / 
▸ Let us begin with the estimate of p_i with respect to the semi-norm |.|_0,2,K. 
From the local estimate of p_i given by (<ref>), we directly get the |.|_0,2,K-semi-norm for each polynomial p_i, (1 ≤ i≤ N), as follows:


    |p_i|_0,2,K≤√((K)) k^n+1,


and by aggregating on all the N basis function p_i, we get:


    ∑_i=1^N|p_i|_0,2,K   ≤   √((K))   k^n+1 (n+k)!/n! k!, 
       ≤   √((K)) (k+n)^n k^n+1,


as n≥ 1 and where we used the value of N corresponding to the dimension of the space P_k(K) given by (<ref>). 
Finally, with (<ref>) we get the first estimate of (<ref>). 

▸ Let us prove now the second estimate of (<ref>) with respect to the semi-norm |.|_1,2,K. 
Due to remark 2.2 in R. Arcangeli and J. L. Gout <cit.>, for each canonical basis function p_i, if k > n/2, we have:


    |p_i|_1,2,K≤1/ρ_K{∫_K[∑_j=1^n|∂ p_i/∂ x^j(x)|]^2dx}^1/2,


where ρ_K is the supremum of the diameters of the inscribed spheres within the n-simplex K.
Moreover, each partial derivative ∂ p_i/∂ x^j can be computed using the chain rule as follows:


    ∂ p_i/∂ x^j = ∑_l=1^n+1∂ p_i/∂λ_l∂λ_l/∂ x^j,


where each partial derivative ∂λ_l/∂ x^j is a constant Λ^(l)_j that does not depend on k,  λ_l being a polynomial of degree at most equal to one. 
So, we have:


    ∂ p_i/∂ x^j = ∑_l=1^n+1Λ^(l)_j∂ p_i/∂λ_l.


Consequently, from (<ref>), we get:


    |p_i|^2_1,2,K   ≤   1/ρ^2_K∫_K[∑_j=1^n|∑_l=1^n+1Λ^(l)_j∂ p_i/∂λ_l|]^2dx
    
    
       ≤   ( nΛ/ρ_K)^2∫_K|∑_l=1^n+1∂ p_i/∂λ_l|^2dx, 
    
       ≤   (K)( n(n+1)Λ  k^n+2/ρ_K)^2,


where we set Λ≡max_1 ≤ j ≤ n 
 1 ≤ l ≤ n+1|Λ^(l)_j| and due to (<ref>).
By aggregating (<ref>) on the N basis functions p_i we finally get:


    ∑_i=1^N|p_i|_1,2,K≤√((K))  n(n+1)Λ/ρ_K  (k+n)^n k^n+2,


which corresponds to the second estimate of (<ref>).




We notice that for applications the condition k > n/2 holds for the dimension n=1 when k≥ 1, but if n=2 or n=3 this requires k≥ 2. Consequently, the case of the finite element P_1 could be still considered by using other results of <cit.>, as we will mention later, (see Theorem <ref>).


The two estimates (<ref>) will now be used to determine the asymptotic behavior of the probability distribution (<ref>) of theorem <ref>.



§ ASYMPTOTIC LIMIT OF THE "SIGMOID" PROBABILITY DISTRIBUTION


As we already mentioned, the probability distribution (<ref>) has been approximated by the stepwise law (<ref>) if one assumes the independency between the events A and B defined by (<ref>) and (<ref>). 
Conversely, here we will study the behavior of the non linear law (<ref>) when q ≡ m-k goes to infinity. This study is not only theoretical. It is clearly related to the well-known question, namely,  in which way high order finite element methods can solve partial differential equations more efficiently than low order methods. More precisely, how large of a polynomial degree is beneficial? Here again, we have chosen to treat the problem via a probabilistic approach, handling the uncertainties (randomness of the data, of the mesh, etc.) by random variables. Note that often in the applications, one considers cases when k=1 or 2, whereas the high order degree m is around 20-25, see for instance <cit.>. 
More precisely, let us give a fixed value of k. 
Then, we define the sequence of functions (𝒫_q(h))_q ∈^⋆ by:


    𝒫_q(h) ≡ Prob{ X^(k+q)(h) ≤ X^(k)(h)},


where Prob{ X^(k+q)(h) ≤ X^(k)(h)}, following (<ref>), is given by:


    Prob{ X^(k+q)(h) ≤ X^(k)(h)} =  |
    [ 1 - 1/2(h/h_q^*)^q    0 < h/h_q^*≤ 1,;     1/2(h_q^*/h)^q        h/h_q^*≥ 1, ].


and where (h^*_q)_q∈^⋆ is the sequence defined by:


    h^*_q≡(C_k/C_k+q)^1/q.


As one can see the critical value h^*_q strongly depends on q, among others, by the constant C_k+q.
To this end, we will firstly determine an estimate of the constant 𝒞_k defined by (<ref>), relatively to the finite element P_k.

This is the purpose of the following theorem:


For k > n/2, let 𝒞_k be the constant introduced in the error estimate (<ref>). Then, the following estimation holds: 

    𝒞_k = 𝒪((k+n)^n k^ n+2/(k-1)! (k-n/2)).






The proof of this theorem is based on the paper of R. Arcangeli and J.L. Gout <cit.>,  itself an extension of the one of P.G. Ciarlet and P.A. Raviart <cit.>. 
To this end, let us firstly recall the conditions of theorem 2.1 of R. Arcangeli and J.L. Gout. 
Let Ω be an open bounded and non empty convex subset of ^n and Γ its Lipschitz boundary. We assume that Σ={a_i}_i=1,N is a P-unisolvent set of points which belong to Ω̅, where P denotes a space of finite dimension such that P_k ⊂ P ⊂ C^k(Ω̅), and P_k the space of polynomial functions of degree less than or equal to k. 
Then, for all u∈ W^k+1,p(Ω) and for all integer ν≥ 0 such that


    k+1 > ν + n/p,


we have:


    |u -Π_h u|_ν,p,Ω   ≤   1/(k-ν)!(k+1-ν-n/p)|u|_k+1,p,Ω h^k+1-ν
        +    (∑_i=1^N|p_i|_ν,p,Ω)/[(Ω)]^1/p k!(k+1-n/p)|u|_k+1,p,Ω h^k+1,


where |.|_ν,p,Ω denotes the usual semi-norm in the Sobolev spaces W^ν,p(Ω), Π_h the classical Lagrange interpolation which consists to interpolate the set of points Σ in ^n by a polynomial function of a given degree k, and (p_i)_i=1,N the unique functions such that


    p_i(Mj)=δ_ij, ∀  M_j ∈Σ, ∀  1≤ i, j≤ N,


where δ_ij denotes the Krönecker's symbol. 
Here, for our objectives, we write (<ref>) in the particular following conditions:


  * Ω=K_μ, where K_μ, (1 ≤μ≤ N_K), is a n-simplex which belongs to a given regular mesh 𝒯_h.

  * u is the exact solution to the variational formulation (VP) defined in (<ref>).

  * The set of points Σ in ^n corresponds to the P_k finite element degrees of freedom defined on the n-simplex K_μ defined by (<ref>) and (<ref>).

  * The interpolation operator Π_h is replaced by Π_K_μ, the local Lagrange interpolation operator.


Then, we choose in (<ref>) p=2, ν=0 and ν=1 which implies that k > n/2 due to (<ref>), or equivalently, k ≥ 2 for a problem set in dimension n≥ 2.
The case of the finite element P_1 in dimension n≥ 2 could also be considered by adapting our theorem with another result from R. Arcangeli and J.L. Gout (see remark 2.3 and theorem 1.1 in <cit.>). 
So, we get the following inequalities:
▸ For ν=0 we have:


    ∀  K_μ∈𝒯_h, 1 ≤μ≤ N_K     :






    |u -Π_K_μu|_0,K_μ   ≤   1/k!(k+1-n/2)|u|_k+1,K_μ h_K_μ^k+1
        +    (∑_i=1^N|p_i|_0,K_μ)/[(K_μ)]^1/2 k!(k+1-n/2)|u|_k+1,K_μ h_K_μ^k+1,


which becomes:


    |u -Π_K_μu|_0,K_μ≤[1 + (k+n)^n k^n+1/k!(k+1-n/2)]|u|_k+1,K_μ h_K_μ^k+1,


due to (<ref>). 


▸ In the same way, for ν=1, we have:


    |u -Π_K_μu|_1,K_μ   ≤   1/(k-1)!(k-n/2)|u|_k+1,K_μ h_K_μ^k
        +    (∑_i=1^N|p_i|_1,K_μ)/[(K_μ)]^1/2k!(k+1-n/2)|u|_k+1,K_μ h_K_μ^k+1.


which leads to:


    |u -Π_K_μu|_1,K_μ   ≤   1/(k-1)!1/(k-n/2)|u|_k+1,K_μ h_K_μ^k
    
        +    n(n+1)Λ/ρ_K_μ(k+n)^n k^n+2/k! (k+1-n/2)|u|_k+1,K_μ h_K_μ^k+1,


due to (<ref>), and finally:


    |u -Π_K_μu|_1,K_μ≤[1+ σ n(n+1) Λ(k+n)^n k^n+2/(k-1)! (k-n/2)] |u|_k+1,K_μ h_K_μ^k,


where we introduced the number σ≥ 1 such that h_K_μ/ρ_K_μ≤σ, ∀ K_μ∈𝒯_h.
Therefore, by the help of (<ref>) and (<ref>), we get the following estimate of the local interpolation error with respect to the H^1-norm, using that the mesh 𝒯_h is regular, and by setting h ≡max_K_μ∈𝒯_h h_K_μ:


    u -Π_K_μ u^2_1,K_μ   ≤   [1 + (k+n)^n k^n+1/k!(k+1-n/2)]^2
    |u|^2_k+1,K_μ h_K_μ^2(k+1)
        +    [1+ σ n(n+1) Λ(k+n)^n k^n+2/(k-1)! (k-n/2)]^2 |u|^2_k+1,K_μ h_K_μ^2k.


Then, we get:


    u -Π_K_μ u_1,K_μ   ≤   [C(Ω̅,σ,Λ,n)(k+n)^n k^n+2/(k-1)! (k-n/2)] |u|_k+1,K_μ h^k,


where we introduced the constant C(Ω̅,σ,Λ,n) defined by:


    C(Ω̅,σ,Λ,n) ≡ 1+ 2 (Ω̅) + σ n(n+1) Λ.


Therefore, by the help of (<ref>), we get for the whole domain Ω the following estimate of the interpolation error:


    u -Π_h u_1,Ω    =    (∑_K_μ∈𝒯_hu -Π_K_μu_1,K_μ^2)^1/2
    
       ≤    C(Ω̅,σ,Λ,n)(k+n)^n k^n+2/(k-1)! (k-n/2)(∑_K_μ∈𝒯_h|u|^2_k+1,K_μ)^1/2h^k,
    
       ≤    C(Ω̅,σ,Λ,n)(k+n)^n k^n+2/(k-1)! (k-n/2) |u|_k+1,Ω h^k.


Then, inequality (<ref>) leads to the estimate (<ref>) if one takes into account the estimate of Céa's lemma <cit.>. Indeed, consider the H^1-norm to measure the difference between the exact solution u to the variational problem (<ref>) and its approximation u_h solution to (<ref>), we have:


    u-u_h_1,Ω≤M/α u-Π_h u_1,Ω,


where M  is the continuity constant and α the ellipticity constant of the bilinear form a(·,·). 
As a consequence, by the help of (<ref>) we obtain that the constant 𝒞_k in (<ref>) satisfies:


    𝒞_k ≤MC(Ω̅,σ,Λ,n)/α (k+n)^n k^n+2/(k-1)! (k-n/2),


which corresponds to (<ref>).



For the sequel, we introduce the constant 𝒞^*_k defined by:


    𝒞_k^* ≡MC(Ω̅,σ,Λ,n)/α (k+n)^n k^n+2/(k-1)! (k-n/2),


and the corresponding h^*_q defined in (<ref>), in which we substitute C_k by the corresponding value of C^*_k, that is:


    h^*_q≡(C^*_k/C^*_k+q)^1/q = (𝒞^*_k |u|_k+1,Ω/𝒞^*_k+q|u|_k+q+1,Ω)^1/q.


As we are interested in the asymptotic behavior of h^*_q defined by (<ref>) when q goes +∞, we will assume that the solution u to the variational problem (VP) belongs to H^r(Ω), (∀  r∈).
We are now in position to propose an estimate of the sequence (h^*_q)_q∈^⋆ defined by (<ref>), when q goes to infinity and the corresponding asymptotic limit of the sequence of functions defined by (<ref>).



Let us assume that solution u of problem (VP) belongs to H^r(Ω), (∀  r∈).
Let also (h^*_q)_q∈^⋆ be the sequence defined by (<ref>) and (/𝒫_q(h))_q∈^⋆ the corresponding sequence of functions defined by (<ref>).
For a fixed value of k, (k > n/2), if


    lim_q→ +∞|u|_k+q+2,Ω/|u|_k+q+1,Ω = l, (l∈^*_+),


then,

    h^*_q q → +∞∼1/e l q,   lim_q→ +∞h^*_q = +∞.


Moreover, the sequence of functions (/𝒫_q(h))_q∈^⋆converges pointwise when q goes to +∞ to the function 𝒫_0 defined on ^*_+ by:


    𝒫_0(h)=  |
    [            1 0 ≤ h < +∞ ,;          1/2     h = +∞ . ].





/


▸ Let us replace the expression of 𝒞^*_k defined by (<ref>) in h^*_q given by (<ref>). Then, we have following asymptotic behavior:


    (h^*_q)^q q → +∞∼(k+n)^n k^n+2/(k-1)! (k-n/2) (q+k-1)!(q+k-n/2)/(/q+k+n)^n (q+k)^n+2 .|u|_k+1,Ω/|u|_k+q+1,Ω.


However, from Stirling's formula, when q goes to +∞, we can specify the equivalent of h^*_q  given by inequality (<ref>):


    (q+k-1)!(q+k-n/2)/(q+k)^n+2(/q+k+n)^n   q → +∞∼   √(2π (q+k-1))(q+k-1/e)^(q+k-1)(q+k-n/2)/(q+k)^n+2(/q+k+n)^n, 
    
       q → +∞∼   √(2 π)(q+k-1)^(q+k-1/2)/e^q+k-11/(/q+k)^2n+1, 
    
       q → +∞∼   √(2 π) (q+k)^q+k-2n-3/2/e^q+k.


Then, (<ref>) in (<ref>) leads to:


    (h^*_q)^q q → +∞∼Θ  e^-(q+k)(q+k)^q+k-2n-3/2 .|u|_k+1,Ω/|u|_k+q+1,Ω,


where we introduced the constant Θ independent of q defined by:


    Θ≡√(2 π)(k+n)^n k^n+2/(k-1)! (k-n/2).


Moreover, as we assume condition (<ref>), if we introduce the two sequences (v_q)_q∈ and (w_q)_q∈ as follows:


    ∀  q∈: v_q≡ln |u|_k+q+1,Ω,  w_q ≡ q,


then, the ratio r_q defined by:


    r_q ≡v_q+1-v_q/w_q+1-w_q = ln(|u|_k+q+2,Ω/|u|_k+q+1,Ω),


has a limit L≡ln  l ∈, when q goes to +∞: lim_q→ +∞r_q = L.
As a consequence, due to Stolz-Cesaro theorem <cit.>, the ratio v_q/w_q also has the same limit L when q goes to +∞:


    lim_q→ +∞v_q/w_q = lim_q→ +∞ln |u|_k+q+1,Ω/q = L,


and, |u|_k+1,Ω being a constant with respect to q,


    lim_q→ +∞(|u|_k+1,Ω/|u|_k+q+1,Ω)^1/q = lim_q→ +∞(1/|u|_k+q+1,Ω)^1/q =e^-L=1/l.



As a consequence, from (<ref>) and (<ref>) we conclude that:

    h^*_q q → +∞∼1/e l q,


and


    lim_q→ +∞ h^*_q = +∞.


▸ /
Let us now examine the convergence pointwise of the sequence of functions (𝒫_q(h))_q ∈^⋆ defined in (<ref>). 
To this end let us, for example, consider a fixed value h_0 such that 0 < h_0< h^*_q. Then, due to (<ref>), we have:


    ∀  0 < h_0 < h^*_q: lim_q→ +∞(h_0/h^*_q)^q = lim_q→ +∞e^ qln(h_0/h^*_q) = 0^+,


and similarly for the second part of (<ref>) corresponding to the case h_0 > h^*_q.
Moreover, when h_0 = h^*_q, 𝒫_q(h^*_q) = 1/2, ∀  q ∈^*.
This enables us to define the pointwise limit function 𝒫_0(h) of 𝒫_q(h)  when q goes to +∞ as:

    𝒫_0(h)=  |
    [            1 0 ≤ h < +∞ ,;          1/2     h = +∞ . ].

Remark that, as h goes to infinity in (<ref>), the limit function 𝒫_0(h) has a discontinuity.  This comes from the interchange of limiting in q and in h is illicit, namely:

    1/2= lim_q → +∞[ lim_h → +∞𝒫_q(h)] ≠lim_h ≠ h^*_qh → h^*_q[  lim_q → +∞𝒫_q(h)]   =
    lim_h ≠ h^*_qh → h^*_q1_]0,+∞[(h) =1 .








/





  * In theorem <ref> we assume that the exact solution u to the second order elliptic variational problem (VP) belongs to H^r(Ω), ∀  r ∈. This is typically the case when the linear form l(.) in (<ref>) is defined by a sufficiently regular function denoted f. For example, if f∈ H^r(Ω) then, for a second order elliptic operator, u belongs to H^r+2(Ω) (see for example <cit.>).


  * Even if condition (<ref>) of theorem <ref> seems restrictive for applications, it is not necessary the case. Take for example the following standard problem:


    (VP){[        u ∈ V; -Δ u = f Ω ,;    u =g ∂Ω , ].


where Ω is the open unit square ]0,1[×]0,1[ in ℝ^2, and f(x,y)=2π^2 sin(π x) cos(π y). 
We readily get that u(x,y)=sin(π x) cos(π y) is the exact solution of (<ref>) in  V≡ H^r(Ω), (∀  r∈), provided that the Dirichlet boundary condition g is defined by:

    {[   g(x,0) = sin(π x)         g(0,y) = 0,; g(x,1) = -sin(π x),         g(1,y) = 0. ].


Then, we obtain that the semi-norm |u|_k,Ω is equal to

    ∀ k ≥ 0: |u|_k,Ω = (√(2))^k-2π^k.


Finally, on can check that condition (<ref>) is satisfied in that case, as we have


    lim_q→ +∞|u|_k+q+2,Ω/|u|_k+q+1,Ω = 1/π√(2).









/




  * Theorem <ref> corresponds to an expected behavior. Indeed, when q=m-k tends to infinity, it claims that the event "P_m is more accurate that P_k" is an almost sure event for all positive values of h. In other words, the higher the distance between m and k, the higher the size of the interval [0,h^*_q] where the event "P_m is more accurate that P_k" is an almost sure event.



  * One can notice that this asymptotic feature is also very intuitive in terms of probability. Indeed, as q=m-k goes to infinity,  for h < h^*_q, the probability such that X^(k+q)(h)≤ X^(k)(h) goes to 1, since the domain of existence of X^(k+q) goes to 0 whereas the one of X^(k) stays fixed and finite. On the contrary, when  h > h^*_q the complementary situation has to be taken into account. Namely, the domain of existence of X^(k+q) goes to infinity in comparison with those of X^(k) which stays again fixed and finite, (see Figure <ref>).








§ DISCUSSION AND CONCLUSION


In this paper we applied to high order finite elements the novel probabilistic approach we developed in <cit.> to evaluate the relative accuracy between two Lagrange finite elements P_k and P_m, (k<m). This new way to evaluate the relative accuracy is based on a geometrical interpretation of the error estimate and by considering the approximation errors as random variables. Therefore, we derived two probabilistic laws, the "two steps" one and  the "sigmoid" one, which describe new features of the relative accuracy between finite elements.


The perspectives of this new approach are not restricted to finite element methods but can be extended to other approximation methods: given a class of numerical schemes and their corresponding error estimates, one is able to order them, not only in terms of asymptotic rate of convergence, but also by evaluating the most probably accurate. 
For the finite elements we considered, we can state the following properties as consequences of Theorem <ref>:





  * For the very small values of h, the probability such that "P_m is more accurate than P_k" probability goes to 1. It corresponds to the classical interpretation of Bramble-Hilbert lemma.


  * Depending on the position of h with respect to the critical value h^* defined by (<ref>), P_k or P_m finite elements are more likely accurate.


  * When h is smaller than h^*, P_m finite elements are not only asymptotically better than P_k finite elements as h becomes small, but they are almost surely more accurate for all of these values of h, with a probability between 0.5 and 1.


  * When h is greater than h^*, P_k finite elements are almost surely more accurate than P_m finite elements, even though k<m, with a probability between 0.5 to 1.


This last property upsets the widespread idea regarding the relative accuracy between P_k and P_m, (k< m), finite elements. It clearly indicates that there exist cases where P_m finite elements surely must be overqualified and a significant reduction of implementation time and execution cost can be obtained without a loss of accuracy. We already observed such a phenomenon by using data mining techniques (see <cit.>, <cit.>, <cit.> and <cit.>). 
However, when the difference between k and m becomes large, high order finite elements are concerned and on can raise the question if P_k finite element is still almost surely more accurate than P_m when h≥ h^*, or at least, on which interval of h this would still be true. 
It is the purpose of the results presented in this paper. To achieve these objectives, we first remark that we need to get asymptotic information regarding the critical value of h^*_q defined in (<ref>) to evaluate the limit of the corresponding sequence of probabilities (/𝒫_q(h))_q∈^⋆ defined in (<ref>). 
As a consequence, we needed to also get an estimate of the asymptotic behavior of the constant 𝒞_k which is involved in the error estimate (<ref>) derived from Bramble Hilbert Lemma.
So, based on the reference work of Arcangeli and Gout <cit.>, we identified the necessity to get several estimates which concerned the canonical basis of polynomials p_i introduced in (<ref>), related to a given finite element P_k . This is the purpose of Lemma <ref>. 
Then, we got an explicit k-dependency of the constant 𝒞_k in Theorem <ref> which enabled us to determine the asymptotic behavior of h^*_q, and accordingly, for the sequence of probabilities (/𝒫_q(h))_q∈^⋆ defined in (<ref>). This is the purpose of Theorem <ref>. 
This theorem claims that when one considers two finite elements P_k and P_m for a fixed value of k, when m goes to infinity, it does not exist anymore any interval or any value of the mesh size h such that the finite element P_k could be almost surely more accurate than P_m. 
In other words, when one implements high order finite elements P_m, the choice of the mesh size h has not to be done to guarantee the better accuracy of the concerned finite element P_m in comparison with another cheaper one P_k, (k<<m).
Acknowledgements:

The authors want to warmly dedicate this research to pay homage to the memory of Professors André Avez and Gérard Tronel who largely promote the passion of research and teaching in mathematics.

100

Arcangeli_Gout R. Arcangeli, J.L. Gout, Sur l'évaluation de l'erreur d'interpolation de Lagrange dans un ouvert de ^n, ESAIM: Mathematical Modelling and Numerical Analysis - Modélisation Mathématique et Analyse Numérique, 10, pp. 5–27 (1976).

AsCh11 F. Assous, J. Chaskalovic, Data mining techniques for scientific computing: Application to asymptotic paraxial approximations to model ultra-relativistic particles, J. Comput. Phys., 230, pp. 4811–4827 (2011).

AsCh13 F. Assous, J. Chaskalovic, Error estimate evaluation in numerical approximations of partial differential equations: A pilot study using data mining methods, C. R. Mecanique 341 (2013) 304–313.

CMAM_2019 J. Chaskalovic, F. Assous, A new probabilistic interpretation of Bramble-Hilbert lemma, accepted in Computational Methods in Applied Mathematics, (2018).

Brezis H. Brezis, Analyse fonctionnelle - Thorie et applications, Masson (1992).

ChaskaPDE J. Chaskalovic, Mathematical and numerical methods for partial differential equations, Springer Verlag, (2013).

AsCh16 J. Chaskalovic, F. Assous, Data mining and probabilistic models for error estimate analysis of finite element method, Maths. And Comp. in Simulation 129 (2016) 50–68.

AsCh17 J. Chaskalovic, F. Assous, Probabilistic approach to characterize quantitative uncertainty in numerical approximations, Maths. Model. and Anal. 22 (2013) 106–120.

Ciarlet P.G. Ciarlet, Basic error estimates for elliptic problems, in Handbook of Numerical Analysis, Vol. II, Eds. P.G. Ciarlet and J. L. Lions, North Holland, (1991).

Ciarlet_Raviart Ciarlet, Raviart, General Lagrange and Hermite interpolation in ^n with applications to finite element methods, Arch. Rat. Mech. Anal. 46 (1972) 177–199.

Mitc15 W.F. Mitchell, How High a Degree is High Enough for High Order Finite Elements? Procedia Computer Science 15 (2015) 246–255.

OviFur O. Furdui, Limits, Series, and Fractional Part Integrals, Springer (2013).

RaTho82 P.A. Raviart, J.M. Thomas, Introduction à l'analyse numérique des équations aux dérivées partielles, Masson (1982).

BrHi70 J.H. Bramble, S.R. Hilbert, Estimation of linear functional on Sobolev spaces with applications to Fourier transforms and spline interpolation, SIAM J. Numer. Anal., 7 (1970), 112–124.

Babu71 I. Babuska, Error bounds for finite element method, Numer. Math. 16 (1971) 322–333 .

StFi73 G. Strang, G.J. Fix, An Analysis of the Finite Element Method, Prentice-Hall Series in Automatic Computation (1973).


