Some New Results Concerning the Classical Bernstein Cubature Formula

Abstract

In this article, we present a solution to the approximation problem of the volume obtained by the integration of a bivariate function on any finite interval $[a,b]\times [c,d]$, as well as on any symmetrical finite interval $[-a,a]\times [-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.

1. Introduction

Let $\mathbb{N}$ be the set of positive integers and ${\mathbb{N}}_0=\mathbb{N}\cup \left\{0\right\}$. Recent studies concerning the classical Bernstein operators associated with any real-valued function $F:[a,b]\to \mathbb{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

$$B_n(F;x)=\frac{1}{{(b-a)}^n}\sum _{k=0}^n\left(\genfrac{}{}{0pt}{}{n}{k}\right){(x-a)}^k{(b-x)}^{n-k}F\left(a+\frac{k(b-a)}{n}\right)$$

(1)

for any $x\in [a,b]$ and any $n\in \mathbb{N}$. If $a:=0$ and $b:=1$ in (1), then

$$B_n(f;x)=\sum _{k=0}^n{p}_{n,k}\left(x\right)f\left(\frac{k}{n}\right)=\sum _{k=0}^n\left(\genfrac{}{}{0pt}{}{n}{k}\right)x^k{(1-x)}^{n-k}f\left(\frac{k}{n}\right),$$

(2)

which are the well-known Bernstein operators [13] associated with any real-valued function $f:[0,1]\to \mathbb{R}$, any $x\in [0,1]$, and any $n\in \mathbb{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 Newton–Cotés quadrature formulas to a new approximation formula based on the classical Bernstein polynomials. In order to increase the validity of the theoretical aspects, we tested a series of numerical examples. By approximation of definite integrals on some intervals $[a,b]$, we observed that the resulting error is proportional with the interval length. Therefore, if the interval $[a,b]$ is large, then it is advisable to divide $[a,b]$ into smaller intervals and use the classical Bernstein quadrature formula on each subinterval. By adding the results, we get the approximation of the definite integral on the whole interval $[a,b]$. Such a scheme is called the classical composite Bernstein quadrature formula, and it can be found in [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]. All the mentioned results were established in the univariate case. In natural logic, we can continue our study in the bivariate case. Our aim is to compute the definite integral (on any finite bivariate interval $[a,b]\times [c,d]$ as well as on any symmetrical finite bivariate interval $[-a,a]\times [-a,a]$) of a function depending on two variables. For instance, let us consider the following two functions:

$$\begin{array}{c}{F}_1:[-5,5]\times [-5,5]\to \mathbb{R},F_1(x,y)=x^2+y^2+cosx+siny,\\ F_2:[-10,10]\times [-10,10]\to \mathbb{R},F_2(x,y)=e^{-(x^2+y^2)},\end{array}$$

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 cubature formula. If the bivariate interval $[a,b]\times [c,d]$ (the bivariate symmetrical interval $[-a,a]\times [-a,a]$) is large, then the classical composite Bernstein cubature formula is suitable for the approximation of a double integral. Numerical examples are given to increase the validity of the theoretical aspects.

2. Auxiliary Results

For any $n\in \mathbb{N}$, let $a\le x_0<x_1<\cdots <x_n\le b$ be some distinct nodes, and let $F:[a,b]\to \mathbb{R}$ be a real-valued function. In many books on numerical analysis, divided differences are defined recursively

$$[x_0;F]=F\left(x_0\right),\dots ,[x_0,x_1,\cdots ,x_n;F]=\frac{1}{x_n-x_0}\left([x_1,\cdots ,x_n;F]-[x_0,\cdots ,x_{n-1};F]\right).$$

If $I,J$ are certain real intervals, $F:I\times J\to \mathbb{R}$ is a real-valued function and $(x_0,y_0),(x_1,y_1)\in I\times J$ ($x_0\ne x_1$, $y_0\ne 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

$$\left[\begin{array}{c}{x}_0,x_1\hfill \\ y_0,y_1\hfill \end{array};F\right]=\frac{F(x_1,y_1)-F(x_0,y_1)-F(x_1,y_0)+F(x_0,y_0)}{(x_1-x_0)(y_1-y_0)}.$$

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,\dots ,x_p\in I$ and $y_0,y_1,\dots ,y_q\in J$ are distinct nodes, then the following recurrence formula

$$\begin{array}{cc}\hfill \left[\begin{array}{c}{x}_0,x_1,\dots ,x_p\hfill \\ y_0,y_1,\dots ,y_q\hfill \end{array};F\right]=\frac{1}{(x_p-x_0)(y_q-y_0)}& \left(\left[\begin{array}{c}{x}_1,\dots ,x_p\hfill \\ y_1,\dots ,y_q\hfill \end{array};F\right]-\left[\begin{array}{c}{x}_0,\dots ,x_{p-1}\hfill \\ y_1,\dots ,y_q\hfill \end{array};F\right]\right.\hfill \\ & \left.-\left[\begin{array}{c}{x}_1,\dots ,x_p\hfill \\ y_0,\dots ,y_{q-1}\hfill \end{array};F\right]+\left[\begin{array}{c}{x}_0,\dots ,x_{p-1}\hfill \\ y_0,\dots ,y_{q-1}\hfill \end{array};F\right]\right)\hfill \end{array}$$

holds (see [16]), for $p,q\in \mathbb{N}$, $p,q\ge 2$ and

$$\left[\begin{array}{c}{x}_0,x_1,\dots ,x_p\hfill \\ y_0,y_1,\dots ,y_q\hfill \end{array};F\right]=\left[\begin{array}{c}{x}_{i_0},x_{i_1},\dots ,x_{i_p}\hfill \\ y_{j_0},y_{j_1},\dots ,y_{j_q}\hfill \end{array};F\right],$$

(3)

where $(i_0,i_1,\dots ,i_p)$, $(j_0,j_1,\dots ,j_q)$ are permutations of $(0,1,\dots ,p)$ and $(0,1,\dots ,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]\times [c,d]$ be a rectangular domain and suppose that $F:S\to \mathbb{R}$ is given. The parametric extensions of the classical Bernstein operators (1) to the space $C\left(S\right)$ are defined for $x,y\in S$ and $n_1,n_2\in \mathbb{N}$ by

$${}_x{B}_{n_1}(F;x,y)=\sum _{k_1=0}^{n_1}\left(\genfrac{}{}{0pt}{}{n_1}{k_1}\right)\frac{{(x-a)}^{k_1}{(b-x)}^{n_1-k_1}}{{(b-a)}^{n_1}}F\left(a+\frac{k_1(b-a)}{n_1},y\right)$$

(4)

and

$${}_y{B}_{n_2}(F;x,y)=\sum _{k_2=0}^{n_2}\left(\genfrac{}{}{0pt}{}{n_2}{k_2}\right)\frac{{(y-c)}^{k_2}{(d-y)}^{n_2-k_2}}{{(d-c)}^{n_2}}F\left(x,c+\frac{k_2(d-c)}{n_2}\right).$$

(5)

For brevity, we can use the following two notations:

$${}_x{B}_{n_1}(F;x,y)=\sum _{k_1=0}^{n_1}{p}_{n_1,k_1}\left(x\right)F\left(t_1\left(k_1\right),y\right),{}_y{B}_{n_2}(F;x,y)=\sum _{k_2=0}^{n_2}{p}_{n_2,k_2}\left(y\right)F\left(x,t_2\left(k_2\right)\right).$$

The tensor product of the parametric extensions (4) and (5) is the classical bivariate Bernstein operators

$$\left({}_x{B}_{n_1}\otimes {}_y{B}_{n_2}\right)(F;x,y)=B_{n_1,n_2}(F;x,y)=\sum _{k_1=0}^{n_1}\sum _{k_2=0}^{n_2}{p}_{n_1,k_1}\left(x\right)·p_{n_2,k_2}\left(y\right)F\left(t_1\left(k_1\right),t_2\left(k_2\right)\right).$$

(6)

The equality

$$F(x,y)=B_{n_1,n_2}(F;x,y)+R_{n_1,n_2}(F;x,y)$$

(7)

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).

3. 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

$$\begin{array}{c}{R}_{n_1,n_2}(F;x,y)=-\frac{(x-a)(b-x)}{n_1}\sum _{k_1=0}^{n_1-1}\sum _{k_2=0}^{n_2}{p}_{n_1-1,k_1}\left(x\right)p_{n_2,k_2}\left(y\right)\left[\begin{array}{ccc}x,& t_1\left(k_1\right),& t_1(k_1+1)\\ & t_2\left(k_2\right)\end{array};F\right]\\ -\frac{(y-c)(d-y)}{n_2}\sum _{k_1=0}^{n_1}\sum _{k_2=0}^{n_2-1}{p}_{n_1,k_1}\left(x\right)p_{n_2-1,k_2}\left(y\right)\left[\begin{array}{ccc}& t_1\left(k_1\right)& \\ y,& t_2\left(k_2\right),& t_2(k_2+1)\end{array};F\right]\\ +\frac{(x-a)(b-x)(y-c)(d-y)}{n_1·n_2}\sum _{k_1=0}^{n_1-1}\sum _{k_2=0}^{n_2-1}{p}_{n_1-1,k_1}\left(x\right)p_{n_2-1,k_2}\left(y\right)\left[\begin{array}{ccc}x,& t_1\left(k_1\right),& t_1(k_1+1)\\ y,& t_2\left(k_2\right),& t_2(k_2+1)\end{array};F\right].\end{array}$$

Proof. 

Starting with the classical bivariate Bernstein approximation formula (7) and based on the fact that bivariate operators (6) preserve constants, we may write

$$\begin{array}{cc}\hfill R_{n_1,n_2}(F;x,y)& =F(x,y)-B_{n_1,n_2}(F;x,y)\hfill \\ & =\sum _{k_1=0}^{n_1}\sum _{k_2=0}^{n_2}{p}_{n_1,k_1}\left(x\right)·p_{n_2,k_2}\left(y\right)·\left(F(x,y)-F\left(t_1\left(k_1\right),t_2\left(k_2\right)\right)\right).\hfill \end{array}$$

Putting an emphasis on the following equality

$$\begin{array}{c}F(x,y)-F\left(t_1\left(k_1\right),t_2\left(k_2\right)\right)=F\left(x,t_2\left(k_2\right)\right)-F\left(t_1\left(k_1\right),t_2\left(k_2\right)\right)+F\left(t_1\left(k_1\right),y\right)\\ -F\left(t_1\left(k_1\right),t_2\left(k_2\right)\right)+F(x,y)-F\left(x,t_2\left(k_2\right)\right)\\ -F\left(t_1\left(k_1\right),y\right)+F\left(t_1\left(k_1\right),t_2\left(k_2\right)\right),\end{array}$$

we get

$$R_{n_1,n_2}(F;x,y)=S_1+S_2+S_3,$$

with

$$\begin{array}{c}{S}_1=\sum _{k_1=0}^{n_1}\sum _{k_2=0}^{n_2}{p}_{n_1,k_1}\left(x\right)·p_{n_2,k_2}\left(y\right)·\left(F\left(x,t_2\left(k_2\right)\right)-F\left(t_1\left(k_1\right),t_2\left(k_2\right)\right)\right)\\ =\sum _{k_1=0}^{n_1}\sum _{k_2=0}^{n_2}{p}_{n_1,k_1}\left(x\right)·p_{n_2,k_2}\left(y\right)·\left(x-t_1\left(k_1\right)\right)\left[\begin{array}{c}x,t_1\left(k_1\right)\\ t_2\left(k_2\right)\end{array};F\right]\end{array}$$

$$\begin{array}{c}=\sum _{k_1=0}^{n_1}\sum _{k_2=0}^{n_2}{p}_{n_1,k_1}\left(x\right)·p_{n_2,k_2}\left(y\right)\frac{(n_1-k_1)(x-a)-k_1(b-x)}{n_1}\left[\begin{array}{c}x,t_1\left(k_1\right)\\ t_2\left(k_2\right)\end{array};F\right]\\ =\sum _{k_1=0}^{n_1-1}\sum _{k_2=0}^{n_2}\left(\genfrac{}{}{0pt}{}{n_1-1}{k_1}\right)\frac{{(x-a)}^{k_1+1}{(b-x)}^{n_1-k_1}}{{(b-a)}^{n_1}}·p_{n_2,k_2}\left(y\right)\left[\begin{array}{c}x,t_1\left(k_1\right)\\ t_2\left(k_2\right)\end{array};F\right]\\ -\sum _{k_1=1}^{n_1}\sum _{k_2=0}^{n_2}\left(\genfrac{}{}{0pt}{}{n_1-1}{k_1-1}\right)\frac{{(x-a)}^{k_1}{(b-x)}^{n_1+1-k_1}}{{(b-a)}^{n_1}}·p_{n_2,k_2}\left(y\right)\left[\begin{array}{c}x,t_1\left(k_1\right)\\ t_2\left(k_2\right)\end{array};F\right]\\ =\frac{(x-a)(b-x)}{b-a}\sum _{k_1=0}^{n_1-1}\sum _{k_2=0}^{n_2}{p}_{n_1-1,k_1}\left(x\right)·p_{n_2,k_2}\left(y\right)\left(\left[\begin{array}{c}x,t_1\left(k_1\right)\\ t_2\left(k_2\right)\end{array};F\right]-\left[\begin{array}{c}x,t_1(k_1+1)\\ t_2\left(k_2\right)\end{array};F\right]\right)\\ =-\frac{(x-a)(b-x)}{b-a}\sum _{k_1=0}^{n_1-1}\sum _{k_2=0}^{n_2}{p}_{n_1-1,k_1}\left(x\right)·p_{n_2,k_2}\left(y\right)\left(t_1(k_1+1)-t_1\left(k_1\right)\right)\left[\begin{array}{ccc}x,& t_1\left(k_1\right),& t_1(k_1+1)\\ & t_2\left(k_2\right)\end{array};F\right]\\ =-\frac{(x-a)(b-x)}{n_1}\sum _{k_1=0}^{n_1-1}\sum _{k_2=0}^{n_2}{p}_{n_1-1,k_1}\left(x\right)p_{n_2,k_2}\left(y\right)\left[\begin{array}{ccc}x,& t_1\left(k_1\right),& t_1(k_1+1)\\ & t_2\left(k_2\right)\end{array};F\right].\end{array}$$

In a similar way, one obtains

$$S_2=-\frac{(y-c)(d-y)}{n_2}\sum _{k_1=0}^{n_1}\sum _{k_2=0}^{n_2-1}{p}_{n_1,k_1}\left(x\right)p_{n_2-1,k_2}\left(y\right)\left[\begin{array}{ccc}& t_1\left(k_1\right)& \\ y,& t_2\left(k_2\right),& t_2(k_2+1)\end{array};F\right]$$

and

$$S_3=\frac{(x-a)(b-x)(y-c)(d-y)}{n_1·n_2}\sum _{k_1=0}^{n_1-1}\sum _{k_2=0}^{n_2-1}{p}_{n_1-1,k_1}\left(x\right)p_{n_2-1,k_2}\left(y\right)\left[\begin{array}{ccc}x,& t_1\left(k_1\right),& t_1(k_1+1)\\ y,& t_2\left(k_2\right),& t_2(k_2+1)\end{array};F\right].$$

□

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)$.

Corollary 1.

If $F\in C^{2,2}\left(S\right)$, then

$$\begin{array}{cc}\hfill R_{n_1,n_2}(F;x,y)=& -\frac{(x-a)(b-x)}{2n_1}{F}^{(2,0)}(\xi ,y)-\frac{(y-c)(d-y)}{2n_2}{F}^{(0,2)}(x,\eta )\hfill \\ & +\frac{(x-a)(b-x)(y-c)(d-y)}{4n_1{n}_2}{F}^{(2,2)}(\xi ,\eta ),\hfill \end{array}$$

for $\xi ,\eta \in (a,b)\times (c,d)$ and

$$\begin{array}{cc}\hfill \left|R_{n_1,n_2}(F;x,y)\right|& \le \frac{(x-a)(b-x)}{2n_1}{∥F^{(2,0)}∥}_{\infty }+\frac{(y-c)(d-y)}{2n_2}{∥F^{(0,2)}∥}_{\infty }\hfill \\ & +\frac{(x-a)(b-x)(y-c)(d-y)}{4n_1{n}_2}{∥F^{(2,2)}∥}_{\infty }\hfill \\ & \le \frac{{(b-a)}^2}{8n_1}{∥F^{(2,0)}∥}_{\infty }+\frac{{(d-c)}^2}{8n_2}{∥F^{(0,2)}∥}_{\infty }+\frac{{(b-a)}^2{(d-c)}^2}{64n_1{n}_2}{∥F^{(2,2)}∥}_{\infty }.\hfill \end{array}$$

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) on $S=[a,b]\times [c,d]$, such that

$$\underset{a}{\overset{b}{\int }}\underset{c}{\overset{d}{\int }}F(x,y)dxdy=\underset{a}{\overset{b}{\int }}\underset{c}{\overset{d}{\int }}{B}_{n_1,n_2}(F;x,y)dxdy+\underset{a}{\overset{b}{\int }}\underset{c}{\overset{d}{\int }}{R}_{n_1,n_2}(F;x,y)dxdy.$$

(8)

Theorem 2.

The approximation of the volume obtained by integration of any real-valued function $F:S\to \mathbb{R}$ can be done with the classical Bernstein cubature formula

$$\underset{a}{\overset{b}{\int }}\underset{c}{\overset{d}{\int }}F(x,y)dxdy\approx \frac{(b-a)(d-c)}{(n_1+1)(n_2+1)}\sum _{k_1=0}^{n_1}\sum _{k_2=0}^{n_2}F\left(a+\frac{k_1(b-a)}{n_1},c+\frac{k_2(d-c)}{n_2}\right).$$

(9)

Proof. 

Using the definition of the classical bivariate Bernstein operators (6), we get

$$\begin{array}{c}\underset{a}{\overset{b}{\int }}\underset{c}{\overset{d}{\int }}{B}_{n_1,n_2}(F;x,y)dxdy=\underset{a}{\overset{b}{\int }}\underset{c}{\overset{d}{\int }}\sum _{k_1=0}^{n_1}\sum _{k_2=0}^{n_2}{p}_{n_1,k_1}\left(x\right)·p_{n_2,k_2}\left(y\right)F\left(t_1\left(k_1\right),t_2\left(k_2\right)\right)dxdy\\ =\sum _{k_1=0}^{n_1}\sum _{k_2=0}^{n_2}F\left(a+\frac{k_1(b-a)}{n_1},c+\frac{k_2(d-c)}{n_2}\right)·{\int }_a^b{p}_{n_1,k_1}\left(x\right)dx·{\int }_c^d{p}_{n_2,k_2}\left(y\right)dy.\end{array}$$

Taking into account the definition and the properties of the Euler’s integral of first kind $\beta (u,v)={\int }_0^1{t}^{u-1}{(1-t)}^{v-1}dt$, it can be shown that

$$\begin{array}{cc}\hfill {\int }_a^b{p}_{n_1,k_1}\left(x\right)dx& =\frac{1}{{(b-a)}^{n_1}}\left(\genfrac{}{}{0pt}{}{n_1}{k_1}\right){\int }_a^b{(x-a)}^{k_1}{(b-x)}^{n_1-k_1}dx\hfill \\ & =\frac{1}{{(b-a)}^{n_1}}\left(\genfrac{}{}{0pt}{}{n_1}{k_1}\right){\int }_0^1{t}^{k_1}{(1-t)}^{n_1-k_1}{(b-a)}^{n_1+1}dt\hfill \\ & =\left(\genfrac{}{}{0pt}{}{n_1}{k_1}\right)(b-a)·\beta (k_1+1,n_1+1-k_1)=\frac{b-a}{n_1+1}\left(\genfrac{}{}{0pt}{}{n_1}{k_1}\right){\left(\genfrac{}{}{0pt}{}{n_1}{k_1}\right)}^{-1}=\frac{b-a}{n_1+1}\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill {\int }_c^d{p}_{n_2,k_2}\left(y\right)dy& =\frac{1}{{(d-c)}^{n_2}}\left(\genfrac{}{}{0pt}{}{n_2}{k_2}\right){\int }_c^d{(y-c)}^{k_2}{(d-y)}^{n_2-k_2}dy\hfill \\ & =\frac{1}{{(d-c)}^{n_2}}\left(\genfrac{}{}{0pt}{}{n_2}{k_2}\right){\int }_0^1{t}^{k_2}{(1-t)}^{n_2-k_2}{(d-c)}^{n_2+1}dt\hfill \\ & =\left(\genfrac{}{}{0pt}{}{n_2}{k_2}\right)(d-c)·\beta (k_2+1,n_2+1-k_2)=\frac{d-c}{n_2+1}\left(\genfrac{}{}{0pt}{}{n_2}{k_2}\right){\left(\genfrac{}{}{0pt}{}{n_2}{k_2}\right)}^{-1}=\frac{d-c}{n_2+1}.\hfill \end{array}$$

Collecting the results demonstrated above and eliminating the remainder term in the relation (8), we get the classical Bernstein cubature formula

$$\begin{array}{cc}\hfill \underset{a}{\overset{b}{\int }}\underset{c}{\overset{d}{\int }}F(x,y)dxdy& \approx \underset{a}{\overset{b}{\int }}\underset{c}{\overset{d}{\int }}{B}_{n_1,n_2}(F;x,y)dxdy\hfill \\ & =\frac{(b-a)(d-c)}{(n_1+1)(n_2+1)}\sum _{k_1=0}^{n_1}\sum _{k_2=0}^{n_2}F\left(a+\frac{k_1(b-a)}{n_1},c+\frac{k_2(d-c)}{n_2}\right).\hfill \end{array}$$

□

Theorem 3.

Assuming that $F\in C^{2,2}\left(S\right)$ and $\xi ,\eta \in (a,b)\times (c,d)$, then

$$\begin{array}{cc}\hfill R_{n_1,n_2}\left[F\right]& =\underset{a}{\overset{b}{\int }}\underset{c}{\overset{d}{\int }}{R}_{n_1,n_2}(F;x,y)dxdy\hfill \\ & =-\frac{{(b-a)}^3(d-c)}{12n_1}{F}^{(2,0)}(\xi ,\eta )-\frac{(b-a){(d-c)}^3}{12n_2}{F}^{(0,2)}(x,\eta )+\frac{{(b-a)}^3{(d-c)}^3}{144n_1{n}_2}{F}^{(2,2)}(\xi ,\eta )\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill \left|R_{n_1,n_2}\left[F\right]\right|& =\left|\underset{a}{\overset{b}{\int }}\underset{c}{\overset{d}{\int }}{R}_{n_1,n_2}(F;x,y)dxdy\right|\hfill \\ & \le \frac{{(b-a)}^3(d-c)}{12n_1}{∥F^{(2,0)}∥}_{\infty }+\frac{(b-a){(d-c)}^3}{12n_2}{∥F^{(0,2)}∥}_{\infty }+\frac{{(b-a)}^3{(d-c)}^3}{144n_1{n}_2}{∥F^{(2,2)}∥}_{\infty }.\hfill \end{array}$$

Proof. 

Both forms of the remainder term presented in Corollary 1 get integrated on S. □

Summarizing the above results, we get the classical bivariate Bernstein formula

$$\underset{a}{\overset{b}{\int }}\underset{c}{\overset{d}{\int }}F(x,y)dxdy=\frac{(b-a)(d-c)}{(n_1+1)(n_2+1)}\sum _{k_1=0}^{n_1}\sum _{k_2=0}^{n_2}F\left(a+\frac{k_1(b-a)}{n_1},c+\frac{k_2(d-c)}{n_2}\right)+R_{n_1,n_2}\left[F\right],$$

(10)

with the remainder term (assuming that $F\in C^{2,2}\left(S\right)$ and $\xi ,\eta \in (a,b)\times (c,d)$)

$$R_{n_1,n_2}\left[F\right]=-\frac{{(b-a)}^3(d-c)}{12n_1}{F}^{(2,0)}(\xi ,y)-\frac{(b-a){(d-c)}^3}{12n_2}{F}^{(0,2)}(x,\eta )+\frac{{(b-a)}^3{(d-c)}^3}{144n_1{n}_2}{F}^{(2,2)}(\xi ,\eta ).$$

Remark 1.

Due to the form of the remainder term, the classical bivariate Bernstein formula (10) has the degree of exactness $(1,1)$. The statement could be checked in a series of examples.

Example 1.

Let us consider the following double integral

$${\int }_{-1}^3{\int }_2^4(7x+5y)dxdy.$$

Applying the classical bivariate Bernstein formula (10) with $n_1=n_2=1$ (the simplest case), we get

$${\int }_{-1}^3{\int }_2^4(7x+5y)dxdy=\frac{(3+1)·(4-2)}{2·2}\sum _{k_1=0}^1\sum _{k_2=0}^{1}F(-1+4k_1,2+2k_2)=176,$$

which is the exact value of the given integral.

Below, we present some particular cases of the classical bivariate Bernstein formula.

Case 1.

If $S=[0,1]\times [0,1]$, we get the bivariate Bernstein formula [22]

$${\int }_0^1{\int }_0^{1}F(x,y)dxdy=\frac{1}{(n_1+1)(n_2+1)}\sum _{k_1=0}^{n_1}\sum _{k_2=0}^{n_2}F\left(\frac{k_1}{n_1},\frac{k_2}{n_2}\right)+R_{n_1,n_2}\left[F\right],$$

with

$$R_{n_1,n_2}\left[F\right]=-\frac{1}{12n_1}{F}^{(2,0)}(\xi ,y)-\frac{1}{12n_2}{F}^{(0,2)}(x,\eta )+\frac{1}{144n_1{n}_2}{F}^{(2,2)}(\xi ,\eta ),$$

for $\xi ,\eta \in (0,1)\times (0,1).$ The upper bound estimation is

$$\left|R_{n_1,n_2}\left[F\right]\right|\le \frac{1}{12n_1}{∥F^{(2,0)}∥}_{\infty }+\frac{1}{12n_2}{∥F^{(0,2)}∥}_{\infty }+\frac{1}{144n_1{n}_2}{∥F^{(2,2)}∥}_{\infty }.$$

Case 2.

If we consider $n_1=n_2=1$ in relation (10), then we get the bivariate trapezoidal formula

$${\int }_a^b{\int }_c^{d}F(x,y)dxdy=\frac{(b-a)(d-c)}{4}\left(F(a,c)+F(a,d)+F(b,c)+F(b,d)\right)+R_{1,1}\left[F\right],$$

with

$$R_{1,1}\left[F\right]=-\frac{{(b-a)}^3(d-c)}{12}{F}^{(2,0)}(\xi ,y)-\frac{(b-a){(d-c)}^3}{12}{F}^{(0,2)}(x,\eta )+\frac{{(b-a)}^3{(d-c)}^3}{144}{F}^{(2,2)}(\xi ,\eta ).$$

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

$$\begin{array}{cc}\hfill B_{1,1}(F;x,y)& =\sum _{k_1=0}^1\sum _{k_2=0}^1{p}_{1,k_1}\left(x\right)p_{1,k_2}\left(y\right)·F\left(a+k_1(b-a),c+k_2(d-c)\right)\hfill \\ & =\frac{(b-x)(d-y)}{(b-a)(d-c)}F(a,c)+\frac{(b-x)(y-c)}{(b-a)(d-c)}F(a,d)+\frac{(x-a)(d-y)}{(b-a)(d-c)}F(b,c)\hfill \\ & +\frac{(x-a)(y-c)}{(b-a)(d-c)}F(b,d)=L_{1,1}(F;x,y).\hfill \end{array}$$

To demonstrate the approximation accuracy of the classical Bernstein cubature formula

$$\underset{a}{\overset{b}{\int }}\underset{c}{\overset{d}{\int }}F(x,y)dxdy\approx \frac{(b-a)(d-c)}{(n_1+1)(n_2+1)}\sum _{k_1=0}^{n_1}\sum _{k_2=0}^{n_2}F\left(a+\frac{k_1(b-a)}{n_1},c+\frac{k_2(d-c)}{n_2}\right)=:I_{n_1,n_2}\left(F\right),$$

(11)

we test it with the Maple mathematical software. For our study, we consider five test functions as follows:

$$\begin{array}{cc}\hfill F_1:[0,0.75]\times [0,0.75]\to \mathbb{R},F_1(x,y)& =e^{2y-x};\hfill \\ \hfill F_2:[1,2]\times [1,2]\to \mathbb{R},F_2(x,y)& =5x^2+3x{y}^2+7y;\hfill \\ \hfill F_3:[1.4,2]\times [1,1.5]\to \mathbb{R},F_3(x,y)& =ln(x+2y);\hfill \\ \hfill F_4:[-1,1]\times [-1,1]\to \mathbb{R},F_4(x,y)& =e^{-(x^2+y^2)};\hfill \\ \hfill F_5:[0,4]\times [0,3]\to \mathbb{R},F_5(x,y)& =e^{-(x+y)}·sin(2x+2y).\hfill \end{array}$$

Let us denote the approximative values of the integrals with ten exact decimals by

$$\begin{array}{c}{I}_{F_1}=\underset{0}{\overset{0.75}{\int }}\underset{0}{\overset{0.75}{\int }}{F}_1(x,y)dxdy,I_{F_2}=\underset{1}{\overset{2}{\int }}\underset{1}{\overset{2}{\int }}{F}_2(x,y)dxdy,I_{F_3}=\underset{1.4}{\overset{2}{\int }}\underset{1}{\overset{1.5}{\int }}{F}_3(x,y)dxdy,\\ I_{F_4}=\underset{-1}{\overset{1}{\int }}\underset{-1}{\overset{1}{\int }}{F}_4(x,y)dxdy,I_{F_5}=\underset{0}{\overset{4}{\int }}\underset{0}{\overset{3}{\int }}{F}_5(x,y)dxdy,\end{array}$$

and the appropriate absolute errors

$$e{I}_{n_1,n_2}\left(F\right)=|I_F-I_{n_1,n_2}\left(F\right)|.$$

We approximate the double integrals and record in

Table 1. The approximation of the given double integrals.

${\mathit{I}}_{\mathit{F}}$	${\mathit{eI}}_{1,1}\left(\mathit{F}\right)$	${\mathit{eI}}_{10,10}\left(\mathit{F}\right)$	${\mathit{eI}}_{50,50}\left(\mathit{F}\right)$	${\mathit{eI}}_{500,500}\left(\mathit{F}\right)$
$I_{F_1}$	$2.164\times {10}^{-1}$	$2.101\times {10}^{-2}$	$4.180\times {10}^{-3}$	$4.175\times {10}^{-4}$
$I_{F_2}$	$1.583$	$1.583\times {10}^{-1}$	$3.166\times {10}^{-2}$	$3.166\times {10}^{-3}$
$I_{F_3}$	$1.971\times {10}^{-3}$	$1.958\times {10}^{-4}$	$3.913\times {10}^{-5}$	$3.911\times {10}^{-6}$
$I_{F_4}$	$1.689$	$2.137\times {10}^{-1}$	$4.474\times {10}^{-2}$	$4.522\times {10}^{-3}$
$I_{F_5}$	$1.378\times {10}^{-1}$	$8.252\times {10}^{-2}$	$2.076\times {10}^{-2}$	$2.162\times {10}^{-3}$

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}\left(S\right)$;
• A double integral could be approximated with a desired precision $\epsilon$, where $\epsilon >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.

