A Probabilistic Approach for Solutions of deterministic PDEs as Well as Their Finite Element Approximations

: A probabilistic approach is developed for the exact solution u to a deterministic partial differential equation as well as for its associated approximation u ( k ) h performed by P k Lagrange ﬁnite element. Two limitations motivated our approach: On the one hand, the inability to determine the exact solution u relative to a given partial differential equation (which initially motivates one to approximating it) and, on the other hand, the existence of uncertainties associated with the numerical approximation u ( k ) h . We, thus, ﬁll this knowledge gap by considering the exact solution u together with its corresponding approximation u ( k ) h as random variables. By a method of consequence, any function where u and u ( k ) h are involved are modeled as random variables as well. In this paper, we focus our analysis on a variational formulation deﬁned on W m , p Sobolev spaces and the corresponding a priori estimates of the exact solution u and its approximation u ( k ) h in order to consider their respective W m , p -norm as a random variable, as well as the W m , p approximation error with regards to P k ﬁnite elements. This will enable us to derive a new probability distribution to evaluate the relative accuracy between two Lagrange ﬁnite elements P k 1 and P k 2 , ( k 1 < k 2 ) . 65N75


Introduction
We recently proposed new perspectives on relative finite element accuracy [1,2] by using a mixed geometrical-probabilistic interpretation of the error estimate in the case of finite element approximation (see for example [3] or [4]) derived from the Bramble-Hilbert lemma [5].
In [6], we further extended the results we had derived in the case of H 1 to Sobolev spaces W m,p , (p = 2). To this end, we had to consider a more general framework, which mainly relied on the Banach-Nečas-Babusȟka (BNB) abstract problem [7] devoted to Banach spaces.
This enabled us to obtain two new probability distributions that estimate the relative accuracy between two Lagrange finite elements P k 1 and P k 2 , (k 1 < k 2 ), by considering it as a random variable.
Thus, we obtained new results that show, among others, which of P k 1 or P k 2 is the most likely accurate, depending on the value of the mesh size h; this value is not considered anymore as proceeding towards zero, as in the standard procedure.
However, while obtaining these probability distributions, we only considered the standard error estimate dedicated to P k Lagrange finite elements, which approximates solution u to a variational problem (BNB), u being formulated in the present case in Sobolev space W m,p (Ω).
In the current study, we enrich the model published in [6] in multiple manners. Indeed, considering the functional framework (BNB) in the case of the W m,p Sobolev spaces, we take into account the available a priori estimates one can deal with with respect to the solution of this kind of problem as well as for its approximation, together with the approximation error which corresponds to P k Lagrange finite elements (see Section 2.2). This will enable us to derive a new probabilistic model applied to the relative accuracy between two finite elements P k 1 and P k 2 , (k 1 < k 2 ). To this end, we also generalize the discrete probabilistic framework we considered in [1,6] by introducing continuous probabilistic formalism based on an appropriate density of probability (see Section 3.1).
The paper is organized as follows. In Section 2, we recall the mathematical problem we considered and introduce the basic definitions of functional tools to consider different estimations in W m,p Sobolev spaces. Section 3.1 is dedicated to the analysis of the relative finite elements accuracy based on a probabilistic approach. In Section 4, we detail the contribution of the a priori estimates and their interactions with the error estimate in the probability distributions we derived in Section 3.1. Concluding remarks follow.

