Some New Results Concerning the Classical Bernstein Cubature Formula

: In this article, we present a solution to the approximation problem of the volume obtained by the integration of a bivariate function on any ﬁnite interval [ a , b ] × [ c , d ] , as well as on any symmetrical ﬁnite interval [ − a , a ] × [ − a , a ] when a double integral cannot be computed exactly. The approximation of various double integrals is done by cubature formulas. We propose a cubature formula constructed on the base of the classical bivariate Bernstein operator. As a valuable tool to approximate any volume resulted by integration of a bivariate function, we use the classical Bernstein cubature formula. Numerical examples are given to increase the validity of the theoretical aspects.


Introduction
Let N be the set of positive integers and N 0 = N ∪ {0}. Recent studies concerning the classical Bernstein operators associated with any real-valued function F : [a, b] → R, where a, b are two real and finite numbers, such that a < b can be found in a series of research papers, see [1][2][3][4][5][6][7][8][9][10][11][12], which shows that there is sustained interest in the study of the classical linear positive operators. These operators are given by for any x ∈ [a, b] and any n ∈ N. If a := 0 and b := 1 in (1), then which are the well-known Bernstein operators [13] associated with any real-valued function f : [0, 1] → R, any x ∈ [0, 1], and any n ∈ N. The Bernstein polynomials (2) opened a new era in the approximation theory beginning with Sergei Natanovich Bernstein's famous proof of the Weierstrass approximation theorem in 1912, and continuing until today with thousands of interesting papers. Some generalizations (approximation of integrable functions, approximation of measurable functions, Bernstein polynomials on an unbounded interval, degenerate Bernstein polynomials), as well as many other applications of the Bernstein polynomials, can be consulted in the excellent book [14]. An exceptional historical perspective is provided in [6] on the introduction and evolution of the Bernstein polynomials. Constructed as an important tool in approximation of functions by polynomials, they became in time more useful in exploiting computers to interactively design polynomial functions. In this context, the Bernstein polynomial basis provides valuable insight into its behavior over a certain finite interval, yielding many useful properties and elegant algorithms that are now being increasingly adopted in other application domains. In the real daily applications, different situations arise where it is necessary to solve certain definite integrals, much more complicated than those presented in the courses of mathematical analysis. In order to solve such integrals in an efficient way, numerous numerical methods called quadrature formulas have been developed. The most known method of numerical integration is obtained by integrating the Lagrange interpolation formula resulting in the class of Newton-Cotés quadrature formulas. Engineers and researchers use numerical integration as a primary tool to get approximate values for definite integrals that cannot be solved analytically. The quadrature formulas based on the linear positive Bernstein type operators represented a theme of intensive study for the distinguished mathematician Dimitrie D. Stancu. His contributions in this area were collected in the second volume of the monograph "Numerical Analysis and Approximation Theory " [15]. Looking deep at this problem of numerical integration based on the Bernstein polynomial, we remarked that it is not fully researched. In [10], we managed to assign a place on the map of the closed  [11]. Although effective in most situations, there are instances when the composite quadrature formulas cannot be applied, as they use equally-spaced nodes. To avoid this fact, we constructed an adaptive method in [10].
having the graphs (see Figure 1 for F 1 and Figure 2 for F 2 ) given below.  In general, the graph of a bivariate function F is a surface having the equation z = F(x, y) as we can see in Figure 3. Therefore, the definite integral of a positive function with two variables is a double integral, which represents the volume of the solid between the surface of the function and the plane that contains its domain (see Figure 3). In this article, we present a solution to the approximation problem of the volume obtained by the integration of a bivariate function when a double integral cannot be computed exactly. The approximation of various double integrals can be done by a few cubature formulas (for instance, the Newton-Cotés cubature formulas) according to the specialty literature. Constructed by means of the bivariate Lagrange polynomial, trapezoidal and Simpson cubature formulas use a fixed number of nodes, resulting in a single possible approximation for a double integral. In order to be more flexible with this fact, we bring to light a cubature formula constructed on the base of the classical bivariate Bernstein operator. As a valuable tool to approximate any volume resulting by integration of a bivariate function, we use the classical Bernstein  Numerical examples are given to increase the validity of the theoretical aspects.