4. 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. With this aim, it is advisable to divide $[a,b]\times [c,d]$ into smaller bivariate intervals and to use the classical Bernstein cubature formula on each bivariate subinterval. By adding the results, we get the approximation of the double integral on the whole interval $[a,b]\times [c,d]$. Such a scheme is called the classical composite Bernstein cubature formula. Now, let us consider the bivariate interval $[a,b]\times [c,d]$ divided in $m_1·m_2$ equally spaced subintervals, such that, for $i=\overline{1,m_1}$ and $j=\overline{1,m_2}$, we get

$$\left[a+\frac{(i-1)(b-a)}{m_1},a+\frac{i(b-a)}{m_1}\right]\times \left[c+\frac{(j-1)(d-c)}{m_2},c+\frac{j(d-c)}{m_2}\right]=:I_i\times I_j.$$

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\in \mathbb{N}$, we consider the parametric extensions

$${}_x{B}_{n_1,m_1}(F;x,y)=\frac{m_1^{n_1}}{{(b-a)}^{n_1}}\sum _{k_1=0}^{n_1}\left(\genfrac{}{}{0pt}{}{n_1}{k_1}\right){\left(x-a-\frac{(i-1)(b-a)}{m_1}\right)}^{k_1}{\left(a+\frac{i(b-a)}{m_1}-x\right)}^{n_1-k_1}F\left(a+\frac{(n_{1}i-n_1+k_1)(b-a)}{m_1·n_1},y\right)$$