Abstract Problem in Banach Spaces and Corresponding Fundamental Results
In this section, we define an abstract framework that will enable us to consider the solution of a variational problem in Banach W m,p Sobolev spaces when p = 2 and its corresponding approximation computed by P k Lagrange finite elements.
In order to conduct this, we follow the presentation of A. Ern and J. L. Guermond [7], where two general Banach spaces W and V are involved where V is reflexive. We also recall the different assumptions needed in order to apply the (BNB) Theorem valid in Banach spaces.
Let u ∈ W be the solution of the variational formulation (VP) defined by the following: where the following is the case: 1.
W and V are two Banach spaces equipped with norms denoted by . W and . V , respectively; moreover, V is reflexive; 2.
a is a continuous bilinear form on W × V, i.e, a ∈ L(W × V; R) : l is a continuous linear form on V, i.e, l ∈ V = L(V; R).
We further make the two following assumptions.
Then, one can prove the (BNB) Theorem ( [7], Theorem 2.6), which claims that the variational problem (VP) has one and only one solution in W and that the following a priori estimate holds.
We also define the approximation u (k) h of u solution to the approximate variational formulation: where we assume that W (k) h and V (k) h are two finite-dimensional subsets of W and V, respectively. Moreover, as observed in [7] (Remark 2.23, p.92), neither condition (BNB1) nor condition (BNB2) imply its discrete counterpart. Then, the well-posedness of (3) is equivalent to the two following discrete conditions.
h ), a direct application of Theorem 2.2 in [7] enables us to write the following a priori estimate.
In the next section, we will apply these results to the particular case where the exact solution u belongs to W m,p Sobolev spaces and the approximation u h is computed by the help of P k Lagrange finite elements.

Application to W m,p Sobolev Spaces and the Corresponding Error Estimate
We introduce an open-bounded subset Ω ⊂ R n exactly recovered by a mesh T h composed by N K n-simplexes K µ , (1 ≤ µ ≤ N K ), which respect the classical rules of regular discretization (see, for example, [4] or [8]). Moreover, we denote by h the mesh size of T h (the largest diameter in the mesh T h ) and the space of polynomials defined on a given n-simplex K µ of degree less than or equal to k (k ≥ 1) by P k (K µ ).
Henceforth, we assume that approximate spaces W h ) and that they are included in the space of functions defined on Ω and composed of polynomials belonging to P k (K µ ), (1 ≤ µ ≤ N K ).
Finally, we also specify the functional framework of the abstract problem (VP) by introducing W m,p Sobolev spaces as follows.
For any integer m ≥ 1 and any 1 ≤ p ≤ +∞, we denote by W m,p (Ω) the Sobolev space of (class of) real-valued functions that, together with all their partial distributional derivatives of order less or equal to m, belongs to L p (Ω): α = (α 1 , α 2 , . . . , α n ) ∈ N n being a multi-index for which its length |α| is given by |α| = α 1 + · · · + α n , and ∂ α u is the the partial derivative of order |α| defined by the following.
We also consider norm . m,p,Ω and semi-norms |.| l,p,Ω , which are, respectively, defined by the following: where . L p denotes the standard norm in L p (Ω).
Then, in order to fulfill the conditions of the (BNB) Theorem, particularly so that W is a Banach space and V a reflexive one, in the sequel of the paper, the following definitions of spaces W and V hold: where m and m are two non zero integers, and p and p two real positive numbers that satisfy p = 2 and p > 1 such that the following is the case.
Regarding these choices, Sobolev's space W m,p (Ω), . m,p,Ω is a Banach space and W m ,p (Ω), . m ,p ,Ω is a reflexive one [9]. Moreover, we have: We can now recall the W m,p a priori error estimate for P k Lagrange finite elements we derived in [6]: where C k is a positive constant independent of h.

Remark 1. Since we noticed that W (k)
h is included in W m,p (Ω), by considering for the topology of W (k) h that was induced from W m,p (Ω) thanks to the triangle inequality, (2)-(4) and (10) result in the following error estimate: where V is the duality of W m ,p (Ω). As one can observe, the right-hand side of (11) contains α (k) h for which its dependency on h is usually unknown, except for particular cases. This is the reason why we will assume it is bounded from below by a positive constant δ independent of h: This uniform boundedness property of α (k) h , crucial to guaranteeing optimal error estimates [7], is valid in multiple cases (see Chapters 4 and 5 in [7] or [10]) but not systematically (see for example first-order PDEs in [7]).
Next, from (11), we obtain the following: where 1 In the sequel, we will denote by β k the following expression.
Finally, we also observe that there exists a critical value k of h defined by the following: such that the error estimate (12) can be split according to the following.
Based on this error estimate, the following section is devoted to deriving a probabilistic model applied to the relative accuracy between two Lagrange finite elements P k 1 and P k 2 , (k 1 < k 2 ), when the mesh size h has a given and fixed value.

