Generalized Beta Prime Distribution Applied to Finite Element Error Approximation

Abstract

In this paper, we propose a new family of probability laws based on the Generalized Beta Prime distribution to evaluate the relative accuracy between two Lagrange finite elements $P_{k_1}$ and $P_{k_2},(k_1<k_2)$. Usually, the relative finite element accuracy is based on the comparison of the asymptotic speed of convergence, when the mesh size h goes to zero. The new probability laws we propose here highlight that there exists, depending on h, cases where the $P_{k_1}$ finite element is more likely accurate than the $P_{k_2}$ element. To confirm this assertion, we highlight, using numerical examples, the quality of the fit between the statistical frequencies and the corresponding probabilities, as determined by the probability law. This illustrates that, when h goes away from zero, a finite element $P_{k_1}$ may produce more precise results than a finite element $P_{k_2}$, since the probability of the event “$P_{k_1}$ is more accurate than $P_{k_2}$” becomes greater than 0.5. In these cases, finite element $P_{k_2}$ is more likely overqualified.

1. Introduction

Recently, we proposed a new vision to consider error estimates applied to finite element approximations (see [1,2,3]).

Indeed, we followed the emerging field called probabilistic numerics whose founding principles and aims consist of providing a statistical treatment of the errors and/or approximations that are made en route to the output of a deterministic numerical method, particularly the discretized solution of a partial differential equation (see, for example, [4,5,6]).

Then, each method that belongs to the field of probabilistic numerics introduces a probability measure on the approximated solution of a given classical numerical method.

One of the main applications of probabilistic numerics leads to the possibility for quantifying uncertainties contained in numerical approximations due to several causes, like numerical errors, randomness of the mesh generator, and so on.

In the same spirit, we introduced a probability measure of the classical error estimate considered as a random variable to derive two probability laws. This allowed us to model the relative accuracy between two Lagrange finite elements $P_{k_1}$ and $P_{k_2},(k_1<k_2),$ analyzed as a random variable.

This new point of view enabled us to complete the classical results, which are limited to the convergence analysis. Indeed, the relative accuracy between two finite elements $P_{k_1}$ and $P_{k_2},(k_1<k_2)$ is usually obtained by comparing the asymptotic rate of convergence when the mesh size h goes to zero. This enables one to conclude that the $P_{k_2}$ finite element is more accurate than $P_{k_1}$, since $h^{k_2}$ goes faster to zero than $h^{k_1}$.

The probability laws we derived confirmed our intuition (see also our previous approaches in [7]), that when h is set to a fixed value, greater than an undetermined threshold, the asymptotic behavior recalled above does not necessarily apply. Indeed, in the error estimates, the upper bound of the approximation error is an unknown constant that depends, among others, on a given semi-norm of the unknown exact solution. Hence, the numerical comparison of error approximations between $P_{k_1}$ and $P_{k_2}$ cannot be achieved.

As a consequence, determining the smallest of the two concerned approximation errors remains an open question.

Starting from there, we first consider a given approximation error as a positive number whose position is unknown within an interval defined by the unknown upper bound of the corresponding error estimate. Then, we consider this position as the result of a random variable. Indeed, the approximation error depends on the approximation, hence, on quantitative uncertainties which are related to the mesh generator process.

In addition to the new insights we obtained from these probability laws, we also implemented practical cases [8], which enabled us to evaluate the quality of the fit between these probability laws and the corresponding statistical frequencies. We show that these two probability laws globally behave like the corresponding statistical frequencies. However, the fit was sometimes not precise enough. Next, we identified the reasons for this lack of accuracy: it basically originates in the excessive rigidity of the probabilistic assumptions we made to derive these laws.

Thus, we developed a new generation of probabilistic models based on the Generalized Beta Prime distribution (GBP). These models enabled us to derive, under probabilistic acceptable hypotheses, the probability law of the relative accuracy between two Lagrange finite elements.

In this paper, we introduce the probabilistic framework built to obtain this law; we justify its use and show that it fits well with several examples. This confirms, in a significant manner, the relevance of considering error estimates as random variables in a suitable probabilistic environment.

The paper is organized as follows. In Section 2, we recall the mathematical problem we consider and a corollary of the Bramble-Hilbert lemma, from which the previous probabilistic laws [8] have been derived. In Section 3, we show a typical numerical result we obtained combining numerical statistics and these probability laws. Then, in Section 4, we derive the new probability law that enable one to evaluate the relative error accuracy between two finite elements $P_{k_1}$ and $P_{k_2},(k_1<k_2)$ relying on the GBP distribution. Finally, in Section 5, we show, using several examples, the high quality of the fit between the GBP probabilistic law and the corresponding statistical frequencies. Concluding remarks follow.

2. Abstract Problem and Finite Element Error Estimates

We consider an open-bounded and non-empty subset $\mathsf{\Omega }$ of ${\mathbb{R}}^n$, and denote by $\mathsf{\Gamma }$ its boundary, assumed to be $C^1-$ piecewise. We also introduce a Hilbert space V endowed with a norm ${∥.∥}_V$, and a bilinear, continuous and $V-$ elliptic form $a(·,·)$ defined on $V\times V$. Finally, $l(·)$ denotes a linear continuous form defined on V.

Let $u\in V$ be the unique solution to the second order elliptic variational formulation:

$$\left\{\begin{array}{c}\mathrm{Find} u\in V\mathrm{solution} \mathrm{to}:\hfill \\ a(u,v)=l\left(v\right),\forall v\in V.\hfill \end{array}\right.$$

(1)

In this paper, we restrict ourselves to the simple case where V is the usual Sobolev space of distributions $H^1\left(\mathsf{\Omega }\right)$. More general cases can be found in [2].

Now, for a given $h>0$, let us introduce $V_h$, a finite-dimensional subset of V, and consider $u_h\in V_h$ an approximation of u, the solution to the approximate variational formulation:

$$\left\{\begin{array}{c}\mathrm{Find} u_h\in V_h\mathrm{solution} \mathrm{to}:\hfill \\ a(u_h,v_h)=l\left(v_h\right),\forall v_h\in V_h.\hfill \end{array}\right.$$

(2)

In what follows, we are interested in evaluating error estimates for the finite element method. Hence, we first assume that domain $\mathsf{\Omega }$ is exactly covered by a regular ${\mathcal{T}}_h$ mesh composed by $N_s$ n-simplexes $K_j,(1\le j\le N_s),$ which respects the classical rules for regular discretization (see, for example, ref. [9] for the bidimensional case or [10,11] in ${\mathbb{R}}^n$). We also denote by $P_k\left(K_j\right)$ the space of polynomial functions defined on a given n-simplex $K_j$ of degree less than or equal to k, ($k\ge$ 1).

Let us first briefly recall the result of [11], which is the starting point of our study. Let ${\parallel .\parallel }_1$ be the classical norm in $H^1\left(\mathsf{\Omega }\right)$ and ${|.|}_{k+1}$ the semi-norm in $H^{k+1}\left(\mathsf{\Omega }\right)$, and h the mesh size, namely the largest diameter of the elements of mesh ${\mathcal{T}}_h$. We then have:

Lemma 1.

Suppose that there exists an integer $k\ge 1$ such that approximation $u_h$ of $V_h$ is a continuous piecewise function composed by polynomials that belong to $P_k\left(K_j\right),(1\le j\le N_s)$.

Then, if the exact solution u to (1) belongs to $H^{k+1}\left(\mathsf{\Omega }\right)$, we have the following error estimate:

$$\parallel u_h{-u\parallel }_1\le {\mathcal{C}}_k{h}^k{\left|u\right|}_{k+1},$$

(3)

where ${\mathcal{C}}_k$ is a positive constant independent of h.

Let us now consider two families of Lagrange finite elements $P_{k_1}$ and $P_{k_2}$ corresponding to a set of values $(k_1,k_2)\in {\mathbb{N}}^2$ such that $0<k_1<k_2$.

The two corresponding inequalities given by (3), assuming that the solution u to (1) belongs to $H^{k_2+1}\left(\mathsf{\Omega }\right)$, are:

$$\begin{array}{ccc}\hfill \parallel u_h^{\left(k_1\right)}{-u\parallel }_{1,\mathsf{\Omega }}& \le & {\mathcal{C}}_{k_1}{h}^{k_1}{\left|u\right|}_{k_1+1,\mathsf{\Omega }},\hfill \end{array}$$

(4)

$$\begin{array}{ccc}\hfill \parallel u_h^{\left(k_2\right)}{-u\parallel }_{1,\mathsf{\Omega }}& \le & {\mathcal{C}}_{k_2}{h}^{k_2}{\left|u\right|}_{k_2+1,\mathsf{\Omega }},\hfill \end{array}$$

(5)

where $u_h^{\left(k_1\right)}$ and $u_h^{\left(k_2\right)}$, respectively, denote the $P_{k_1}$ and $P_{k_2}$ Lagrange finite element approximations of u.

Now, if we consider a given mesh for the finite element of $P_{k_2}$, which contains that of $P_{k_1}$, then, for the particular class of problems where the variational formulation (1) is equivalent to a minimization formulation (see for example [9]), one can show that the approximation error of $P_{k_2}$ is always lower than that of $P_{k_1}$, and that $P_{k_2}$ is more accurate than $P_{k_1}$ for all values of mesh size h.

Next, for a given mesh size value of h, we consider two independent meshes for $P_{k_1}$ and $P_{k_2}$ built by a mesh generator (see Figure 1).

Usually, in order to compare the relative accuracy between these two finite elements, one asymptotically considers inequalities (4) and (5) and concludes that, when h goes to zero, $P_{k_2}$ is more accurate than $P_{k_1}$, since $h^{k_2}$ goes faster to zero than $h^{k_1}$.

However, for a given application, h has a fixed value and this method of comparison is no more valid. Therefore, we propose a solution to determine the relative accuracy between the finite elements $P_{k_1}$ and $P_{k_2},(k_1<k_2)$, for a given value of h, for which two independent meshes are considered.

To this end, let us set:

$${\beta }_{k_1}={\mathcal{C}}_{k_1}{h}^{k_1}{\left|u\right|}_{k_1+1,\mathsf{\Omega }}\mathrm{and}{\beta }_{k_2}={\mathcal{C}}_{k_2}{h}^{k_2}{\left|u\right|}_{k_2+1,\mathsf{\Omega }},$$

(6)

that allows us to rewrite (4) and (5) as:

$$\begin{array}{ccc}\hfill \parallel u_h^{\left(k_1\right)}{-u\parallel }_{1,\mathsf{\Omega }}& \le & {\beta }_{k_1},\hfill \end{array}$$

(7)

$$\begin{array}{ccc}\hfill \parallel u_h^{\left(k_2\right)}{-u\parallel }_{1,\mathsf{\Omega }}& \le & {\beta }_{k_2}.\hfill \end{array}$$

(8)

Now, as explained in [3], there is no a priori available information to surely or better specify the relative position of the approximation errors $\parallel u_h^{\left(k_1\right)}{-u\parallel }_{1,\mathsf{\Omega }}$ and $\parallel u_h^{\left(k_2\right)}{-u\parallel }_{1,\mathsf{\Omega }}$, which respectively belong to the interval $[0,{\beta }_{k_1}]$ and $[0,{\beta }_{k_2}]$, where ${\beta }_{k_i}$ are unknown since ${\mathcal{C}}_{k_i}$ and ${\left|u\right|}_{k_i+1,\mathsf{\Omega }}$ are unknown quantities.

In [3], we explained why quantitative uncertainties have to be handled using finite element computations. This originates in the way the mesh grid generator builds the mesh used to compute the approximation $u_h^{\left(k_i\right)},(i=1,2)$; this leads to a partial non-control of the mesh, even for a given maximum mesh size. Indeed, most meshes are generated by a mesh generator based on a Delaunay-Voronoi method. As a consequence, the corresponding point generation, and/or the treatment of ambiguities inevitably contains some randomness. Therefore, the generated grid is a priori random, and so is the corresponding approximation $u_h^{\left(k_i\right)}(i=1,2)$.

For these reasons, we introduced in [3] a probabilistic framework that enables one to consider the possible values of the norm $\parallel u_h^{\left(k\right)}{-u\parallel }_{1,\mathsf{\Omega }}$ as a random variable, which is defined as follows:

• For a fixed value of the mesh size h, a random trial corresponds to the grid generation and the associated approximation $u_h^{\left(k\right)}$.
• The probability space $\mathsf{\Omega }$, therefore, contains all the possible results for a given random trial, namely, all of the possible grids that the mesh generator may construct, or equivalently, all of the corresponding associated approximations $u_h^{\left(k\right)}$.

Then, for a fixed value of k, we define the random variable $X^{\left(k\right)}$ as follows:

$$\begin{array}{ccc}\hfill X^{\left(k\right)}:& \mathsf{\Omega }& \to [0,{\beta }_k]\hfill \end{array}$$

(9)

$$\begin{array}{ccc}& \omega \equiv u_h^{\left(k\right)}& {\displaystyle \mapsto X^{\left(k\right)}\left(\omega \right)=X^{\left(k\right)}\left(u_h^{\left(k\right)}\right)={\parallel u_h^{\left(k\right)}-u\parallel }_{1,\mathsf{\Omega }}.}\hfill \end{array}$$

(10)

In what follows, for simplicity, we will set: $X^{\left(k\right)}\left(u_h^{\left(k\right)}\right)\equiv X^{\left(k\right)}\left(h\right)$.

Therefore, our goal is now to evaluate the probability of the event:

$$\left\{\parallel u_h^{\left(k_2\right)}{-u\parallel }_{1,\mathsf{\Omega }}\le {\parallel u_h^{\left(k_1\right)}-u\parallel }_{1,\mathsf{\Omega }}\right\}\equiv \left\{X^{\left(k_2\right)}\left(h\right)\le X^{\left(k_1\right)}\left(h\right)\right\},$$

(11)

which will enable us to estimate the more likely accuracy of the two finite elements $P_{k_1}$ and $P_{k_2}$, $(k_1<k_2)$.

Now, regarding the absence of information concerning the more or less likely values of the norm $\parallel u_h^{\left(k_i\right)}{-u\parallel }_{1,\mathsf{\Omega }},(i=1,2)$ in the interval $[0,{\beta }_{k_i}],(i=1,2)$, we assumed in [1], that each random variable $X^{\left(k_i\right)},(i=1,2)$ had a uniform distribution on the interval $[0,{\beta }_{k_i}]$, and that they were also independent.

This probabilistic framework enabled us to obtain, in a more general context (see Theorem 3.1 in [1]), the density of the probability of the random variable Z defined by $Z=X^{\left(k_2\right)}-X^{\left(k_1\right)}$ and, consequently, ${\mathcal{P}}_{k_1,k_2}\equiv Prob\left\{X^{\left(k_2\right)}\le X^{\left(k_1\right)}\right\}$. This probability corresponds to the value of the entire cumulative distribution function $F_Z\left(z\right)$ at $z=0$, defined by:

$$F_Z\left(z\right)={\int }_{-\infty }^z{f}_Z\left(z\right)dz,$$

(12)

with $f_Z\left(z\right)$ denoting the density function of the random variable Z.

For the purpose of the present work, by performing elementary adaptations of Theorem 3.3 in [1], we get the corresponding probability law given by:

$${\displaystyle {\mathcal{P}}_{k_1,k_2}\left(h\right)=|\begin{array}{}\left(13\right)& 1-\frac{1}{2}{\left(\frac{h}{h_{k_1,k_2}^{*}}\right)}^{k_2-k_1}\mathrm{if}0\le h\le h_{k_1,k_2}^{*},\left(14\right)& \frac{1}{2}{\left(\frac{h_{k_1,k_2}^{*}}{h}\right)}^{k_2-k_1}\mathrm{if}h\ge h_{k_1,k_2}^{*},\end{array}}$$

where $h_{k_1,k_2}^{*}$ is defined by:

$${\displaystyle h_{k_1,k_2}^{*}\equiv {\left(\frac{{\mathcal{C}}_{k_1}{\left|u\right|}_{k_1+1,\mathsf{\Omega }}}{{\mathcal{C}}_{k_2}{\left|u\right|}_{k_2+1,\mathsf{\Omega }}}\right)}^{\frac{1}{k_2-k_1}}.}$$

(15)

The shape of this law is similar to a “sigmoid” curve, as one can see in Figure 2.

We already remarked in [3] that the probability law (13) and (14) can asymptotically—when $k_2-k_1$ goes to infinity—lead to the limit situation we defined as the “two-steps” model (see Figure 2).

Nevertheless, we also proved in [3] that under suitable probabilistic assumptions, one can directly derive this “two-steps” probabilistic law for all non-zero integers $k_1$ and $k_2$ to finally obtain the following probability law:

$${\displaystyle {\mathcal{P}}_{k_1,k_2}\left(h\right)=\left|\begin{array}{cc}1\hfill & \mathrm{if}0<h<h_{k_1,k_2}^{*},\hfill \\ 0\hfill & \mathrm{if}h>h_{k_1,k_2}^{*}.\hfill \end{array}\right.}$$

(16)

The next section is dedicated to the analysis of the fit between statistical data and the above probability laws (13), (14) and (16).

3. Numerical Statistics and Probability-Law Comparison

Here, we compare the theoretical probability law (13) and (14) and the statistics obtained from two examples of approximate variational formulations (2).

As shown in [8], the “two-steps” law (16) provides good-quality fits in several numerical cases. However, this law is a bit rough and cannot really follow the variations of the convexity of the statistical data curve (see below). For this reason, we also tested the accuracy of the fit with the “sigmoid” law given by (13) and (14).

In order to numerically check the accuracy of the fit between the probability law (13) and (15) and the corresponding statistical data produced by numerical simulations, we considered two finite elements $P_{k_1}$ and $P_{k_2}$ ($k_1<k_2$).

One of the main ingredients of the method is the computation of $h_{k_1,k_2}^{*}$, as defined by (15), which is evaluated using a maximum likelihood estimator; see for instance [12].

In our case, this principle is applied as follows: for a given finite element $P_{k_1}$, we consider N different meshes with the same (maximum) mesh size h. Then, we compute:

$${\displaystyle \underset{N,h}{max}\frac{\parallel u_h^{\left(k_1\right)}{-u\parallel }_1}{h^{k_1}}}$$

(17)

which constitutes, using estimate (4), the maximum likelihood estimator for ${\mathcal{C}}_{k_1}{\left|u\right|}_{k_1+1}$. Indeed, due to inequality (4), quantity $\frac{X^{k_1}\left(h\right)}{h^{k_1}}$ is also a uniform random variable whose support is $[0,{\mathcal{C}}_{k_1}|u{|}_{k_1+1}]$.

Then, one can show [13] that for a given uniform random variable Y whose support is $[0,\theta ]$, $\theta$ being an unknown real parameter, the maximum likelihood estimator $\widehat{\theta }$ is given by:

$$\widehat{\theta }=max(Y_1,\dots ,Y_N),$$

where $(Y_1,\dots ,Y_N)$ is a sample built with independent and identically distributed random variables ${\left(Y_i\right)}_{i=1,N}$, with the same distribution as Y.

In our case, this implies that (17) is the maximum likelihood estimator for $\mathcal{C}{k}_1{\left|u\right|}_{k_1+1}$, since N and h each take a finite number of values.

Doing the same for another finite element $P_{k_2},(k_1<k_2)$, we obtain that the estimator for $h_{k_1,k_2}^{*}$, denoted $\widehat{h_{k_1,k_2}^{*}}$, is defined by:

$$\widehat{h_{k_1,k_2}^{*}}={\left({\displaystyle \frac{{\displaystyle \underset{N,h}{max}\frac{\parallel u_h^{\left(k_1\right)}{-u\parallel }_1}{h^{k_1}}}}{{\displaystyle \underset{N,h}{max}\frac{\parallel u_h^{\left(k_2\right)}{-u\parallel }_1}{h^{k_2}}}}}\right)}^{1/k_2-k_1}$$

(18)

Then, one can easily compute the probability law given by (13) and (14), when one replaces $h_{k_1,k_2}^{*}$ by its estimator $\widehat{h_{k_1,k_2}^{*}}$.

Finally, in order to numerically check the validity of this probability law, we now compare the probabilities computed by (13), (14) and (18), with the corresponding statistical frequencies computed on the N meshes, for each fixed value h of the mesh size, as follows.

We considered for the two finite elements $P_{k_1}$ and $P_{k_2}$ ($k_1<k_2$), the same number N of different meshes with the same (maximum) mesh size h. From there, we computed the approximate solution $u_h^{\left(k_1\right)}$ and $u_h^{\left(k_2\right)}$, and tested if $\parallel u_h^{\left(k_2\right)}{-u\parallel }_1\le {\parallel u_h^{\left(k_1\right)}-u\parallel }_1$. Then, we repeated the same process for different values of h, either lower or greater than $\widehat{h_{k_1,k_2}^{*}}$. This gives, as a function of h, the percentage of cases where the approximation error of $P_{k_2}$ is lower than the approximation error of $P_{k_1}$.

In all cases, we used package FreeFem++ [14] to compute the $P_{k_1}$ and $P_{k_2}$ finite element approximations. Basically, the mesh generator of FreeFem++ is based on the Delaunay-Voronoi algorithm, and allows us to fix the maximum mesh size h for N different meshes of the same domain $\mathsf{\Omega }$. Since FreeFEM++ is based on updated and modern methods of meshing, we are confident that using another mesh generator will lead to equivalent results.

Now, in order to motivate the new probability law of the next section, we recall typical results obtained in [8]. To this end, let us look at the numerical tests we implemented to analyze the relative accuracy between two finite elements.

Hence, we consider the following classical elliptic problem, with obvious notations: Find $u\in H^1\left(\mathsf{\Omega }\right)$ solution to:

$$\left\{\begin{array}{c}-\mathsf{\Delta }u=f \mathrm{in} \mathsf{\Omega },\\ u=g \mathrm{on} \partial \mathsf{\Omega },\end{array}\right.$$

(19)

where, for simplicity, $\mathsf{\Omega }$ is the open unit square in ${\mathbb{R}}^2$: $\mathsf{\Omega }=]0,1[\times ]0,1[$. The associated variational formulation, which is analogous to (1), can readily be derived. According to the choice of f and g, we considered, as examples, a stiff problem, where the solution exhibits rapid variations or, alternatively, a very smooth problem.

More precisely, when the right-hand side $f(x,y)$ is the Runge-like function, (see [15] or [16]), given by:

$${\displaystyle f(x,y)=\frac{1}{1+\alpha x^2}\frac{1}{1+\alpha y^2},}$$

$\alpha$ being a real parameter, we derived in [8], the exact solution of (19) according to a suitable choice of the Dirichlet condition, parametered by g.

Then, considering 100 meshes for each value of h and $\alpha =3000$ in the Runge function, we plot in Figure 3 the statistical frequencies corresponding to “$P_3$ more accurate that $P_1$” (blue line), together with the corresponding “sigmoid” probability law (orange line).

As one can see, the fit between the statistics and the probability law is not satisfactory. The same gaps were observed for other simulations with different sets of parameter values, and for several pairs of finite elements $P_{k_1}$ and $P_{k_2}$.

We now illustrate the probabilistic law for a very smooth solution to the variational problem (1). To build such a case, we choose the right-hand side f of problem (19) as:

$f(x,y)=2{\pi }^{2}sin\left(\pi x\right)cos\left(\pi y\right)$, so that $u(x,y)=sin\left(\pi x\right)cos\left(\pi y\right)$ is the exact solution of the problem, provided that the Dirichlet boundary condition g is taken as the trace of $u(x,y)$ on the boundary $\partial \mathsf{\Omega }$ (for more details, see [8]).

As previously, we first compute $\widehat{h_{k_1,k_2}^{*}}$ defined by (15), then we compute the probabilistic model recalled above. After that, we compare these results to the statistical frequencies.

As an example, we consider the finite elements $P_2$ and $P_3$, the corresponding results being depicted in Figure 4.

As one can see, the comparison with the sigmoid law in (13) and (14) gives only a global trend and is not accurate.

Here, as for the previous Runge example, this lack of accuracy will be improved by the new probabilistic laws derived in the next section. We will show how to enrich the “sigmoid” model, which actually only depends on one parameter, namely, $h_{k_1,k_2}^{*}$.

4. A New Probability Law Based on GBP

This section is devoted to the new probabilistic law we derived in order to evaluate the relative error accuracy between two finite elements $P_{k_1}$ and $P_{k_2},(k_1<k_2)$. To motivate this new approach, let us make some remarks:

1.

In the previous works that enabled us to derive the “Sigmoid” probability law, we had made some assumptions (i.e., uniform distribution and independency). Since we now aim at obtaining a better fit between the probabilistic law and the statistical data, we will now relax the hypothesis of uniformity we had made on the densities of the random variables $X^{\left(k_i\right)}\left(h\right),(i=1,2)$, (see [1,2,3]).

2.

In order to obtain more degrees of freedom for the shape of these densities, let us first consider the random variable $Z=X^{\left(k_2\right)}-X^{\left(k_1\right)}$. Our goal is to obtain, for the cumulative distribution function $F_Z$ defined by (12) at point $z=0$, a curve whose shape looks like a “Sigmoid”. We thus have to enrich our modeling process by adding more degrees of freedom, the “Sigmoid” probability law (13) and (14) containing only one parameter $h_{k_1,k_2}^{*}$. For this reason, we now consider a density $f_Z$ for the random variable Z, whose associated cumulative distribution function $F_Z$ at point $z=0$ will include two exogenous parameters; these will be statistically estimated.

Keeping in mind these remarks, we begin by introducing the probability density function $f_X$ of the normalized Beta random variable X (see [17]), defined by:

$${\displaystyle \forall x\in [0,1]:f_X(x;p,q)\equiv \frac{x^{p-1}{(1-x)}^{q-1}}{{\displaystyle {\int }_0^1{u}^{p-1}{(1-u)}^{q-1}du}}=\frac{\mathsf{\Gamma }(p+q)}{\mathsf{\Gamma }\left(p\right)\mathsf{\Gamma }\left(q\right)}{x}^{p-1}{(1-x)}^{q-1},}$$

(20)

with p and q, two parameters of shape that belong to ${\mathbb{R}}_{+}^{*}$, and $\mathsf{\Gamma }(.)$, the classical Gamma function.

Among the numerous properties of the Beta distribution, let us mention the one that particularly motivates us into considering it. Depending on the two parameters p and q, the shapes of the corresponding cumulative distributions are very rich. In particular, they include the shape of the “Sigmoid” curve we are looking for to fit the statistics we get in this numerical example.

However, one cannot directly apply the Beta distribution to get the probability law we are looking for. More specifically, two properties have to be taken into account:

1.

The support of the Beta density $f_X$ of the random variable X, denoted $SuppX$, is included in $[0,1]$. Nevertheless, that of the random variable Z is $[-{\beta }_{k_1},{\beta }_{k_2}]$, since $Z=X^{\left(k_2\right)}-X^{\left(k_1\right)}$ and $Supp{X}^{\left(k_i\right)}\subset [0,{\beta }_{k_i}],(i=1,2)$. This will make us use a suitable transformation of the density $f_X$ to guarantee the correct support of the density $f_Z$ of Z.

2.

We are looking for a probability law for the event ${\displaystyle \left\{X^{\left(k_2\right)}\le X^{\left(k_1\right)}\right\}}$ as a function of h, that belongs to $[0,+\infty [$. Consequently, we have to apply another transformation to the density $f_Z$ to guarantee this property for the support of h.

To achieve these transformations, we establish the following results:

Lemma 2.

Let Z be the random variable defined by $Z\equiv X^{\left(k_2\right)}-X^{\left(k_1\right)}$, where $X^{\left(k_i\right)}(i=1,2)$ are defined by (10). Let $X\hookrightarrow B(p,q)$ be the Beta distribution parameterized by two given shape parameters $(p,q)\in {\mathbb{R}}_{+}^{*2}$. Assume that Z is given by:

$$Z=-{\beta }_{k_1}+({\beta }_{k_1}+{\beta }_{k_2})X.$$

(21)

Then, the probability density function $f_Z$ of Z can be written as:

$${\displaystyle f_Z\left(z\right)=\frac{\mathsf{\Gamma }(p+q)}{\mathsf{\Gamma }\left(p\right)\mathsf{\Gamma }\left(q\right)}\frac{{\beta }_{k_1}^{p-1}{\beta }_{k_2}^{q-1}}{{({\beta }_{k_1}+{\beta }_{k_2})}^{p+q-1}}{\left(1+\frac{z}{{\beta }_{k_1}}\right)}^{p-1}{\left(1-\frac{z}{{\beta }_{k_2}}\right)}^{q-1}{𝟙}_{[-{\beta }_{k_1},{\beta }_{k_2}]}\left(z\right),}$$

(22)

where ${𝟙}_{[-{\beta }_{k_1},{\beta }_{k_2}]}$ is the indicator function of the interval $[-{\beta }_{k_1},{\beta }_{k_2}]$.

Proof. 

As we already noticed, the support of $f_Z$ is clearly within $[-{\beta }_{k_1},{\beta }_{k_2}]$ as soon as those of $X^{\left(k_i\right)}$ are in $[0,{\beta }_{k_i}]$.

Now, let us evaluate the cumulative distribution function $F_Z\left(z\right)$ of the random variable Z as defined by (21):

$$\begin{array}{ccc}\hfill {\displaystyle F_Z\left(z\right)}& =& {\displaystyle Prob\left\{Z\le z\right\}=Prob\left\{-{\beta }_{k_1}+({\beta }_{k_1}+{\beta }_{k_2})X\le z\right\}}\hfill \end{array}$$

(23)

$$\begin{array}{ccc}& =& Prob\left\{{\displaystyle X\le \frac{z+{\beta }_{k_1}}{{\beta }_{k_1}+{\beta }_{k_2}}}\right\}=\frac{\mathsf{\Gamma }(p+q)}{\mathsf{\Gamma }\left(p\right)\mathsf{\Gamma }\left(q\right)}{\int }_0^{\frac{z+{\beta }_{k_1}}{{\beta }_{k_1}+{\beta }_{k_2}}}{u}^{p-1}{(1-u)}^{q-1}du.\hfill \end{array}$$

(24)

Then, by deriving $F_Z\left(z\right)$ as given by (24), we obtain the expression (22) of the density $f_Z$. □

Remark 1.

Let us now give an interpretation of this result. As mentioned above, in our previous works (see, for example, [1,3]), we assumed the random variables $X^{\left(k_i\right)}$ to be uniformly distributed on their support $[0,{\beta }_{k_i}]$. This enabled us to derive the density $f_Z$ of the random variable Z (see Theorem 3.1 in [1]).

Here, Lemma 2 can be interpreted regarding the new assumption implicitly made on the random variables $X^{\left(k_i\right)}$. Indeed, if we rewrite Z as:

$$Z={\beta }_{k_2}X-({\beta }_{k_1}-{\beta }_{k_1}X),$$

by setting:

$$X^{\left(k_1\right)}\equiv {\beta }_{k_1}-{\beta }_{k_1}Xand{X}^{\left(k_2\right)}={\beta }_{k_2}X,$$

(25)

we observe that the support of $X^{\left(k_i\right)}$ belongs to $[0,{\beta }_{k_i}]$, on the one hand, and that the difference $X^{\left(k_2\right)}-X^{\left(k_1\right)}$ is equal to Z, on the other hand.

In other words, the choice to write Z given by (21) as a dimensional Beta distribution on $[-{\beta }_{k_1},{\beta }_{k_2}]$ (which corresponds to altering the location and scale of the standard Beta distribution), accordingly modifies the hypothesis of uniformity of the two variables $X^{\left(k_i\right)}$ to the one of a non-dimensional Beta distribution.

We are now capable to derive the new probability distribution ${\mathcal{P}}_{k_1,k_2}\left(h\right)$ in order to evaluate which, between two Lagrange finite elements $P_{k_1}$ and $P_{k_2},(k_1<k_2)$, is the more likely accurate. This is the purpose of the following Theorem.

Theorem 1.

Let $X^{\left(k_i\right)}(i=1,2),$ be the two random non-dimensional Beta distribution variables defined by (25) and Z the corresponding random variable defined by (21).

Then, ${\mathcal{P}}_{k_1,k_2}\left(h\right)$ is the cumulative distribution function of a Generalized Beta Prime random variable H, whose density of probability $f_H$, is defined by four parameters $(p,q,k_2-k_1,h_{k_1,k_2}^{*})$, and we have:

$${\displaystyle {\mathcal{P}}_{k_1,k_2}\left(h\right)=Prob\{H\ge h\}={\int }_h^{+\infty }{f}_H(s;q,p,k_2-k_1,h_{k_1,k_2}^{*})ds,}$$

(26)

where:

$${\displaystyle f_H(s;q,p,k_2-k_1,h_{k_1,k_2}^{*})=\frac{\mathsf{\Gamma }(p+q)}{\mathsf{\Gamma }\left(p\right)\mathsf{\Gamma }\left(q\right)}\frac{(k_2-k_1)}{h_{k_1,k_2}^{*}}{\left(\frac{s}{h_{k_1,k_2}^{*}}\right)}^{q(k_2-k_1)-1}{\left[1+{\left(\frac{s}{h_{k_1,k_2}^{*}}\right)}^{k_2-k_1}\right]}^{-p-q}.}$$

(27)

Proof. 

On the one hand, we are looking for a random variable H whose $[0,+\infty [$, which is associated to ${\mathcal{P}}_{k_1,k_2}\left(h\right)$ equal to $F_Z\left(0\right)$. On the other hand, (24) was derived by the help of the dimensionless Beta distribution $B(p,q)$, whose support belongs to $[0,1]$. Hence, we set the following change of variable in (24):

$${\displaystyle s=\frac{u}{1-u},}$$

and we get:

$${\displaystyle F_Z\left(z\right)=\frac{\mathsf{\Gamma }(p+q)}{\mathsf{\Gamma }\left(p\right)\mathsf{\Gamma }\left(q\right)}{\int }_0^{\frac{z+{\beta }_{k_1}}{{\beta }_{k_2}-z}}{s}^{p-1}{(1+s)}^{-p-q}ds.}$$

(28)

Now, taking $z=0$ in (28), and using ${\beta }_{k_i}(i=1,2)$, given by (6) and $h_{k_1,k_2}^{*}$ by (15), we obtain:

$$\begin{array}{ccc}\hfill {\displaystyle F_Z\left(0\right)={\mathcal{P}}_{k_1,k_2}\left(h\right)}& =& {\displaystyle \frac{\mathsf{\Gamma }(p+q)}{\mathsf{\Gamma }\left(p\right)\mathsf{\Gamma }\left(q\right)}{\int }_0^{\frac{{\beta }_{k_1}}{{\beta }_{k_2}}}{s}^{p-1}{(1+s)}^{-p-q}ds,}\hfill \end{array}$$

(29)

$$\begin{array}{ccc}& =& \frac{\mathsf{\Gamma }(p+q)}{\mathsf{\Gamma }\left(p\right)\mathsf{\Gamma }\left(q\right)}{\int }_0^{{(h_{k_1,k_2}^{*}/h)}^{k_2-k_1}}{s}^{p-1}{(1+s)}^{-p-q}ds.\hfill \end{array}$$

(30)

A last change of variables, ${\displaystyle t=\frac{1}{s}}$, in (30) leads to:

$${\displaystyle F_Z\left(0\right)=\frac{\mathsf{\Gamma }(p+q)}{\mathsf{\Gamma }\left(p\right)\mathsf{\Gamma }\left(q\right)}{\int }_{{(h/h_{k_1,k_2}^{*})}^{k_2-k_1}}^{+\infty }{t}^{q-1}{(1+t)}^{-p-q}dt.}$$

(31)

Finally, by considering the probability of the complementary event considered in (31) as a function of h (namely, $1-{\mathcal{P}}_{k_1,k_2}\left(h\right)$), and after derivation with respect to h, we get the density $f_H$ defined in (27), which is a GBP probability density function parameterized by $(q,p,k_2-k_1,h_{k_1,k_2}^{*})$, see for example [17].

In other words, we have obtained:

$$Prob\left\{X^{\left(k_1\right)}\le X^{\left(k_2\right)}\right\}=1-Prob\left\{X^{\left(k_2\right)}\le X^{\left(k_1\right)}\right\}\left(h\right)={\int }_o^h{f}_H(s;q,p,k_2-k_1,h_{k_1,k_2}^{*})ds,$$

(32)

and ${\displaystyle {\mathcal{P}}_{k_1,k_2}\left(h\right)}$ in (26) is deduced by complementarity. □

Remark 2.

1.

From (26) and (27), or equivalently from (30), we can obtain the asymptotic behavior of the probability of the event ${\displaystyle \left\{X^{\left(k_1\right)}\le X^{\left(k_2\right)}\right\}}$ as the mesh size h goes to 0.

Clearly, $Prob\left\{X^{\left(k_1\right)}\le X^{\left(k_2\right)}\right\}$ goes to 0 with h. In other words, we have found with the new probabilistic law the usual result, which claims that the finite element $P_{k_2}$ is more accurate than $P_{k_1}$, since from (4) and (5), $h^{k_2}$ goes faster to zero than $h^{k_1}$ when $k_1<k_2$.

Here, by (30), the same property is expressed in terms of probability, namely: the event ${\displaystyle \left\{X^{\left(k_1\right)}\le X^{\left(k_2\right)}\right\}}$ occurs almost never, or equivalently, the event ${\displaystyle \left\{X^{\left(k_2\right)}\le X^{\left(k_1\right)}\right\}}$ is an almost surely one, since its probability is equal to 1.

2.

From (26), using the positivity of the density $f_H$, we conclude that ${\displaystyle {\mathcal{P}}_{k_1,k_2}\left(h\right)}$ is a decreasing function of h. This property was already observed with the “Sigmoid” model in [1,3] where uniformity of the random variables $X^{\left(k_i\right)},(i=1,2),$ was assumed.

In other words, this property claims that the event ${\displaystyle \left\{P_{k_2} is more accurate than P_{k_1}\right\}}$ is true only for small h, and the more h increases the less likely this property is. Moreover, when h becomes large, asymptotically the event ${\displaystyle \left\{P_{k_2} is more accurate than P_{k_1}\right\}}$ occurs almost never.

5. Numerical Comparisons

This section is devoted to the numerical comparison of the GBP probabilistic law (26), as derived in Theorem 1, with the statistical frequencies, computed as explained in Section 3.

Concerning the mesh generation, as explained above, we have considered N different meshes with the same (maximum) mesh size h. In these conditions, if we denote by ${\mathcal{T}}_h$ one of the N meshes, each other mesh is constructed as a random perturbation of ${\mathcal{T}}_h$ obtained by moving and/or adding internal vertices, in such a way that h remains constant, the number of points on the boundary being constant.

Let us also explain how we proceeded to determine the four parameters $(q,p,k_2-k_1,h_{k_1,k_2}^{*})$ of density $f_H$. First of all, we considered and fixed two finite elements, for example $P_1$ and $P_3$, which corresponds to $k_2-k_1=2$.

Regarding the three other parameters, $p,q$ and $h_{k_1,k_2}^{*}$, we have implemented a least squares optimization (see for instance [18]) to minimize the norm $L^2$ of the difference between the statistical frequencies and the corresponding values, computed by the GBP model. This was computed using the Excel solver.

The results we obtained are depicted in Figure 5 for $N=500$ different meshes. As one can see, the least squares algorithm found the optimal parameters p, q and $h_{k_1,k_2}^{*}$ such that the fit between the statistical frequencies and the GBP law (26) is very satisfactory.

This is clearly due to the the richness and flexibility of this distribution, which provides two supplementary degrees of freedom, p and q, with respect to the “Sigmoid” law. This also validates a posteriori the approach presented in Section 4, which describes the randomness values of the approximation errors $\parallel u_h^{\left(k_i\right)}{-u\parallel }_{1,\mathsf{\Omega }},(i=1,2),$ within their respective interval $[0,{\beta }_{k_i}],(i=1,2)$.

Another example that enables one to appreciate the accuracy of the fit may be obtained by comparing the results implemented with the two Lagrange finite elements $P_1$ and $P_4$, also for $N=500$. In this case, as one can see in Figure 6, the “Sigmoid” law only describes the global trend of the statistical frequencies, whereas the GBP law perfectly fits the corresponding data.

As a last illustration, let us present the results obtained when comparing the $P_2$ and $P_3$ finite elements. The difference in the quality of the fit between the “Sigmoid” law and the GBP law becomes even more important than in the previous case. Indeed, in this case, $k_2-k_1$ is equal to one and the “Sigmoid” law becomes a linear function of h when $h\le h_{k_1,k_2}^{*}$ in formula (13). In addition, the GBP law once again fit the statistical frequencies very well (see Figure 7).

6. Conclusions

In this paper, we derived a new family of probabilistic laws to compare the accuracy between two Lagrange finite elements $P_{k_1}$ and $P_{k_2}(k_1<k_2)$. Based on a new probabilistic approach introduced in [1], we extended these previous results to improve the fit between the statistical frequencies obtained by implementing practical cases and the corresponding probabilities.

We recall, in Section 3, the main problems we were confronted with regarding the accuracy of the previous probabilistic laws we derived in [1,2,3], by identifying the gaps we observed with the corresponding statistical frequencies.

Afterwards, in Section 4 we derived the new probabilistic law based on the Generalized Beta Prime law. In Section 5, we showed the very good fit we obtained between the numerical results and the new probability law. This corrected the main observed deficiencies described above with simpler laws.

Finally, the new probabilistic law, together with our statistical validation, confirmed the validity of considering the approximation errors as random variables defined in a well-adapted probabilistic framework. Of course, this probabilistic approach is not restricted to Lagrange $P_k$ finite elements and can be extended to any other finite elements, like rectangular $Q_k$. It can also be generalized to any types of approximation in the following sense: given a class of numerical schemes and their corresponding approximation error estimates, we should be able to order them not only in terms of asymptotic rate of convergence but also by evaluating the most probably accurate one.