and

$${}_y{B}_{n_2,m_2}(F;x,y)=\frac{m_2^{n_2}}{{(d-c)}^{n_2}}\sum _{k_2=0}^{n_2}\left(\genfrac{}{}{0pt}{}{n_2}{k_2}\right){\left(y-c-\frac{(j-1)(d-c)}{m_2}\right)}^{k_2}{\left(c+\frac{j(d-c)}{m_2}-y\right)}^{n_2-k_2}F\left(x,c+\frac{(n_{2}j-n_2+k_2)(d-c)}{m_2·n_2}\right).$$

The tensor product of the parametric extensions defined above is $\left({}_x{B}_{n_1,m_1}\otimes {}_y{B}_{n_2,m_2}\right)$$(F;x,y)=:\mathcal{B}(F;x,y)$, having the following expression:

$$\begin{array}{cc}\hfill \mathcal{B}(F;x,y)=& \frac{m_1^{n_1}·m_2^{n_2}}{{(b-a)}^{n_1}·{(d-c)}^{n_2}}\sum _{k_1=0}^{n_1}\sum _{k_2=0}^{n_2}\left(\genfrac{}{}{0pt}{}{n_1}{k_1}\right)\left(\genfrac{}{}{0pt}{}{n_2}{k_2}\right){\left(x-a-\frac{(i-1)(b-a)}{m_1}\right)}^{k_1}{\left(a+\frac{i(b-a)}{m_1}-x\right)}^{n_1-k_1}\hfill \\ & ·{\left(y-c-\frac{(j-1)(d-c)}{m_2}\right)}^{k_2}{\left(c+\frac{j(d-c)}{m_2}-y\right)}^{n_2-k_2}F\left(a+\frac{(n_{1}i-n_1+k_1)(b-a)}{m_1·n_1},c+\frac{(n_{2}j-n_2+k_2)(d-c)}{m_2·n_2}\right).\hfill \end{array}$$