The Probabilistic Analysis of the Relative Finite Elements Accuracy
In [6], we proposed two probability distributions that enabled us to obtain an evaluation of the more likely W m,p accurate between two Lagrange finite elements P k 1 and P k 2 , (k 1 < k 2 ). These distributions were essentially derived by considering the error estimate (10).
In the present paper, we generalize these probability distributions by introducing two new inputs, which are the following: 1.

2.
The probabilistic approach we develop is an extension of those we considered in [1,6]. More precisely, by the help of an ad hoc density probability, we derive the probability distribution of a suitable random variable so that we can compare two approxima- , considered as random variables, as will be introduced now.

Random Solution and Random Approximations of Deterministic Partial Differential Equation
The purpose of this section is based of the following fundamental remark: the solution u to the variational problem (VP), except for particular cases, is totally unknown (being impossible to calculate it analytically); this motivates the numerical schemes one will choose to implement.
This inability to determine, in most cases, the exact solution u is mainly due to the complexity of the involved PDE's operator; it indeed depends on complex combinations of integrals, partial derivatives and boundary conditions, as well as on the bent geometrical shape of the domain of integration Ω. All of these ingredients hence participate in the incapability to analytically determine exact solution u as their relationship with u is inextricable and unknown.
As a consequence, this lack of knowledge and information regarding the dependency between these ingredients and solution u motivates us to consider u as a random variable, as well as any function of u. This paper is dedicated to the W m,p approximation error of h , considered as a random variable. In this frame, we view solution u with respect to variational formulation (VP) defined by (1) in the same manner as it normally is viewed to consider the trajectory and contact point with the ground of any solid body that is thrown, i.e., as random. Indeed, in this case, due to the lack of information concerning the initial conditions of the trajectory of the body, the solution of the concerned inverse kinematic operator is inaccessible and is, thus, observed as a random variable.
In the case analyzed in this paper, the situation is much worse. Indeed, we investigate a general variational formulation (VP) where the analog of the inverse kinematic operator is too complex to enable us to analytically determine the corresponding solution of (VP) by any mathematical expression. It is one of the reasons that motivate us to view solution u and its approximation u (k) h as random variables, since the corresponding approximate operator conserves the complexity of the original one described above.

The Probabilistic Distribution of the Relative Finite Elements Accuracy
In the previous section, we expressed the reason for why we consider solution u and its approximation u (k) h as random variables. In order to complete the description of the randomness feature of approximation error ,Ω , we also remark that since the manner a given mesh grid generator will produce any mesh is random, the corresponding approximation u (k) h is random too. For all of these reasons, based on the error estimate (12) and (13), we can only affirm that the value of approximation error u As a consequence, we decide to view norm u − u (k) h m,p,Ω as a random variable defined as follows.
Let k, m and p be fixed. We introduce the random variable X (k) m,p defined by the following.
h plays the role of the usual probability space introduced in this context. Now, regarding the absence of information concerning the more likely or less likely Equation (18)  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 ) ∈ N 2 such that 0 < k 1 < k 2 .
The two corresponding inequalities given by (12) and (13), assuming that solution u to (VP) belong to H k 2 +1 (Ω), are as follows: where u h , respectively, denote P k 1 and P k 2 Lagrange finite element approximations of u and β k i , (i = 1, 2), as defined by (13).