Auxiliary Results
For any n ∈ N, let a ≤ x 0 < x 1 < · · · < x n ≤ b be some distinct nodes, and let F : [a, b] → R be a real-valued function. In many books on numerical analysis, divided differences are defined recursively If I, J are certain real intervals, F : I × J → R is a real-valued function and (x 0 , y 0 ), (x 1 , y 1 ) ∈ I × J (x 0 = x 1 , y 0 = y 1 ), then the bivariate divided differences of F with respect to the nodes (x 0 , y 0 ), (x 1 , y 1 ) are defined in [16], using the method of parametric extensions by Other equivalent definitions for univariate and bivariate divided differences can be found in the excellent monograph [17]. In the definition of the divided differences, the number of abscissas in general is not equal to the number of ordinates. If x 0 , x 1 , . . . , x p ∈ I and y 0 , y 1 , . . . , y q ∈ J are distinct nodes, then the following recurrence formula holds (see [16]), for p, q ∈ N, p, q ≥ 2 and x 0 , x 1 , . . . , x p y 0 , y 1 , . . . , y q ; F = x i 0 , x i 1 , . . . , x i p y j 0 , y j 1 , . . . , y j q ; F , where (i 0 , i 1 , . . . , i p ), (j 0 , j 1 , . . . , j q ) are permutations of (0, 1, . . . , p) and (0, 1, . . . , q), respectively. The extension of the divided differences in two dimensions that are used in this paper is a particular case of both T. Popoviciu's definitions [18] and Popoviciu's particular case [18] on networks of an M.A. Màrchàud type is defined in [19]. In [20], W.J. Gordon introduced the basic notions of the algebraic theory of multivariate functions approximation, a theory that was studied and further developed by F.J. Delvos and W. Schempp [21]. The method of parametric extension is a procedure for constructing linear operators on the spaces of multivariate functions starting from the linear operators defined on the spaces of univariate functions. Let S = [a, b] × [c, d] be a rectangular domain and suppose that F : S → R is given. The parametric extensions of the classical Bernstein operators (1) to the space C(S) are defined for x, y ∈ S and n 1 , n 2 ∈ N by and For brevity, we can use the following two notations: The tensor product of the parametric extensions (4) and (5) is the classical bivariate Bernstein operators The equality F(x, y) = B n 1 ,n 2 (F; x, y) + R n 1 ,n 2 (F; x, y) is called the classical bivariate Bernstein approximation formula, and R n 1 ,n 2 are the bivariate remainder operators associated with the classical bivariate Bernstein operators (6).

The Classical Bernstein Cubature Formula
We firstly establish some new results concerning the remainder term of the classical bivariate Bernstein approximation formula (7). The following two results play an important role in building a valuable approximation formula for certain double integrals. Theorem 1. The representation of the remainder term associated with the classical bivariate Bernstein operators (6) is given by Proof. Starting with the classical bivariate Bernstein approximation formula (7) and based on the fact that bivariate operators (6) preserve constants, we may write Putting an emphasis on the following equality In a similar way, one obtains As a direct consequence of the above theorem, we get the simplified form of the remainder term in the approximation formula (7), as well as an upper bound estimation. In the following lines, F (o,p) are the partial derivatives of order (o, p).
Proof. The form of the remainder obtained in the above theorem is being processed.
Applying the mean value theorem on the bivariate divided differences leads us to the equality. Furthermore, by using modulus and the fact that partial derivatives of the function F are bounded on S, we get the inequalities.
Based on the results proved above, we can integrate the classical bivariate Bernstein approximation formula (7) Proof. Using the definition of the classical bivariate Bernstein operators (6) Taking into account the definition and the properties of the Euler's integral of first kind Collecting the results demonstrated above and eliminating the remainder term in the relation ( Theorem 3. Assuming that F ∈ C 2,2 (S) and ξ, η ∈ (a, b) × (c, d), then Proof. Both forms of the remainder term presented in Corollary 1 get integrated on S.
Applying the classical bivariate Bernstein formula (10) with n 1 = n 2 = 1 (the simplest case), we get which is the exact value of the given integral.

Remark 2.
The veracity of the above statement (Case 2) is based on the fact that the classical bivariate Bernstein polynomial of the first degree is equal to the bivariate Lagrange polynomial of the first degree To demonstrate the approximation accuracy of the classical Bernstein cubature formula we test it with the Maple mathematical software. For our study, we consider five test functions as follows: Let us denote the approximative values of the integrals with ten exact decimals by and the appropriate absolute errors eI n 1 ,n 2 (F) = I F − I n 1 ,n 2 (F) .
We approximate the double integrals and record in Table 1 the absolute errors obtained by applying the classical Bernstein cubature formula (11). In Table 1, e − p means 10 −p . Analyzing the results recorded in Table 1, we remark how unconstrained the classical Bernstein cubature formula (11) is. This is given by the possibility to approximate the double definite integrals using a certain number of nodes. In the second column of Table 1, we used the simplest form of the classical Bernstein cubature formula (n 1 = n 2 = 1), which is in fact the trapezoidal cubature formula, see Case 2. In the next columns, we obtained better approximations of the double integrals due to the freedom to increase the number of nodes. Based on the results presented above, we are able to put an emphasis on some advantages of the use of the classical bivariate Bernstein formula (10): • The freedom to increase the number of nodes (in particular, a number of nodes greater than 2 2 ); • In evaluation of the remainder, we request a single condition for the integrated bivariate function, namely to belong to the space C 2,2 (S); • A double integral could be approximated with a desired precision ε, where ε > 0 is a real number.
Remark 3. The above advantages are nonexistent for the well-known Newton-Cotés bivariate formulas, as any of them has a fixed number of nodes, also resulting in a single possible approximation for a double integral.