The classical bivariate Bernstein approximation formula on the interval $I_i\times I_j$, $i=\overline{1,m_1}$$\mathrm{and}j=\overline{1,m_2}$ becomes

$$F(x,y)=\mathcal{B}(F;x,y)+\mathcal{R}(F;x,y),$$

(12)

with $\mathcal{R}\left(F\right)$ the bivariate remainder operators associated with the tensor product $\mathcal{B}\left(F\right)$. Taking Corollary 1 into account and supposing that $F\in C^{2,2}\left(I_i\times I_j\right)$, $i=\overline{1,m_1},j=\overline{1,m_2}$, we can present an upper bound estimation of the remainder term given by

$$\begin{array}{cc}\hfill \left|\mathcal{R}\right(F;x,y\left)\right|\le & \frac{\left(x-a-\frac{(i-1)(b-a)}{m_1}\right)\left(a+\frac{i(b-a)}{m_1}-x\right)}{2n_1}{∥F^{(2,0)}∥}_{\infty }+\frac{\left(y-c-\frac{(j-1)(d-c)}{m_2}\right)\left(c+\frac{j(d-c)}{m_2}-y\right)}{2n_2}{∥F^{(0,2)}∥}_{\infty }\hfill \\ & +\frac{\left(x-a-\frac{(i-1)(b-a)}{m_1}\right)\left(a+\frac{i(b-a)}{m_1}-x\right)\left(y-c-\frac{(j-1)(d-c)}{m_2}\right)\left(c+\frac{j(d-c)}{m_2}-y\right)}{4n_1{n}_2}{∥F^{(2,2)}∥}_{\infty }.\hfill \end{array}$$