Remark 2.
If one considers a given mesh for finite element P k 2 that contains that of P k 1 , then for the particular class of problems where (VP) is equivalent to a minimization formulation (MP) (see for example [11]), one can show that the approximation error of the P k 2 finite element is always smaller than that of P k 1 , and P k 2 is more accurate than P k 1 for all values of the mesh size h.
Therefore, in order to avoid this situation, for a given value of h, we consider two independent meshes built by a mesh generator for P k 1 and P k 2 . Now, usually, in order to compare the accuracy between these two finite elements, one asymptotically considers inequalities (19) and (20) to conclude that, when h proceeds to zero, P k 2 is more accurate than P k 1 , since h k 2 proceeds faster to zero than h k 1 .
However, for a given application, h has a given and fixed value; thus, this method of comparison is not valid anymore. For this reason, our purpose is to determine the relative accuracy between two finite elements P k 1 and P k 2 , (k 1 < k 2 ) for a fixed value of h corresponding to two independent meshes.
Moreover, since we chose to consider two random variables X (k i ) m,p , (i = 1, 2), as uniformly distributed on their respective interval of values [0, β k i ], (i = 1, 2), we also assume that they are independent. This assumption is, once again, the result of the lack of information which led us to model the relationship between these two variables as independent, since any knowledge is available for more precisely localizing the value of u − u m,p,Ω was known and vice versa. By the following result, we establish the density of probability of the random variable Z defined by the following: , be the two uniform and independent random variables defined by (16) and (17) (13). Then, random variable Z defined on R has the following density of probability.
Proof. Let us now remark that since the support of two random variables X (k i ) m,p , (i = 1, 2), is [0, β k i ], the support of density f Z is, therefore, [−β k 2 , β k 1 ], which corresponds to (21).
Let us consider the case where β k 1 ≤ β k 2 .
If f X (k i ) m,p (x i ) denotes the density of probability defined on R associated to each random variable X (k i ) m,p , since we assume that they are independent variables, density f Z (z) is given by the following: where f Furthermore, due to the definition (27) of the density for each variable X (k i ) m,p , (i = 1, 2), the integrand of (26) can be expressed as follows: As a consequence, , which leads one to consider the five following cases corresponding to the significant relative positions between intervals [0, β k 2 ] and [−z, −z + β k 1 ]: 1.
Let us assume that −z = ∅ and f Z (z) = 0, which is again a result of (21);
The following case concerns the values of z such that −z ≥ 0 and −z and we obtain (23) since the following is the case. 4.
The other cases corresponding to β k 1 ≥ β k 2 can be deduced by using the same arguments.
From Theorem 1, one can infer the entire cumulative distribution function F Z (z) defined by the following.
However, since we are interested in determining the more likely finite element between P k 1 and P k 2 , (k 1 < k 2 ), we focus the following corollary to the value of F Z (0), which corresponds to Prob X  (13). Then, we have: Proof.
-Let us consider the case where β k 1 ≤ β k 2 . Then, by definition (33) of the entire cumulative distribution function F Z (z), at z = 0, we have the following: where we used the values of the density f Z (z) given by (22) and (23). -In the same manner, when β k 1 ≥ β k 2 , we have the following.
We can now explicate the probability distribution of event X given by (34) and (35) as a function of the mesh size h.
To this end, we remark that, since each β k i (i = 1, 2) has two possible values depending on the relative position between h and k i defined by (14), the probability distribution we are looking for must be splitted in the corresponding cases as well.
Let us, hence, introduce constants C k i , (i = 1, 2), defined by the following: and the specific value h * k 1 ,k 2 of h which corresponds to the intersection of the curves ϕ k i , (i = 1, 2), defined by ϕ k i (h) ≡ C k i h k i +1−m (see Figure 1). Then, we have the following. (39) We notice that h * k 1 ,k 2 and k i , (i = 1, 2), strongly depend on m and p, since C k i and β k i depend on these two parameters as well. As a consequence, in the following theorem, the different formulas of Prob X will contain this dependency on m and p.
m,p , (i = 1, 2), be the two uniform and independent random variables defined by (16) and (17). Then, P k 1 , is determined by the following.
If h * k 1 ,k 2 ≥ max( k 1 , k 2 ), then k 1 ≤ k 2 and we have: Proof. In order to establish the proof of Theorem 2, we will consider a geometrical interpretation of error estimate (19) and (20). These two inequalities can indeed be geometrically viewed with the help of the relative position between the two curves ϕ k i (h), (i = 1, 2), introduced above and the horizontal line defined by ψ(h) ≡ l V α * (see Figure 1). Then, depending on the position of the horizontal line ψ(h) with particular value C k i h * k 1 ,k 2 , (i = 1, 2), we have to consider the two following cases: Thus, let us consider the first case when In this case, we showed in [12] two numerical examples based on PDE's formulation where exact solutions may be computed to illustrate probabilistic law (47) and (48).
In particular, we highlighted the statistical methodology that enabled us to determine a statistical estimator of the threshold h * k 1 ,k 2 given by (39). Based on this study, one may adapt it in order to obtain applications that will illustrate the new probabilistic value we derived in Theorem 2.
Moreover, in Theorem 2, the contribution of a priori estimates (2) and (4) modifies the probability distribution given by (47) and (48), since the horizontal line ψ(h) = l V α * interferes with the two polynomials C k i h k i +1−m , (i = 1, 2) (see Figure 1). This phenomenon will be analyzed in the following section.