Remark 4.
In particular, a special attention goes to the results shown in Table 1, the last column. Although a large number of nodes were used, one remarks a slow approximation of the integrals (exception I F 3 ) due to: • The classical bivariate Bernstein polynomials that interpolate the functions only at the end points of intervals.

•
The convex and concave shape of surfaces with large functional variation.
In order to avoid the disadvantages pointed out in Remark 4 and to reduce the number of nodes necessary for computing an approximation of a given double integral, we present the next section.

The Classical Composite Bernstein Cubature Formula
As we already said, our scope is to reduce the number of nodes necessary for computing an approximation of a given double integral.
In the following lines, the symbols n 1 , n 2 will refer to the polynomial degree, while the numbers m 1 , m 2 will be related to the grid (interval discretization). For n 1 , n 2 , m 1 , m 2 ∈ N, we consider the parametric extensions The tensor product of the parametric extensions defined above is x B n 1 ,m 1 ⊗ y B n 2 ,m 2 (F; x, y) =: B(F; x, y), having the following expression: The classical bivariate Bernstein approximation formula on the interval I i × I j , i = 1, m 1 and j = 1, m 2 becomes F(x, y) = B(F; x, y) + R(F; x, y), with R(F) the bivariate remainder operators associated with the tensor product B(F). Taking Corollary 1 into account and supposing that F ∈ C 2,2 I i × I j , i = 1, m 1 , j = 1, m 2 , we can present an upper bound estimation of the remainder term given by Integrating the classical bivariate Bernstein approximation formula (12) on I i × I j , i = 1, m 1 , and j = 1, m 2 , we get Summing the relations (14) for i = 1, m 1 and j = 1, Proof. Let i = 1, m 1 and j = 1, m 2 be fixed. Using the definition of the tensor product B(F), we have For brevity, we designate Proof. Starting with the relation (15), we have Using the relation (13), it follows that Summarizing the above results, we get the classical bivariate Bernstein composite formula with an upper bound estimation of the remainder term (assuming F ∈ C 2,2 (S)) Below, we present some particular cases of the classical bivariate Bernstein composite formula.
Remark 5. The corrected upper bound estimation (19) of the remainder term in the bivariate Bernstein composite formula was established in [1].

Numerical Examples
In order to validate the approximation accuracy of the classical composite Bernstein F( , •) =: I m 1 ,m 2 ,n 1 ,n 2 (F), (20) we test it with the Maple mathematical software. For our study, we reconsider the five test functions presented in the third section. We approximate the double integrals, and, in the next tables, we record the absolute errors eI m 1 ,m 2 ,n 1 ,n 2 (F) = I F − I m 1 ,m 2 ,n 1 ,n 2 (F) , obtained by applying the classical composite Bernstein cubature formula (20). In the following tables, e − p means 10 −p .
Analyzing carrefully the results recorded in the above five tables (Tables 2-6), we can remark the following aspects: 1.
In the third column of each table, we applied the simplest form of the classical composite Bernstein cubature formula (n 1 = n 2 = 1), which is in fact the composite trapezoidal cubature formula, see Case 4.

2.
Once the number of nodes increased, in the next columns of each table, we obtained better approximation than the previous one.
It is clear that the classical bivariate Bernstein composite formula (18) enjoys some privilegies: • The freedom to increase the number of nodes on all bivariate subintervals; • A single condition is required in the evaluation of the remainder; • A double integral could be approximated with a desired precision by taking a certain number of bivariate subintervals into account as well as a certain number of nodes on all subintervals.     The implementation on the computer can be done using a simple procedure; • The computational cost is limited to summing the values of functions; • A single condition for the integrated function is requested in the evaluation of the remainder; • A composite cubature formula can be used if the approximation is not good enough.
In this research, we were successful in assigning the classical Bernstein cubature formula a right place on the map of cubature formulas. Unconstrained and more flexible than the Newton-Cotés cubature formulas, the classical Bernstein cubature formula can be used in almost all real daily applications offering accurate solutions. Although effective in most situations, there will be instances when the implementation of the classical Bernstein cubature formula would be inappropriate, as a small step size is used uniformly across the entire interval of integration to ensure the overall accuracy. Such an approach does not take into account that some regions of the surface have large functional variations that require more attention than other regions of the surface. Under those circumstances and as a new direction of research, it would be useful to introduce a method that adjusts the step size to be smaller over regions of the surfaces where a larger functional variation occurs. An efficient problem solving technique for these types of issues should predict the amount of functional variation and adapt the step size as necessary.
Funding: This research received no external funding.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.
Data Availability Statement: All data generated and analysed during the study are included in this published article.