(13)

Integrating the classical bivariate Bernstein approximation formula (12) on $I_i\times I_j$, $i=\overline{1,m_1}$, and $j=\overline{1,m_2}$, we get

$$\underset{I_i\times I_j}{\int \int }F(x,y)dxdy=\underset{I_i\times I_j}{\int \int }\mathcal{B}(F;x,y)dxdy+\underset{I_i\times I_j}{\int \int }\mathcal{R}(F;x,y)dxdy.$$

(14)

Summing the relations (14) for $i=\overline{1,m_1}$ and $j=\overline{1,m_2}$, one obtains

$$\underset{a}{\overset{b}{\int }}\underset{c}{\overset{d}{\int }}F(x,y)dxdy=\sum _{i=1}^{m_1}\sum _{j=1}^{m_2}\underset{I_i\times I_j}{\int \int }\mathcal{B}(F;x,y)dxdy+\sum _{i=1}^{m_1}\sum _{j=1}^{m_2}\underset{I_i\times I_j}{\int \int }\mathcal{R}(F;x,y)dxdy.$$

(15)

Theorem 4.

The definite double integral of any real-valued function $F:[a,b]\times [c,d]\to \mathbb{R}$ could be approximated with the formula

$$\underset{a}{\overset{b}{\int }}\underset{c}{\overset{d}{\int }}F(x,y)dxdy\approx \frac{(b-a)(d-c)}{m_1·m_2(n_1+1)(n_2+1)}\sum _{i=1}^{m_1}\sum _{j=1}^{m_2}\sum _{k_1=0}^{n_1}\sum _{k_2=0}^{n_2}F\left(a+\frac{(n_{1}i-n_1+k_1)(b-a)}{m_1·n_1},c+\frac{(n_{2}j-n_2+k_2)(d-c)}{m_2·n_2}\right).$$