Discussion
First of all, let us plot the shape of the two cases of the probability distribution we derived in Theorem 2 and those we obtained in [1,6] that we recalled in (47) and (48).
As one can observe, the distribution (47) and (48) plotted in Figure 4 resembles more of the second case of Theorem 2 plotted in Figure 3 than its first case, as plotted in Figure 2.
In fact, in the latter case, when k 1 ≤ k 2 , line ψ(h) = l V α * interacts with two polynomials C k i h k i +1−m , (i = 1, 2), before critical value h * k 1 ,k 2 . As a consequence, its contribution correspondingly takes place in the expression of the probability distribution for h ≥ k 1 .
On the contrary, when k 1 ≥ k 2 , the contribution of line ψ(h) = l V α * only has to be taken into account after the value of h * k 1 ,k 2 . However, in this case, when h ≥ k 2 , once again line ψ(h) interacts with two polynomials C k i h k i +1−m , (i = 1, 2), and the equivalent distribution we found in (47) and (48) is cut to provide the corresponding one in (45) and (46) plotted in Figure 3.  This shows the importance of the a priori estimates (2) and (4), which are not considered in the classical point of view limited to the asymptotic behavior of error estimate (10), when the mesh size h proceeds to zero. The reason for this is the fact that the right hand sides of (2) and (4) do not depend on h; it does not bring any influence for the desired asymptotic behavior.
The goal of the probabilistic approach we developed in this paper is to evaluate, for any fixed value of the mesh size h, the relative accuracy between two Lagrange finite elements P k 1 and P k 2 , (k 1 ≤ k 2 ). The probability distributions of Theorem 2 adjust those that we had found in [1,6].
Indeed, probability distribution (47) and (48) claim that finite element P k 1 is more likely to be accurate than P k 2 with a high level of probability when h ≥ h * k 1 ,k 2 . Actually, this probability can reach the value of one, as it can be observed in Figure 4 and P k 1 could be almost surely more accurate than P k 2 . Thus, distributions (40)-(42) and (43)-(46) limit the value of this probability. More precisely, if k 1 ≤ k 2 , then the probability such that P k 1 is more likely accurate than P k 2 will never be greater than 0.5 (see Figure 2), and if k 1 ≥ k 2 , then for h ≥ k 2 the corresponding probability will decrease to the limit value of 0.5, even if it increases with values greater than 0.5 for h values between h * k 1 ,k 2 and k 2 (see Figure 3).

The Relative Accuracy between Two Finite Elements
In this paper, we presented a generalized probability distribution that takes into account the three available standard estimates one can confront with respect to solution u to a variational formulation together with its approximation u (k) h performed by Lagrange finite elements P k .
Those include the a priori estimates for solution u and its approximation u (k) h carried out by finite elements, on the one hand, and the error estimate, on the other hand.
As we observed in the previous paragraph, if the contribution of a priori estimates (2) and (4) are not considered when one limits the analysis of the relative accuracy between two finite elements to the asymptotic rate of convergence when the mesh size h proceeds to zero, one must deal with these two estimates in order to conclude about the more likely accurate one when h has a fixed value.
Moreover, regarding the probabilistic approach we developed in this paper, the a priori estimates (2) and (4) brought a more realistic behavior of distributions (40)-(46) in comparison with (47) and (48). Indeed, the probability given by (40)-(46) cannot reach a value of one since it was quite a surprise in (47) and (48) to obtain this theoretical value. This result meant that event "P k 1 is more accurate than P k 2 " is an almost certain event when k 1 < k 2 .