(16)

Proof. 

Let $i=\overline{1,m_1}$ and $j=\overline{1,m_2}$ be fixed. Using the definition of the tensor product $\mathcal{B}\left(F\right)$, we have

$$\begin{array}{c}\underset{I_i\times I_j}{\int \int }\mathcal{B}(F;x,y)dxdy=\frac{m_1^{n_1}·m_2^{n_2}}{{(b-a)}^{n_1}·{(d-c)}^{n_2}}\sum _{k_1=0}^{n_1}\sum _{k_2=0}^{n_2}\left(\genfrac{}{}{0pt}{}{n_1}{k_1}\right)\left(\genfrac{}{}{0pt}{}{n_2}{k_2}\right)·F(\star ,\circ )\\ ·\underset{I_i\times I_j}{\int \int }{\left(x-a-\frac{(i-1)(b-a)}{m_1}\right)}^{k_1}{\left(a+\frac{i(b-a)}{m_1}-x\right)}^{n_1-k_1}·{\left(y-c-\frac{(j-1)(d-c)}{m_2}\right)}^{k_2}{\left(c+\frac{j(d-c)}{m_2}-y\right)}^{n_2-k_2}dxdy\\ =\frac{m_1^{n_1}·m_2^{n_2}}{{(b-a)}^{n_1}·{(d-c)}^{n_2}}\sum _{k_1=0}^{n_1}\sum _{k_2=0}^{n_2}\left(\genfrac{}{}{0pt}{}{n_1}{k_1}\right)\left(\genfrac{}{}{0pt}{}{n_2}{k_2}\right)·F(\star ,\circ )·{\int }_0^1\frac{t^{k_1}{(1-t)}^{n_1-k_1}{(b-a)}^{n_1+1}}{m_1^{n_1+1}}dt·{\int }_0^1\frac{r^{k_2}{(1-r)}^{n_2-k_2}{(d-c)}^{n_2+1}}{m_2^{n_2+1}}dr\\ =\frac{(b-a)(d-c)}{m_1·m_2}\sum _{k_1=0}^{n_1}\sum _{k_2=0}^{n_2}\left(\genfrac{}{}{0pt}{}{n_1}{k_1}\right)\left(\genfrac{}{}{0pt}{}{n_2}{k_2}\right)·F(\star ,\circ )·\beta (k_1+1,n_1+1-k_1)·\beta (k_2+1,n_2+1-k_2)\\ =\frac{(b-a)(d-c)}{m_1·m_2(n_1+1)(n_2+1)}\sum _{k_1=0}^{n_1}\sum _{k_2=0}^{n_2}\left(\genfrac{}{}{0pt}{}{n_1}{k_1}\right){\left(\genfrac{}{}{0pt}{}{n_1}{k_1}\right)}^{-1}\left(\genfrac{}{}{0pt}{}{n_2}{k_2}\right){\left(\genfrac{}{}{0pt}{}{n_2}{k_2}\right)}^{-1}·F(\star ,\circ )=\frac{(b-a)(d-c)}{m_1·m_2(n_1+1)(n_2+1)}\sum _{k_1=0}^{n_1}\sum _{k_2=0}^{n_2}F(\star ,\circ ).\end{array}$$

Summing after all values of $i=\overline{1,m_1}$ and $j=\overline{1,m_2}$, we get

$$\underset{a}{\overset{b}{\int }}\underset{c}{\overset{d}{\int }}F(x,y)dxdy\approx \frac{(b-a)(d-c)}{m_1{m}_2(n_1+1)(n_2+1)}\sum _{i=1}^{m_1}\sum _{j=1}^{m_2}\sum _{k_1=0}^{n_1}\sum _{k_2=0}^{n_2}F\left(a+\frac{(n_{1}i-n_1+k_1)(b-a)}{m_1·n_1},c+\frac{(n_{2}j-n_2+k_2)(d-c)}{m_2·n_2}\right).$$

□

For brevity, we designate

$$\mathcal{R}\left[F\right]:=\sum _{i=1}^{m_1}\sum _{j=1}^{m_2}\underset{I_i\times I_j}{\int \int }\mathcal{R}(F;x,y)dxdy={\int }_a^b{\int }_c^d\mathcal{R}(F;x,y)dxdy.$$

Theorem 5.

If $F\in C^{2,2}\left(S\right)$, then

$$\left|\mathcal{R}\left[F\right]\right|\le \frac{{(b-a)}^3(d-c)}{12m_1^2{n}_1}{∥F^{(2,0)}∥}_{\infty }+\frac{(b-a){(d-c)}^3}{12m_2^2{n}_2}{∥F^{(0,2)}∥}_{\infty }+\frac{{(b-a)}^3{(d-c)}^3}{144m_1^2{m}_2^2{n}_1{n}_2}{∥F^{(2,2)}∥}_{\infty }.$$

(17)

Proof. 

Starting with the relation (15), we have

$$\begin{array}{cc}\hfill \left|\mathcal{R}\right[F\left]\right|& =\left|{\int }_a^b{\int }_c^{d}F(x,y)dxdy-{\int }_a^b{\int }_c^d\mathcal{B}(F;x,y)dxdy\right|\hfill \\ & =\left|\sum _{i=1}^{m_1}\sum _{j=1}^{m_2}\underset{I_i\times I_j}{\int \int }F(x,y)dxdy-\sum _{i=1}^{m_1}\sum _{j=1}^{m_2}\underset{I_i\times I_j}{\int \int }\mathcal{B}(F;x,y)dxdy\right|\hfill \\ & \le \sum _{i=1}^{m_1}\sum _{j=1}^{m_2}\underset{I_i\times I_j}{\int \int }\left|\mathcal{R}(F;x,y)\right|dxdy.\hfill \end{array}$$

Using the relation (13), it follows that

$$\begin{array}{c}\left|\mathcal{R}\left[F\right]\right|\le \sum _{i=1}^{m_1}\sum _{j=1}^{m_2}\underset{I_i\times I_j}{\int \int }\frac{\left(x-a-\frac{(i-1)(b-a)}{m_1}\right)\left(a+\frac{i(b-a)}{m_1}-x\right)}{2n_1}{∥F^{(2,0)}∥}_{\infty }dxdy\end{array}$$

$$\begin{array}{c}+\sum _{i=1}^{m_1}\sum _{j=1}^{m_2}\underset{I_i\times I_j}{\int \int }\frac{\left(y-c-\frac{(j-1)(d-c)}{m_2}\right)\left(c+\frac{j(d-c)}{m_2}-y\right)}{2n_2}{∥F^{(0,2)}∥}_{\infty }dxdy\\ +\sum _{i=1}^{m_1}\sum _{j=1}^{m_2}\underset{I_i\times I_j}{\int \int }\frac{\left(x-a-\frac{(i-1)(b-a)}{m_1}\right)\left(a+\frac{i(b-a)}{m_1}-x\right)\left(y-c-\frac{(j-1)(d-c)}{m_2}\right)\left(c+\frac{j(d-c)}{m_2}-y\right)}{4n_1{n}_2}{∥F^{(2,2)}∥}_{\infty }dxdy\\ =\sum _{i=1}^{m_1}\sum _{j=1}^{m_2}\frac{1}{2n_1}{∥F^{(2,0)}∥}_{\infty }{\int }_0^{1}t(1-t)\frac{{(b-a)}^3}{m_1^3}dt\underset{c+\frac{(j-1)(d-c)}{m_2}}{\overset{c+\frac{j(d-c)}{m_2}}{\int }}dy+\sum _{i=1}^{m_1}\sum _{j=1}^{m_2}\frac{1}{2n_2}{∥F^{(0,2)}∥}_{\infty }\underset{a+\frac{(i-1)(b-a)}{m_1}}{\overset{a+\frac{i(b-a)}{m_1}}{\int }}dx{\int }_0^{1}r(1-r)\frac{{(d-c)}^3}{m_2^3}dr\\ +\sum _{i=1}^{m_1}\sum _{j=1}^{m_2}\frac{1}{4n_1{n}_2}{∥F^{(2,2)}∥}_{\infty }{\int }_0^{1}t(1-t)\frac{{(b-a)}^3}{m_1^3}dt{\int }_0^{1}r(1-r)\frac{{(d-c)}^3}{m_2^3}dr\\ =\frac{{(b-a)}^3(d-c)}{12m_1^2{n}_1}{∥F^{(2,0)}∥}_{\infty }+\frac{(b-a){(d-c)}^3}{12m_2^2{n}_2}{∥F^{(0,2)}∥}_{\infty }+\frac{{(b-a)}^3{(d-c)}^3}{144m_1^2{m}_2^2{n}_1{n}_2}{∥F^{(2,2)}∥}_{\infty }.\end{array}$$

□

Summarizing the above results, we get the classical bivariate Bernstein composite formula

$$\underset{a}{\overset{b}{\int }}\underset{c}{\overset{d}{\int }}F(x,y)dxdy=\frac{(b-a)(d-c)}{m_1·m_2(n_1+1)(n_2+1)}\sum _{i=1}^{m_1}\sum _{j=1}^{m_2}\sum _{k_1=0}^{n_1}\sum _{k_2=0}^{n_2}F\left(\star ,\circ \right)+\mathcal{R}\left[F\right],$$

(18)

with an upper bound estimation of the remainder term (assuming $F\in C^{2,2}\left(S\right)$)

$$\left|\mathcal{R}\left[F\right]\right|\le \frac{{(b-a)}^3(d-c)}{12m_1^2{n}_1}{∥F^{(2,0)}∥}_{\infty }+\frac{(b-a){(d-c)}^3}{12m_2^2{n}_2}{∥F^{(0,2)}∥}_{\infty }+\frac{{(b-a)}^3{(d-c)}^3}{144m_1^2{m}_2^2{n}_1{n}_2}{∥F^{(2,2)}∥}_{\infty }.$$

Below, we present some particular cases of the classical bivariate Bernstein composite formula.

Case 3.

If $S=[0,1]\times [0,1]$, then we get the bivariate Bernstein composite formula [2]

$${\int }_0^1{\int }_0^{1}F(x,y)dxdy=\frac{1}{m_1·m_2(n_1+1)(n_2+1)}\sum _{i=1}^{m_1}\sum _{j=1}^{m_2}\sum _{k_1=0}^{n_1}\sum _{k_2=0}^{n_2}F\left(\frac{n_{1}i-n_1+k_1}{m_1·n_1},\frac{n_{2}j-n_2+k_2}{m_2·n_2}\right)+\mathcal{R}\left[F\right],$$

with

$$\left|\mathcal{R}\left[F\right]\right|\le \frac{1}{12m_1^2{n}_1}{∥F^{(2,0)}∥}_{\infty }+\frac{1}{12m_2^2{n}_2}{∥F^{(0,2)}∥}_{\infty }+\frac{1}{144m_1^2{m}_2^2{n}_1{n}_2}{∥F^{(2,2)}∥}_{\infty }.$$

(19)

Remark 5.

The corrected upper bound estimation (19) of the remainder term in the bivariate Bernstein composite formula was established in [1].

Case 4.

If we consider $n_1=n_2=1$ in the relation (18), then we get the bivariate trapezoidal composite formula

$$\begin{array}{cc}\hfill {\int }_a^b{\int }_c^{d}F(x,y)dxdy=& \frac{(b-a)(d-c)}{4m_1·m_2}\sum _{i=1}^{m_1}\sum _{j=1}^{m_2}\left(F\left(a+\frac{(i-1)(b-a)}{m_1},c+\frac{(j-1)(d-c)}{m_2}\right)+F\left(a+\frac{(i-1)(b-a)}{m_1},c+\frac{j(d-c)}{m_2}\right)\right.\hfill \\ & \left.+F\left(a+\frac{i(b-a)}{m_1},c+\frac{(j-1)(d-c)}{m_2}\right)+F\left(a+\frac{i(b-a)}{m_1},c+\frac{j(d-c)}{m_2}\right)\right)+\mathcal{R}\left[F\right],\hfill \end{array}$$

with

$$\left|\mathcal{R}\left[F\right]\right|\le \frac{{(b-a)}^3(d-c)}{12m_1^2}{∥F^{(2,0)}∥}_{\infty }+\frac{(b-a){(d-c)}^3}{12m_2^2}{∥F^{(0,2)}∥}_{\infty }+\frac{{(b-a)}^3{(d-c)}^3}{144m_1^2{m}_2^2}{∥F^{(2,2)}∥}_{\infty }.$$

5. Numerical Examples

In order to validate the approximation accuracy of the classical composite Bernstein cubature formula

$$\underset{a}{\overset{b}{\int }}\underset{c}{\overset{d}{\int }}F(x,y)dxdy\approx \frac{(b-a)(d-c)}{m_1{m}_2(n_1+1)(n_2+1)}\sum _{i=1}^{m_1}\sum _{j=1}^{m_2}\sum _{k_1=0}^{n_1}\sum _{k_2=0}^{n_2}F\left(\star ,\circ \right)=:I_{m_1,m_2,n_1,n_2}\left(F\right),$$