Generalization to Other Numerical Methods
Finally, we would like to mention that the probabilistic approach we proposed in this study is not restricted to the finite elements method but may be extended to other types of approximation: given a class of numerical schemes and their corresponding approximation 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 one.
It is, for example, the case of the numerical integration where the composite quadrature error has a mathematical structure, which looks similar to error estimate (10) we considered in the present study.
More precisely, as an example, for a composite quadrature of order k on a given interval is a given function, the corresponding composite quadrature error can be written [13] as follows: where h denotes the size of the equally spaced panels that discretized the interval [a, b], (λ i , x i ), (i = 0, N), are (2N + 2) given numbers such that all x i belong to [a, b], and C k is a constant independent of h that mainly depends on f and k.
As a consequence, due to the similar mathematical structure between (10) and (49), with the same arguments we introduced to compare the accuracy between two Lagrange finite elements P k 1 and P k 2 , (k 1 < k 2 ), one can evaluate the probability of the more accurate numerical composite quadrature associated to two different parameters k 1 and k 2 , (k 1 ≤ k 2 ) for a fixed value of h.
This will make sense because, normally, to bypass the lack of information associated with the unknown value of the left hand side of (49) in interval [0, C k h k+1 ], only the asymptotic convergence rate comparison is concerned for evaluating the relative accuracy between two numerical quadratures of order k 1 and k 2 , (k 1 < k 2 ).
Nevertheless, this procedure is no longer available when one wants to compare two composite quadratures in the case when the size of the equally panels h is fixed, as it is for any application. Thus, the probabilistic approach we propose here could be a relevant alternative.
The same consideration may be developed in order to compare the accuracy between two numerical schemes of order k 1 and k 2 , which are performed to approximate the exact solution of ordinary differential equations. To make precise these ideas, let us consider solution y of the first order ordinary differential initial value problem defined on a given interval [t 0 , t 0 + T]: (CP) y (t) = f (t, y(t)), t ∈ [t 0 , t 0 + T], y(t 0 ) = y 0 , where y 0 is given. Let us also restrict ourselves by considering one-step numerical methods to approximate function y solution to problem (CP). Namely, if we introduce a constant mesh size h = t n+1 − t n , where (t n ) n=0,N denotes the sequence of (N + 1) values of t within the interval [t 0 , t 0 + T], then the corresponding one-step numerical scheme is given by the following: (CP) h y n+1 = y n + hΦ(t n , y n ), n > 0, where Φ is a given function "sufficiently" smooth that characterizes the numerical scheme (CP) h . Moreover, (CP) h is called a numerical scheme of order k if [13]: |y(t n+1 ) − y(t n ) − hΦ(t n , y(t n ))| ≤ C k h k , where C k is a constant independent of h that depends on y, Φ and k (see [14] for the dependency on k). Thus, when considering two one-step numerical methods of order k 1 and k 2 , (k 1 < k 2 ), defined by two functions Φ 1 and Φ 2 , and due to the similar structure between (10) and (52), one would be able to evaluate the probability of the more accurate scheme with the same arguments we implemented when comparing accuracy between P k 1 and P k 2 Lagrange finite elements.
In summary, when one wants to evaluate the relative accuracy between two numerical methods that belong to a given family of approximations, the probabilistic approach we propose in this study essentially depends on the ability to determine the constant C k , which appears in different corresponding approximation errors (10), (49) or (52).
Indeed, for each of these errors of approximation, the complexity of the constant C k will suggest appropriate investigations. In the current study, we were devoted to the finite elements method, and we pointed out the important role of the a priori estimates (2) and (4), which are usually neglected, since they do not bring any asymptotic information/behavior when mesh size h proceeds to zero.
As a consequence, in order to render probability distribution (40)-(46) derived in Theorem 2, one will have to consider appropriate techniques that will make the determination (or at least the approximation) of different constants that are involved in (40)-(46) possible. This is mainly due to the contribution of the three estimates (2), (4) and (10), which results in (12), namely, k 1 , k 2 and h * k 1 ,k 2 .
Funding: This research received no external funding.
Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Not applicable.