(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

$$e{I}_{m_1,m_2,n_1,n_2}\left(F\right)=|I_F-I_{m_1,m_2,n_1,n_2}\left(F\right)|,$$

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 (

Table 2. The approximation of the double integral ${\int }_0^{0.75}{\int }_0^{0.75}{e}^{2y-x}dxdy$.

${\mathit{m}}_1$	${\mathit{m}}_2$	${\mathit{eI}}_{{\mathit{m}}_1,{\mathit{m}}_2,1,1}\left(\mathit{F}\right)$	${\mathit{eI}}_{{\mathit{m}}_1,{\mathit{m}}_2,5,10}\left(\mathit{F}\right)$	${\mathit{eI}}_{{\mathit{m}}_1,{\mathit{m}}_2,10,5}\left(\mathit{F}\right)$	${\mathit{eI}}_{{\mathit{m}}_1,{\mathit{m}}_2,50,50}\left(\mathit{F}\right)$
2	2	$5.389\times {10}^{-2}$	$6.427\times {10}^{-3}$	$9.628\times {10}^{-3}$	$1.068\times {10}^{-3}$
5	10	$3.446\times {10}^{-3}$	$5.165\times {10}^{-4}$	$5.165\times {10}^{-4}$	$6.886\times {10}^{-5}$
10	5	$7.312\times {10}^{-3}$	$7.741\times {10}^{-4}$	$1.419\times {10}^{-3}$	$1.461\times {10}^{-4}$
25	25	$3.444\times {10}^{-4}$	$4.133\times {10}^{-5}$	$6.199\times {10}^{-5}$	$6.884\times {10}^{-6}$
50	50	$8.611\times {10}^{-5}$	$1.033\times {10}^{-5}$	$1.550\times {10}^{-5}$	$1.718\times {10}^{-6}$

,

Table 3. The approximation of the double integral ${\int }_1^2{\int }_1^2\left(5x^2+3x{y}^2+7y\right)dxdy$.

${\mathit{m}}_1$	${\mathit{m}}_2$	${\mathit{eI}}_{{\mathit{m}}_1,{\mathit{m}}_2,1,1}\left(\mathit{F}\right)$	${\mathit{eI}}_{{\mathit{m}}_1,{\mathit{m}}_2,5,10}\left(\mathit{F}\right)$	${\mathit{eI}}_{{\mathit{m}}_1,{\mathit{m}}_2,10,5}\left(\mathit{F}\right)$	${\mathit{eI}}_{{\mathit{m}}_1,{\mathit{m}}_2,50,50}\left(\mathit{F}\right)$
2	2	$3.958\times {10}^{-1}$	$6.041\times {10}^{-2}$	$5.833\times {10}^{-2}$	$7.916\times {10}^{-3}$
5	10	$4.083\times {10}^{-2}$	$7.416\times {10}^{-3}$	$4.833\times {10}^{-3}$	$8.166\times {10}^{-4}$
10	5	$3.833\times {10}^{-2}$	$4.666\times {10}^{-3}$	$6.833\times {10}^{-3}$	$7.666\times {10}^{-4}$
25	25	$2.533\times {10}^{-3}$	$3.866\times {10}^{-4}$	$3.733\times {10}^{-4}$	$5.066\times {10}^{-5}$
50	50	$6.333\times {10}^{-4}$	$9.666\times {10}^{-5}$	$9.333\times {10}^{-5}$	$1.266\times {10}^{-5}$

,

Table 4. The approximation of the double integral ${\int }_{1.4}^2{\int }_1^{1.5}ln(x+2y)dxdy$.

${\mathit{m}}_1$	${\mathit{m}}_2$	${\mathit{eI}}_{{\mathit{m}}_1,{\mathit{m}}_2,1,1}\left(\mathit{F}\right)$	${\mathit{eI}}_{{\mathit{m}}_1,{\mathit{m}}_2,5,10}\left(\mathit{F}\right)$	${\mathit{eI}}_{{\mathit{m}}_1,{\mathit{m}}_2,10,5}\left(\mathit{F}\right)$	${\mathit{eI}}_{{\mathit{m}}_1,{\mathit{m}}_2,50,50}\left(\mathit{F}\right)$
2	2	$4.917\times {10}^{-4}$	$6.210\times {10}^{-5}$	$8.519\times {10}^{-5}$	$9.817\times {10}^{-6}$
5	10	$3.526\times {10}^{-5}$	$5.608\times {10}^{-6}$	$4.972\times {10}^{-6}$	$7.054\times {10}^{-7}$
10	5	$6.300\times {10}^{-5}$	$6.820\times {10}^{-6}$	$1.208\times {10}^{-5}$	$1.258\times {10}^{-6}$
25	25	$3.145\times {10}^{-6}$	$3.975\times {10}^{-7}$	$5.457\times {10}^{-7}$	$5.820\times {10}^{-8}$
50	50	$7.872\times {10}^{-7}$	$9.950\times {10}^{-8}$	$1.371\times {10}^{-7}$	$1.870\times {10}^{-8}$

,

Table 5. The approximation of the double integral ${\int }_{-1}^1{\int }_{-1}^1{e}^{-\left(x^2+y^2\right)}dxdy$.

${\mathit{m}}_1$	${\mathit{m}}_2$	${\mathit{eI}}_{{\mathit{m}}_1,{\mathit{m}}_2,1,1}\left(\mathit{F}\right)$	${\mathit{eI}}_{{\mathit{m}}_1,{\mathit{m}}_2,5,10}\left(\mathit{F}\right)$	${\mathit{eI}}_{{\mathit{m}}_1,{\mathit{m}}_2,10,5}\left(\mathit{F}\right)$	${\mathit{eI}}_{{\mathit{m}}_1,{\mathit{m}}_2,50,50}\left(\mathit{F}\right)$
2	2	$3.598\times {10}^{-1}$	$5.585\times {10}^{-2}$	$5.585\times {10}^{-2}$	$7.504\times {10}^{-3}$
5	10	$3.669\times {10}^{-2}$	$6.618\times {10}^{-3}$	$4.409\times {10}^{-3}$	$7.358\times {10}^{-4}$
10	5	$3.669\times {10}^{-2}$	$4.409\times {10}^{-3}$	$6.618\times {10}^{-3}$	$7.358\times {10}^{-4}$
25	25	$2.344\times {10}^{-3}$	$3.517\times {10}^{-4}$	$3.517\times {10}^{-4}$	$4.676\times {10}^{-5}$
50	50	$5.861\times {10}^{-4}$	$8.790\times {10}^{-5}$	$8.790\times {10}^{-5}$	$1.186\times {10}^{-5}$

and

Table 6. The approximation of the double integral ${\int }_0^4{\int }_0^3{e}^{-\left(x+y\right)}·sin(2x+2y)dxdy$.

${\mathit{m}}_1$	${\mathit{m}}_2$	${\mathit{eI}}_{{\mathit{m}}_1,{\mathit{m}}_2,1,1}\left(\mathit{F}\right)$	${\mathit{eI}}_{{\mathit{m}}_1,{\mathit{m}}_2,5,10}\left(\mathit{F}\right)$	${\mathit{eI}}_{{\mathit{m}}_1,{\mathit{m}}_2,10,5}\left(\mathit{F}\right)$	${\mathit{eI}}_{{\mathit{m}}_1,{\mathit{m}}_2,50,50}\left(\mathit{F}\right)$
2	2	$2.078\times {10}^{-1}$	$1.736\times {10}^{-2}$	$1.029\times {10}^{-2}$	$3.031\times {10}^{-3}$
5	10	$5.753\times {10}^{-4}$	$3.319\times {10}^{-4}$	$1.679\times {10}^{-4}$	$4.267\times {10}^{-5}$
10	5	$8.876\times {10}^{-4}$	$3.830\times {10}^{-5}$	$7.770\times {10}^{-5}$	$1.336\times {10}^{-5}$
25	25	$1.623\times {10}^{-5}$	$1.317\times {10}^{-6}$	$1.192\times {10}^{-6}$	$1.344\times {10}^{-7}$
50	50	$3.079\times {10}^{-6}$	$4.204\times {10}^{-7}$	$3.534\times {10}^{-7}$	$4.990\times {10}^{-8}$

), we can remark the following aspects:

• 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.
• 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.

6. Conclusions

Concerning the classical bivariate Bernstein formula, we may conclude the following:

• It is a valuable tool to approximate any volume resulted by the integration of a given function on the interval $[a,b]\times [c,d]$ (on the symmetrical interval $[-a,a]\times [-a,a]$) and which cannot be solved analytically;
• The simplest form (by taking only the minimum number of nodes) is the bivariate trapezoidal formula;
• Its use has certain advantages in comparison with the well-known Newton–Cotés cubature formulas;
• 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.