1 Introduction

In many fields of mathematics several types of differences are used: finite difference, \(\Delta\)-divided difference, \(\nabla\)-difference etc. From the point of view of application there are various implementations of these differences in fields such as numerical analysis, statistics, vector calculus and physics [3,5,6,19]. We chose the \(\nabla\)-operator due to its wide application. In Section 2, we obtain discrete identities for function \(f(y_i, z_j)\) and sequence \(a_{ij}\) involving \(\nabla\)-divided differences. In Section 3, we state general integral identities by different methods for differentiable functions of higher order for two variables of the Popoviciu type. These identities are a generalization of various established results; their application is also given. The present article is related to our previous article [1].

Let us recall some useful definitions and significant results from [1,8,10,11,16]. Throughout this article, I is an interval in \(\mathbb{R}\). Also, we use the following notations for \(\mathbb{R} = (-\infty, \infty)\), \(\mathbb{R}_* = [0, \infty)\) and \(\mathbb{R}_+ = (0, \infty)\).

Definition 1 The m-th order divided difference of a function \(f: I \to \mathbb{R}\), at distinct elements \(y_i, y_{i+1}, ..., y_{i+m} \in I = [a, b] \subset \mathbb{R}\), where \(i \in \mathbb{N}\), is stated as:

\[[y_j; f] = f(y_j), \text{ where } j \in \{i, i+1, ..., i+m\}\]

\[[y_i, y_{i+1}, \dots, y_{i+m}; f] = \frac{[y_{i+1}, \dots, y_{i+m}; f] - [y_i, \dots, y_{i+m-1}; f]}{y_{i+m} - y_i}.\]

This can be written as:

\[[y_i, ..., y_{i+m}; f] = \sum_{k=0}^m \frac{f(y_{i+k})}{\prod_{\substack{i=1, j \neq i+k \ (y_{i+k}-y_j)}}^{i+m}}\]

Remark 2 Some important remarks:

• 1. Let us denote \([y_i, ..., y_{i+m}; f]\) by \(\Delta_{(m)} f(y_i)\). The value \([y_i, ..., y_{i+m}; f]\) is independent of element order \(y_i, y_{i+1}, ..., y_{i+m}\).
• 2. We can extend this definition by including a case in which some elements or all elements coincide by supposing that \(y_i \le \dots \le y_{i+m}\) (see [16]) and letting

\[[y_i, ..., y_{i+m}; f] = \frac{f^{(m)}(y_i)}{m!},\]

provided that \(f^{(m)}(y_i)\) exists.

Definition 3 A function \(f: I \to \mathbb{R}\), is called m-convex or m-th order convex if the inequality \(\Delta_{(m)} f(y_i) \ge 0\) holds \(\forall (m+1)\) different points \(y_i, ..., y_{i+m} \in I\).

Further, if an m-th order derivative of the function exists, then the function is convex of order m if and only if \(f^{(m)} \ge 0\).

Definition 4 Any function \(f: I \to \mathbb{R}\) is known as \(\nabla\)-convex of order m or m-convex if \(\forall\) (m+1) different points \(y_i, y_{i+1}, ..., y_{i+m} \in I\) and we have \(\nabla_{(m)} f(y_i) = (-1)^m \Delta_{(m)} f(y_i) \geq 0\).

Further, if an m-th order derivative of the function exists, then the function is \(\nabla\)-convex of order m if and only if \((-1)^m f^{(m)} \ge 0\).

Definition 5 Let \(E = \{y_1, y_2, ..., y_M\} \subset \mathbb{R}\). Any function \(f: E \to \mathbb{R}\) is called a discrete m-convex function if inequality \([y_i, ..., y_{i+m}; f] \ge 0\) holds \(\forall (m+1)\) different points \(y_i, ..., y_{i+m} \in E\).

We extend the aforementioned definitions up to order (m, n). For this purpose, let us denote \(I \times J = [a, b] \times [c, d] \subset \mathbb{R}^2\).

Definition 6 Let \(f: I \times J \to \mathbb{R}\) be a function, then the divided difference of order (m, n) of f at different elements \(y_i, ..., y_{i+m} \in I, z_j, ..., z_{j+n} \in J\) for some \(i, j \in N\) is stated as \(\Delta_{(m,n)} f(y_i, z_j) = [y_i, ..., y_{i+m}; [z_j, ..., z_{j+n}; f]].\)

Definition 7 A function \(f: I \times J \to \mathbb{R}\) is known as (m, n)-convex if \(\forall\) different elements \(y_i, ..., y_{i+m} \in I\) and \(z_i, ..., z_{i+n} \in J\), and we have \(\Delta_{(m,n)} f(y_i, z_i) \ge 0\).

Further that function f is (m,n)-convex iff \(f_{(m,n)} \ge 0\), given that a partial derivative \(\frac{\partial^{m+n} f}{\partial v^m \partial z^n}\) denoted by \(f_{(m,n)}\) exists.

Definition 8 Let \(E = \{y_1, y_2, ..., y_M\}, F = \{z_1, z_2, ..., z_N\} \subset \mathbb{R}\). Any function \(f: E \times F \to R\) is known as a discrete (m, n)-convex function if inequality \([y_i, ..., y_{i+m}; [z_j, ..., z_{j+n}; f]] \ge 0\) holds \(\forall\) (m+1) different points \(y_i, ..., y_{i+m} \in E\) and (n+1) different points \(z_j, ..., z_{j+n} \in F\).

Definition 9 The finite difference of the function \(f: I \times J \to \mathbb{R}\) of order (m,n), where \(h, k \in R\) and \(y \in I, z \in J\), is stated as:

\[\begin{split} & \Delta_{h,k}^{m,n} f(y,z) = \Delta_{h}^{m} (\Delta_{k}^{n} f(y,z)) = \Delta_{k}^{n} (\Delta_{h}^{m} f(y,z)) \\ &= \sum_{i=0}^{m} \sum_{j=0}^{n} (-1)^{m-i+n-j} \binom{m}{i} \binom{n}{j} f(y+ih,z+jk), \end{split}\] where \(y+ih,z+jk \in I,J\) respectively and \(i \in \{0,1,2,...,m\}\) and \(j \in \{0,1,2,...,n\}\). Moreover, a function \(f:I \times J \to \mathbb{R}\) is called (m,n)-convex if \(\Delta_{h,k}^{m,n}f(y,z) \geq 0 \ \forall \ y \in I, \ z \in J\).

Definition 10 The finite difference and divided difference of (m,n) order of a sequence \((a_{ij})\) are stated as \(\Delta^{m,n}a_{ij} = \Delta^{m,n}_{1,1}f(y_i,z_j)\) and \(\Delta_{(m,n)}a_{ij} = \Delta_{(m,n)}f(y_i,z_j)\), respectively, where \(i \in \{1,2,3,...,m\}\), \(j \in \{1,2,3,...,n\}\). If \(y_i = i\), \(z_j = j\), then \(f:\{1,...,m\} \times \{1,...,n\} \to \mathbb{R}\) is a function \(f(i,j) = a_{ij}\). Moreover, a sequence \((a_{ij})\) is called (m,n)-convex if \(\Delta^{m,n}a_{ij} \geq 0\) for \(m,n \geq 0\) and \(i,j \in \{1,2,3,...\}\).

Further, in the present article we use the following notation for some real sequence \((a_m)\), \(m \in \mathbb{N}\) and \(n \in \{2,3,...\}\).

\[\nabla^{(1)}a_m=\nabla a_m=a_m-a_{m+1}, \ \ \nabla^{(n)}a_m=\nabla(\nabla^{(n-1)}a_m).\]

Also for m distinct real numbers \(y_i\), \(i \in \{1, ..., m\}\) and \(n \ge 0\).

\[(y_k - y_i)^{\{n+1\}} = (y_k - y_i)(y_{k-1} - y_i) \dots (y_{k-n} - y_i),\]
\[(y_k - y_i)^{\{0\}} = 1.\]

Definition 11 We say the function \(f: I \times J \to \mathbb{R}\) is \((m,n) - \nabla\) -convex if inequality \(\nabla_{(m,n)} f(y_i,z_j) = (-1)^{m+n} \Delta_{(m,n)} f(y_i,z_j) \ge 0\) holds \(\forall\) different points \(y_i, \dots, y_{i+m} \in I, z_j, \dots, z_{j+n} \in J\).

Let us first give a brief explanation of the format of our article. Section 1 gave an introduction and some preliminaries. In Section 2, we present identities for the sums \(\sum_{i=1}^{M} \sum_{j=1}^{N} p_{ij} f(y_i, z_j)\) for two dimensions using the \(\nabla\)-operator. Section 3 is devoted to integral case \(\iint P(y, z) f(y, z) dy dz\) for higher-order

differentiable functions of two variables using different techniques. In Section 4 we present an application of the mean value theorem and the last section contains the conclusion.

2 Discrete Case for Function and Sequence of Two Dimensions

Under the given heading, we give the identity of \(\sum_{i=1}^{M} \sum_{j=1}^{N} p_{ij} f(y_i, z_j)\) (see [1]), which involves \(\nabla\) divided differences. We also give the identity of a sequence that involves \(\nabla\) divided differences. Moreover, we split this sequence into two as a special case by using \(a_{ij} = a_i b_j\).

Theorem 12 Let \(p_{ij} \in R\) and \(f: I_1 \times J_1 \to \mathbb{R}\) be a discrete function, where \(i \in \{1,2,3,...,M-1,M\}\) and \(j \in \{1,2,3,...,N-1,N\}\). Then the following identity holds:

\[\sum_{i=1}^{M} \sum_{j=1}^{N} p_{ij} f(y_i, z_j) = \sum_{k=0}^{n-1} \sum_{t=0}^{m-1} \sum_{s=1}^{M-t} \sum_{r=1}^{N-k} p_{sr}(z_N - z_r)^{\{k\}} (y_M - y_s)^{\{t\}} \nabla_{(t,k)} f(y_{M-t}, z_{N-k}) + \sum_{k=0}^{n-1} \sum_{t=1}^{M-m} \sum_{s=1}^{t} \sum_{r=1}^{N-k} p_{sr}(z_N - z_r)^{\{k\}} (y_{t+m-1} - y_s)^{\{m-1\}} \nabla_{(m,k)} f(y_t, z_{N-k}) (y_{t+m} - y_t) + \sum_{k=1}^{N-n} \sum_{t=0}^{m-1} \sum_{s=1}^{M-t} \sum_{r=1}^{k} p_{sr}(z_{k+n-1} - z_r)^{\{n-1\}} (y_M - y_s)^{\{t\}} \nabla_{(t,n)} f(y_{M-t}, z_k) (z_{k+n} - z_k) + \sum_{k=1}^{N-n} \sum_{t=1}^{M-m} \sum_{s=1}^{t} \sum_{r=1}^{k} p_{sr}(z_{k+n-1} - z_r)^{\{n-1\}} (y_{t+m-1} - y_s)^{\{m-1\}} \nabla_{(m,n)} f(y_t, z_k) \times (y_{k+m} - y_k) (z_{k+m} - z_k)\] where \((y_i, z_i) \in I_1 \times J_1\) are distinct points.

Corollary 13 Let ∈ ℝ and be a sequence, where ∈ {1,2,3, … , } and ∈ {1,2,3, …, }. Then the following identity holds:

\[\sum_{i=1}^{M} \sum_{j=1}^{N} p_{ij} a_{ij} \\ = \sum_{k=0}^{n-1} \sum_{t=0}^{m-1} \sum_{s=1}^{M-t} \sum_{r=1}^{N-k} p_{sr} \frac{(M-s)^{\{t\}}}{t!} \frac{(N-r)^{\{k\}}}{k!} \nabla_{(t,k)} a_{(M-t,N-k)} \\ + \sum_{k=0}^{n-1} \sum_{t=1}^{M-m} \sum_{s=1}^{t} \sum_{r=1}^{N-k} p_{sr} \frac{(t-s+m-1)^{\{m-1\}}}{(m-1)!} \frac{(N-r)^{\{k\}}}{k!} \\ + \sum_{k=1}^{N-n} \sum_{t=0}^{m-1} \sum_{s=1}^{M-t} \sum_{r=1}^{k} p_{sr} \frac{(M-s)^{\{t\}}}{t!} \\ + \sum_{k=1}^{N-n} \sum_{t=1}^{M-m} \sum_{s=1}^{t} \sum_{r=1}^{k} p_{sr} \frac{(M-r)^{\{t-1\}}}{(n-1)!} \nabla_{(t,n)} a_{(M-t,k)} \\ + \sum_{k=1}^{N-n} \sum_{t=1}^{M-m} \sum_{s=1}^{t} \sum_{r=1}^{k} p_{sr} \frac{(k-r+n-1)^{\{n-1\}}}{(n-1)!} \\ \frac{(t-s+m-1)^{\{m-1\}}}{(m-1)!} \nabla_{(m,n)} a_{(t,k)}.\]

Remark 14 We can easily obtain the proof of the above corollary by the same method as was used in Theorem 12 of this article [1].

Remark 15 If we simply put = in Corollary 13, then we obtain the same result for two and sequences as follows:

Corollary 16 Let ∈ ℝ, : ↦ and : ↦ be two sequences, where ∈ {1,2,3, … , } and ∈ {1,2,3, ⋯ , }, then:

\[\sum_{i=1}^{M} \sum_{j=1}^{N} p_{ij} a_i b_j\] \[= \sum_{k=0}^{n-1} \sum_{t=0}^{m-1} \sum_{s=1}^{M-t} \sum_{r=1}^{N-k} p_{sr} \frac{(N-r)^{\{k\}}}{k!} \nabla_{(k)} b_{(N-k)} \frac{(M-s)^{\{t\}}}{t!}\]

\[\begin{split} & \nabla_{(t)} a_{(M-t)} \\ & + \sum_{k=0}^{n-1} \sum_{t=1}^{M-m} \sum_{s=1}^{t} \sum_{r=1}^{N-k} p_{sr} \frac{(N-r)^{\{k\}}}{k!} \nabla_{(k)} b_{(N-k)} \\ & \qquad \qquad \frac{(t-s+m-1)^{\{m-1\}}}{(m-1)!} \nabla_{(m)} a_{(t)} \\ & \qquad \qquad + \sum_{k=1}^{N-n} \sum_{t=0}^{m-1} \sum_{s=1}^{M-t} \sum_{r=1}^{k} p_{sr} \frac{(k-r+n-1)^{\{n-1\}}}{(n-1)!} \\ & \qquad \qquad \nabla_{(n)} b_{(k)} \frac{(M-s)^{\{t\}}}{t!} \nabla_{(t)} a_{(M-t)} \\ & \qquad \qquad + \sum_{k=1}^{N-n} \sum_{t=1}^{M-m} \sum_{s=1}^{t} \sum_{r=1}^{k} p_{sr} \frac{(k-r+n-1)^{\{n-1\}}}{(n-1)!} \\ & \qquad \qquad \nabla_{(n)} b_{(k)} \frac{(t-s+m-1)^{\{m-1\}}}{(m-1)!} \times \nabla_{(m)} a_{(t)}. \end{split}\]

3 Integral Case for a Function of Two Variables

Suppose y and z are real continuous variables defined on \(I \times J = [a, b] \times [c, d]\), and \(m, n, M, N \in \mathbb{N} \cup \{0\}\). Throughout this section we will use:

\[f_{(0,0)} = f, f_{(1,0)} = \frac{\partial f}{\partial y'}, f_{(0,1)} = \frac{\partial f}{\partial z'},\] \[f_{(1,1)} = \frac{\partial^2 f}{\partial y \partial z} = \frac{\partial^2 f}{\partial z \partial y'}, f_{(i,j)} = \frac{\partial^{i+j} f}{\partial y^i \partial z^j} = \frac{\partial^{i+j} f}{\partial z^j \partial y^{i'}}\]

Let \(f: I \times J \to R\) and \(P_{i,j}\) both be integrable functions. We introduce the following notations:

\[P_{(1,1)}(y,z) = \int_{y}^{b} \int_{z}^{d} P(s,t) dt ds,\]
\[P_{(m+1,n+1)}(y,z) = \int_{y}^{b} \int_{z}^{d} P_{(m,n)}(s,t) dt ds\] and

\[P_{(m+1,n+1)}(y,z) = \int_{y}^{b} \int_{z}^{d} P(s,t) \frac{(y-s)^{m}}{m!} \frac{(z-t)^{n}}{n!} dt ds.\]

To prove our next theorem we use the two-dimensional induction method as follows:

Definition 17 Let (, ) denote a statement involving two variables and . Suppose

• 1. (0,0) is true.
• 2. If (, 0) is true for some non-negative integer , then ( + 1,0) is also true.
• 3. If (, ) holds for some non-negative integers and , then (, + 1) is also true.

Therefore, (, ) are true ∀ non-negative integers and . This process is called two-dimensional induction.

Remark 18 Here it should be noted that in the above definition of twodimensional induction we used ℕ ∪ {0} instead of an ordinary set of natural number ℕ that starts from 1, see for reference [2, p.~12].

The following lemma is a special case of Theorem 20.

Lemma 19 Let have the continuous partial derivatives (0,1) , (1,0) and (1,1) and , : × ⟶ ℝ be both integrable functions, then we have:

\[\begin{split} &\int_{c}^{d} \int_{a}^{b} P(y,z) f(y,z) dy dz \\ &= \int_{d}^{c} \int_{b}^{a} P(s,t) f(s,t) ds dt \\ &= f_{(0,0)}(b,d) P_{(1,1)}(b,d) + \int_{b}^{a} f_{(1,0)}(s,d) P_{(1,1)}(s,d) ds \\ &+ \int_{d}^{c} f_{(0,1)}(b,t) P_{(1,1)}(b,t) dt + \int_{d}^{c} \int_{b}^{a} f_{(1,1)}(s,t) P_{(1,1)}(s,t) ds dt. \end{split}\]

Now we state the main integral identity of this section using higher-order derivatives. We prove the following result in two different ways, first by twodimensional induction and then by Taylor expansion, recalling Theorem 4.2 from [10].

Theorem 20 Let have continuous partial derivatives (,) and , : × → ℝ both be integrable functions, where ∈ {0,1,2, ⋯ , , + 1}, ∈ {0,1,2, … , , + 1}. Then:

\[\int_{a}^{b} \int_{c}^{d} P(y,z) f(y,z) dz dy\] \[= \sum_{i=0}^{M} \sum_{j=0}^{N} \int_{a}^{b} \int_{c}^{d} P(y,z) \frac{(b-y)^{i}}{i!} \frac{(d-z)^{j}}{j!} (-1)^{i+j}\] (3.1)

\[f_{(i,j)}(b,d)dz dy + \sum_{j=0}^{N} \int_{a}^{b} \int_{c}^{s} \int_{c}^{d} P(y,z) \frac{(s-y)^{M}}{M!} \frac{(d-z)^{j}}{j!} (-1)^{M+j+1} f_{(M+1,j)}(s,d)dz dy ds + \sum_{i=0}^{M} \int_{c}^{d} \int_{a}^{b} \int_{c}^{t} P(y,z) \frac{(b-y)^{i}}{i!} \frac{(t-z)^{N}}{N!} (-1)^{i+N+1} f_{(i,N+1)}(b,t)dz dy dt + \int_{a}^{b} \int_{c}^{d} \int_{a}^{s} \int_{c}^{t} P(y,z) \frac{(s-y)^{M}}{M!} \frac{(t-z)^{N}}{N!} (-1)^{M+N} f_{(M+1,N+1)}(s,t)dz dy dt ds.\]

Proof. (Method I) First we claim that:

\[\int_{d}^{c} \int_{b}^{a} P(y,z)f(y,z)dydz\] \[= \sum_{i=0}^{M} \sum_{j=0}^{N} f_{(i,j)}(b,d)P_{(i+1,j+1)}(b,d)\] \[+ \sum_{j=0}^{N} \int_{b}^{a} f_{(M+1,j)}(s,d)P_{(M+1,j+1)}(s,d)ds\] \[+ \sum_{i=0}^{M} \int_{d}^{c} f_{(i,N+1)}(b,t)P_{(i+1,N+1)}(b,t)dt\] \[+ \int_{d}^{c} \int_{b}^{a} P_{(M+1,N+1)}(s,t)f_{(M+1,N+1)}(s,t)dsdt\] (3.2)

At this stage we show this equality by applying Definition 17 of twodimensional induction, where first we consider a fundamental case, i.e. = = 0.

\[\begin{split} &\int_{c}^{d} \int_{a}^{b} P(y,z) f(y,z) dy dz = \int_{d}^{c} \int_{b}^{a} P(s,t) f(s,t) ds dt \\ &= f_{(0,0)}(b,d) P_{(1,1)}(b,d) + \int_{b}^{a} f_{(1,0)}(s,d) P_{(1,1)}(s,d) ds \\ &+ \int_{d}^{c} f_{(0,1)}(b,t) P_{(1,1)}(b,t) dt + \int_{d}^{c} \int_{b}^{a} f_{(1,1)}(s,t) P_{(1,1)}(s,t) ds dt, \end{split}\]

We will prove that the above is a consequence of Lemma 19.

Now we have to suppose that our hypothesis is true for = (any arbitrary) and = 0, i.e.

\[\int_{d}^{c} \int_{b}^{a} P(y,z)f(y,z)dydz\] \[= \sum_{i=0}^{m} f_{(i,0)}(b,d)P_{(i+1,1)}(b,d)\] \[+ \int_{b}^{a} f_{(m+1,0)}(s,d)P_{(m+1,1)}(s,d)ds\] \[+ \sum_{i=0}^{m} \int_{d}^{c} f_{(i,1)}(b,t)P_{(i+1,1)}(b,t)dt\] \[+ \int_{d}^{c} \int_{b}^{a} f_{(m+1,1)}(s,t)P_{(m+1,1)}(s,t)dsdt,\] (3.3)

We will show that this is also true for = + 1 and = 0, i.e., the below identity holds:

\[\int_{d}^{c} \int_{b}^{a} P(y,z)f(y,z)dydz\] \[= \sum_{i=0}^{m+1} f_{(i,0)}(b,d)P_{(i+1,1)}(b,d)\] \[+ \int_{b}^{a} f_{(m+2,0)}(s,d)P_{(m+2,1)}(s,d)ds\] \[+ \sum_{i=0}^{m+1} \int_{d}^{c} f_{(i,1)}(b,t)P_{(i+1,1)}(b,t)dt\] \[+ \int_{d}^{c} \int_{b}^{a} f_{(m+2,1)}(s,t)P_{(m+2,1)}(s,t)dsdt.\] (3.4)

To prove (3.4), we consider the second term of (3.3):

\[\int_{b}^{a} f_{(m+1,0)}(s,d) P_{(m+1,1)}(s,d) ds = \int_{b}^{a} \left( f_{(m+1,0)}(b,d) + \int_{b}^{s} f_{(m+2,0)}(\theta,d) d\theta \right) P_{(m+1,1)}(s,d) ds\]

\[= f_{(m+1,0)}(b,d) \int_b^a P_{(m+1,1)}(s,d) ds + \int_b^a \int_b^s f_{(m+2,0)}(\theta,d) P_{(m+1,1)}(s,d) d\theta ds,\] by using the Fubini theorem and interchanging \(\theta \leftrightarrow s\):

\[= f_{(m+1,0)}(b,d)P_{(m+2,1)}(b,d) + \int_{b}^{a} \int_{\theta}^{a} f_{(m+2,0)}(\theta,d)P_{(m+1,1)}(s,d)dsd\theta = f_{(m+1,0)}(b,d)P_{(m+2,1)}(b,d) + \int_{b}^{a} f_{(m+2,0)}(s,d) \int_{s}^{a} P_{(m+1,1)}(\theta,d)d\theta ds = f_{(m+1,0)}(b,d)P_{(m+2,1)}(b,d) + \int_{b}^{a} f_{(m+2,0)}(s,d)P_{(m+2,1)}(s,d)ds.\](3.5)

Now consider the fourth term of (3.3):

\[\int_{d}^{c} \int_{b}^{a} f_{(m+1,1)}(s,t) P_{(m+1,1)}(s,t) ds dt = \int_{d}^{c} \int_{b}^{a} \left( f_{(m+1,1)}(b,t) + \int_{b}^{s} f_{(m+2,1)}(\theta,t) d\theta \right) P_{(m+1,1)}(s,t) ds dt = \int_{d}^{c} f_{(m+1,1)}(b,t) \int_{b}^{a} P_{(m+1,1)}(s,t) ds dt + \int_{d}^{c} \int_{b}^{a} \int_{b}^{s} f_{(m+2,1)}(\theta,t) P_{(m+1,1)}(s,t) d\theta ds dt = \int_{d}^{c} f_{(m+1,1)}(b,t) P_{(m+2,1)}(b,t) dt + \int_{d}^{c} \int_{b}^{a} \int_{\theta}^{a} f_{(m+2,1)}(\theta,t) P_{(m+1,1)}(s,t) ds d\theta dt = \int_{d}^{c} f_{(m+1,1)}(b,t) P_{(m+2,1)}(b,t) dt + \int_{d}^{c} \int_{b}^{a} f_{(m+2,1)}(s,t) \int_{s}^{a} P_{(m+1,1)}(\theta,t) d\theta ds dt = \int_{d}^{c} f_{(m+1,1)}(b,t) P_{(m+2,1)}(b,t) dt + \int_{d}^{c} \int_{b}^{a} f_{(m+2,1)}(s,t) P_{(m+2,1)}(s,t) ds dt.\] (3.6)

Substituting the values of Eq. (3.5) and Eq. (3.6) in Eq. (3.3), we obtain the equal Eq. (3.4), which means our first hypothesis is true.

Again, set the hypothesis for M = m (any arbitrary) and N = n, i.e.

\[\int_{d}^{c} \int_{b}^{a} P(y,z) f(y,z) dy dz = \sum_{i=0}^{m} \sum_{j=0}^{n} f_{(i,j)}(b,d) P_{(i+1,j+1)}(b,d)\] (3.7)

\[+\sum_{j=0}^{n} \int_{b}^{a} f_{(m+1,j)}(s,d) P_{(m+1,j+1)}(s,d) ds +\sum_{i=0}^{m} \int_{d}^{c} f_{(i,n+1)}(b,t) P_{(i+1,n+1)}(b,t) dt +\int_{d}^{c} \int_{b}^{a} P_{(m+1,n+1)}(s,t) f_{(m+1,n+1)}(s,t) ds dt.\]

Now we show that it is also valid for M = m and N = n + 1:

\[\int_{d}^{c} \int_{b}^{a} P(y,z) f(y,z) dy dz \qquad (3.8)\] \[= \sum_{i=0}^{m} \sum_{j=0}^{n+1} f_{(i,j)}(b,d) P_{(i+1,j+1)}(b,d) + \sum_{j=0}^{n+1} \int_{b}^{a} f_{(m+1,j)}(s,d) P_{(m+1,j+1)}(s,d) ds + \sum_{i=0}^{m} \int_{d}^{c} f_{(i,n+2)}(b,t) P_{(i+1,n+2)}(b,t) dt + \int_{d}^{c} \int_{b}^{a} f_{(m+1,n+2)}(s,t) P_{(m+1,n+2)}(s,t) ds dt.\]

To prove (3.8) we consider third term of (3.7):

\[\Sigma_{i=0}^{m} \int_{d}^{c} f_{(i,n+1)}(b,t) P_{(i+1,n+1)}(b,t) dt = \Sigma_{i=0}^{m} \int_{d}^{c} \left( f_{(i,n+1)}(b,d) + \int_{d}^{t} f_{(i,n+2)}(b,\phi) d\phi \right) P_{(i+1,n+1)}(b,t) dt = \Sigma_{i=0}^{m} \left( f_{(i,n+1)}(b,d) \int_{d}^{c} P_{(i+1,n+1)}(b,t) dt + \int_{d}^{c} \int_{d}^{t} f_{(i,n+2)}(b,\phi) P_{(i+1,n+1)}(b,t) d\phi dt \right) = \Sigma_{i=0}^{m} \left( f_{(i,n+1)}(b,d) P_{(i+1,n+2)}(b,d) + \int_{d}^{c} \int_{\phi}^{c} f_{(i,n+2)}(b,\phi) P_{(i+1,n+1)}(b,t) dt d\phi \right) = \Sigma_{i=0}^{m} \left( f_{(i,n+1)}(b,d) P_{(i+1,n+2)}(b,d) + \int_{d}^{c} f_{(i,n+2)}(b,t) \int_{t}^{c} P_{(i+1,n+1)}(b,\phi) d\phi dt \right) = \Sigma_{i=0}^{m} \left( f_{(i,n+1)}(b,d) P_{(i+1,n+2)}(b,d) + \int_{d}^{c} f_{(i,n+2)}(b,t) P_{(i+1,n+2)}(b,d) + \int_{d}^{c} f_{(i,n+2)}(b,t) P_{(i+1,n+2)}(b,t) dt \right).\] (3.9)

Now we consider the 4th term of (3.7):

\[\int_{d}^{c} \int_{b}^{a} f_{(m+1,n+1)}(s,t) P_{(m+1,n+1)}(s,t) ds dt\] \[= \int_{b}^{a} \int_{d}^{c} f_{(m+1,n+1)}(s,t) P_{(m+1,n+1)}(s,t) dt ds\]

\[= \int_{b}^{a} \int_{d}^{c} \left( f_{(m+1,n+1)}(s,d) + \int_{d}^{t} f_{(m+1,n+2)}(s,\phi) d\phi \right) P_{(m+1,n+1)}(s,t) dt ds\] \[= \int_{b}^{a} f_{(m+1,n+1)}(s,d) \int_{d}^{c} P_{(m+1,n+1)}(s,t) dt ds\] \[+ \int_{b}^{a} \int_{d}^{c} \int_{d}^{t} f_{(m+1,n+2)}(s,\phi) P_{(m+1,n+1)}(s,t) d\phi dt ds\] \[= \int_{b}^{a} f_{(m+1,n+1)}(s,d) P_{(m+1,n+2)}(s,d) ds\] \[+ \int_{b}^{a} \int_{d}^{c} \int_{\phi}^{c} f_{(m+1,n+2)}(s,\phi) P_{(m+1,n+1)}(s,t) dt d\phi ds\] \[= \int_{b}^{a} f_{(m+1,n+1)}(s,d) P_{(m+1,n+2)}(s,d) ds\] \[+ \int_{b}^{a} \int_{d}^{c} f_{(m+1,n+2)}(s,t) \int_{t}^{c} P_{(m+1,n+1)}(s,\phi) d\phi dt ds\] \[= \int_{b}^{a} f_{(m+1,n+1)}(s,d) P_{(m+1,n+2)}(s,d) ds\] \[+ \int_{b}^{a} \int_{d}^{c} f_{(m+1,n+2)}(s,t) P_{(m+1,n+2)}(s,t) dt ds. \tag{3.10}\]

Substituting Eq. (3.9) and Eq. (3.10) in Eq. (3.7), we obtain (3.8), which proves our result by induction.

To complete our proof, we have to apply the notations that were introduced at the beginning of this section. In Eq. (3.2), we consider the terms and find the values one by one, i.e.

\[\begin{split} & \sum_{i=0}^{M} \sum_{j=0}^{N} f_{(i,j)}(b,d) P_{(i+1,j+1)}(b,d) \\ & = \sum_{i=0}^{M} \sum_{j=0}^{N} \int_{b}^{a} \int_{c}^{c} f_{(i,j)}(b,d) P(y,z) \frac{(z-d)^{j}}{j!} \frac{(y-b)^{i}}{i!} dz dy \\ & = \sum_{i=0}^{M} \sum_{j=0}^{N} \int_{a}^{b} \int_{c}^{d} (-1)^{i+j} f_{(i,j)}(b,d) P(y,z) \frac{(d-z)^{j}}{j!} \frac{(b-y)^{i}}{i!} dz dy, \\ & \sum_{j=0}^{N} \int_{b}^{a} f_{(M+1,j)}(s,d) P_{(M+1,j+1)}(s,d) ds \\ & = \sum_{j=0}^{N} \int_{b}^{a} f_{(M+1,j)}(s,d) \left( \int_{s}^{a} \int_{d}^{c} P(y,z) \frac{(y-s)^{M}}{M!} \frac{(z-d)^{j}}{j!} dz dy \right) ds \\ & = \sum_{j=0}^{N} \int_{a}^{b} \int_{a}^{s} \int_{c}^{d} (-1)^{M+1+j} f_{(M+1,j)}(s,d) P(y,z) \frac{(s-y)^{M}}{M!} \frac{(d-z)^{j}}{j!} dz dy ds, \\ & \sum_{i=0}^{M} \int_{d}^{c} f_{(i,N+1)}(b,t) P_{(i+1,N+1)}(b,t) dt \\ & = \sum_{i=0}^{M} \int_{d}^{c} f_{(i,N+1)}(b,t) \left( \int_{b}^{a} \int_{t}^{c} P(y,z) \frac{(y-b)^{i}}{i!} \frac{(z-t)^{N}}{N!} dz dy \right) dt \\ & = \sum_{i=0}^{M} \int_{a}^{b} \int_{c}^{d} \int_{c}^{t} (-1)^{i+N+1} f_{(i,N+1)}(b,t) P(y,z) \frac{(b-y)^{i}}{i!} \frac{(t-z)^{N}}{N!} dz dy dt, \end{split}\] and

\[\int_{b}^{a} \int_{d}^{c} P_{(M+1,N+1)}(s,t) f_{(M+1,N+1)}(s,t) ds dt = \int_{b}^{a} \int_{d}^{c} f_{(M+1,N+1)}(s,t) \left( \int_{s}^{a} \int_{t}^{c} P(y,z) \frac{(z-t)^{N}}{N!} \frac{(y-s)^{M}}{M!} dz dy \right) ds dt = \int_{a}^{b} \int_{c}^{d} \int_{a}^{s} \int_{c}^{t} (-1)^{M+N} f_{(M+1,N+1)}(s,t) P(y,z) \frac{(s-y)^{M}}{M!} \frac{(t-z)^{N}}{N!} dz dy dt ds.\]

Put all these values in (3.2) to get identity (3.1).

Proof. (Method II) Let F(z) = f(y, z), i.e., we consider f(y, z) as a function of z, where y is fixed. Then F may be written as a Taylor expansion:

\[f(y,z) = F(z)\] \[= \sum_{j=0}^{N} F^{(j)}(d) \frac{(z-d)^{j}}{j!} + \int_{d}^{z} F^{(N+1)}(t) \frac{(z-t)^{N}}{N!} dt\] \[= \sum_{j=0}^{N} (-1)^{j} f_{(0,j)}(y,d) \frac{(d-z)^{j}}{j!} + \int_{z}^{d} (-1)^{N+1} f_{(0,N+1)}(y,t) \frac{(t-z)^{N}}{N!} dt,\] where we use \(F^{(j)}(d) = f_{(0,j)}(y,d)\) and \(F^{(N+1)}(t) = f_{(0,N+1)}(y,t)\).

Multiplying the above equation by P(y, z) and integrating it by z over the limit c to d, we get:

\[\int_{c}^{d} P(y,z)f(y,z)dz = \sum_{j=0}^{N} (-1)^{j} f_{(0,j)}(y,d) \int_{c}^{d} P(y,z) \frac{(d-z)^{j}}{j!} dz + \int_{c}^{d} \left( \int_{z}^{d} P(y,z) (-1)^{N+1} f_{(0,N+1)}(y,t) \frac{(t-z)^{N}}{N!} dt \right) dz.\] (3.11)

Now we represent the functions \(y \mapsto f_{(0,j)}(y,d)\) and \(y \mapsto f_{(0,N+1)}(y,t)\) by using Taylor expansion:

\[\begin{split} &f_{(0,j)}(y,d) \\ &= \sum_{i=0}^{M} (-1)^{i} f_{(i,j)}(b,d) \frac{(b-y)^{i}}{i!} + \int_{y}^{b} (-1)^{M+1} f_{(M+1,j)}(s,d) \frac{(s-y)^{M}}{M!} ds, \\ &f_{(0,N+1)}(y,t) \\ &= \sum_{i=0}^{M} (-1)^{i} f_{(i,N+1)}(b,t) \frac{(b-y)^{i}}{i!} + \\ &\int_{y}^{b} (-1)^{M+1} f_{(M+1,N+1)}(s,t) \frac{(s-y)^{M}}{M!} ds. \end{split}\]

Using the above formulae in Eq. (3.11), we have:

\[\begin{split} &\int_{c}^{d} P(y,z)f(y,z)dz \\ &= \sum_{j=0}^{N} (-1)^{j} \left( \sum_{i=0}^{M} (-1)^{i} f_{(i,j)}(b,d) \frac{(b-y)^{i}}{i!} \right. \\ &+ \int_{y}^{b} (-1)^{M+1} f_{(M+1,j)}(s,d) \frac{(s-y)^{M}}{M!} ds \right) \int_{c}^{d} P(y,z) \frac{(d-z)^{j}}{j!} dz \\ &+ \int_{c}^{d} \left( \int_{z}^{d} P(y,z)(-1)^{N+1} \left( \sum_{i=0}^{M} (-1)^{i} f_{(i,N+1)}(b,t) \frac{(b-y)^{i}}{i!} \right. \\ &+ \int_{y}^{b} (-1)^{M+1} f_{(M+1,N+1)}(s,t) \frac{(s-y)^{M}}{M!} ds \right) \frac{(t-z)^{N}}{N!} dt \right) dz \\ &= \sum_{j=0}^{N} \left( \sum_{i=0}^{M} (-1)^{i+j} f_{(i,j)}(b,d) \frac{(b-y)^{i}}{i!} \right) \int_{c}^{d} P(y,z) \frac{(d-z)^{j}}{j!} dz \\ &+ \sum_{j=0}^{N} \left( \int_{y}^{b} (-1)^{M+1+j} f_{(M+1,j)}(s,d) \frac{(s-y)^{M}}{M!} ds \right) \int_{c}^{d} P(y,z) \frac{(d-z)^{j}}{j!} dz \\ &+ \int_{c}^{d} \int_{z}^{d} P(y,z) \left( \sum_{i=0}^{M} (-1)^{i+N+1} f_{(i,N+1)}(b,t) \frac{(b-y)^{i}}{i!} \right) \frac{(t-z)^{N}}{N!} dt dz \\ &+ \int_{c}^{d} \int_{z}^{d} \left( \int_{y}^{b} P(y,z)(-1)^{M+N} f_{(M+1,N+1)}(s,t) \frac{(s-y)^{M}}{M!} ds \right) \frac{(t-z)^{N}}{N!} dt dz. \end{split}\]

Now integrate P(y, z)f(y, z) by y over the limit a to b to obtain:

\[\begin{split} & \int_{a}^{b} \int_{c}^{d} P(y,z) f(y,z) dz \, dy \\ & = \int_{a}^{b} \left[ \sum_{j=0}^{N} \left( \sum_{i=0}^{M} (-1)^{i+j} f_{(i,j)}(b,d) \frac{(b-y)^{i}}{i!} \right) \int_{c}^{d} P(y,z) \frac{(d-z)^{j}}{j!} dz \right] dy \\ & + \int_{a}^{b} \left[ \sum_{j=0}^{N} \left( \int_{y}^{b} (-1)^{M+1+j} f_{(M+1,j)}(s,d) \frac{(s-y)^{M}}{M!} ds \right) \int_{c}^{d} P(y,z) \frac{(d-z)^{j}}{j!} dz \right] dy \\ & + \int_{a}^{b} \left[ \int_{c}^{d} \int_{z}^{d} P(y,z) \left( \sum_{i=0}^{M} (-1)^{i+N+1} f_{(i,N+1)}(b,t) \frac{(b-y)^{i}}{i!} \right) \frac{(t-z)^{N}}{N!} dt \, dz \right] dy \\ & + \int_{a}^{b} \left[ \int_{c}^{d} \int_{z}^{d} \left( \int_{y}^{b} P(y,z) (-1)^{M+N} f_{(M+1,N+1)}(s,t) \frac{(s-y)^{M}}{M!} ds \right) \frac{(t-z)^{N}}{N!} dt \, dz \right] dy. \end{split}\]

We change the order of summation in the first summand and use linearity of integral to obtain:

\[\sum_{i=0}^{M} \sum_{j=0}^{N} \int_{a}^{b} \int_{c}^{d} P(y,z)(-1)^{i+j} f_{(i,j)}(b,d) \frac{(b-y)^{i}}{i!} \frac{(d-z)^{j}}{i!} dz dy.\]

The second summand is rewritten as:

\[\int_{a}^{b} \left[ \sum_{j=0}^{N} \left( \int_{y}^{b} (-1)^{M+1+j} f_{(M+1,j)}(s,d) \frac{(s-y)^{M}}{M!} ds \right) \int_{c}^{d} P(y,z) \frac{(d-z)^{j}}{j!} dz \right] dy = \int_{a}^{b} \left[ \sum_{j=0}^{N} \left( \int_{y}^{b} \int_{c}^{d} P(y,z) \frac{(d-z)^{j}}{j!} (-1)^{M+1+j} f_{(M+1,j)}(s,d) \frac{(s-y)^{M}}{M!} dz ds \right) \right] dy\]

\[= \sum_{j=0}^{N} \int_{a}^{b} \int_{y}^{b} \int_{c}^{d} P(y,z)(-1)^{M+1+j} f_{(M+1,j)}(s,d) \frac{(s-y)^{M}}{M!} \frac{(d-z)^{j}}{j!} dz ds dy\] \[= \sum_{j=0}^{N} \int_{a}^{b} \int_{a}^{s} \int_{c}^{d} P(y,z)(-1)^{M+1+j} f_{(M+1,j)}(s,d) \frac{(s-y)^{M}}{M!} \frac{(d-z)^{j}}{j!} dz dy ds.\]

We use the Fubini theorem for variables s and y in the last step. Let us first change variable y from a to b and then change the variable s from y to b. After changing the integration order, s is changed a to b while y is changed a to s. In the same way the third summand is rewritten as:

\[\begin{split} &\int_{a}^{b} \left[ \int_{c}^{d} \int_{z}^{d} P(y,z) \left( \sum_{i=0}^{M} (-1)^{i+N+1} f_{(i,N+1)}(b,t) \frac{(b-y)^{i}}{i!} \right) \frac{(t-z)^{N}}{N!} dt \ dz \right] dy \\ &= \sum_{i=0}^{M} \int_{a}^{b} \int_{c}^{d} \int_{z}^{d} P(y,z) (-1)^{i+N+1} f_{(i,N+1)}(b,t) \frac{(b-y)^{i}}{i!} \frac{(t-z)^{N}}{N!} dt \ dz \ dy \\ &= \sum_{i=0}^{M} \int_{a}^{b} \int_{c}^{d} \int_{c}^{t} P(y,z) (-1)^{i+N+1} f_{(i,N+1)}(b,t) \frac{(b-y)^{i}}{i!} \frac{(t-z)^{N}}{N!} dz \ dt \ dy \\ &= \sum_{i=0}^{M} \int_{c}^{d} \int_{a}^{b} \int_{c}^{t} P(y,z) (-1)^{i+N+1} f_{(i,N+1)}(b,t) \frac{(b-y)^{i}}{i!} \frac{(t-z)^{N}}{N!} dz \ dy \ dt. \end{split}\]

In the above, use the Fubini theorem twice, first changing t and z, then changing t and y. Therefore, the last summand is rewritten as:

\[\int_{a}^{b} \left[ \int_{c}^{d} \int_{z}^{d} \left( \int_{y}^{b} P(y,z)(-1)^{M+N} f_{(M+1,N+1)}(s,t) \frac{(s-y)^{M}}{M!} ds \right) \frac{(t-z)^{N}}{N!} dt dz \right] dy = \int_{a}^{b} \int_{c}^{d} \int_{z}^{d} \int_{y}^{b} P(y,z)(-1)^{M+N} f_{(M+1,N+1)}(s,t) \frac{(s-y)^{M}}{M!} \frac{(t-z)^{N}}{N!} ds dt dz dy = \int_{a}^{b} \int_{c}^{d} \int_{a}^{s} \int_{c}^{t} P(y,z)(-1)^{M+N} f_{(M+1,N+1)}(s,t) \frac{(s-y)^{M}}{M!} \frac{(t-z)^{N}}{N!} dz dy dt ds.\]

In above, use the Fubini theorem various times, but first change t and z, then z and s, then s and t, then s and y, then t and y. We obtain the required identity by using all the above results:

\[\begin{split} &\int_{a}^{b} \int_{c}^{d} P(y,z) f(y,z) dz \, dy \\ &= \sum_{i=0}^{M} \sum_{j=0}^{N} \int_{a}^{b} \int_{c}^{d} P(y,z) (-1)^{i+j} f_{(i,j)}(b,d) \frac{(b-y)^{i}}{i!} \frac{(d-z)^{j}}{j!} dz \, dy \\ &+ \sum_{j=0}^{N} \int_{a}^{b} \int_{a}^{s} \int_{c}^{d} P(y,z) (-1)^{M+1+j} f_{(M+1,j)}(s,d) \frac{(s-y)^{M}}{M!} \frac{(d-z)^{j}}{j!} dz \, dy \, ds \\ &+ \sum_{i=0}^{M} \int_{c}^{d} \int_{a}^{b} \int_{c}^{t} P(y,z) (-1)^{i+N+1} f_{(i,N+1)}(b,t) \frac{(b-y)^{i}}{i!} \frac{(t-z)^{N}}{N!} dz \, dy \, dt \\ &+ \int_{a}^{b} \int_{c}^{d} \int_{a}^{s} \int_{c}^{t} P(y,z) (-1)^{M+N} f_{(M+1,N+1)}(s,t) \frac{(s-y)^{M}}{M!} \frac{(t-z)^{N}}{N!} dz \, dy \, dt \, ds. \end{split}\]

We present the necessary and sufficient conditions by using the results of the previous theorem, where \(\Lambda(f) \ge 0\) holds \(\forall (M+1,N+1) - \nabla\) -convex function (see [1]).

Theorem 21 Let the suppositions of Theorem 20 be true. Then the following inequality holds:

\[\Lambda(f) = \int_a^b \int_c^d P(y, z) f(y, z) dz \, dy \ge 0, \tag{3.12}\] for all \((M + 1, N + 1) - \nabla\) —convex function f on \(I \times J\), iff

\[\int_{a}^{b} \int_{c}^{d} P(y, z) \frac{(b-y)^{i}}{i!} \frac{(d-z)^{j}}{j!} dz dy = 0, \quad i \in \{0, ..., M\}; \quad j \in \{0, ..., N\}\] (3.13)

\[\int_{a}^{s} \int_{c}^{d} P(y, z) \frac{(s - y)^{M}}{M!} \frac{(d - z)^{j}}{j!} dz dy = 0, \quad j \in \{0, ..., N\}; \quad \forall \ s \in [a, b]\] (3.14)

\[\int_{a}^{b} \int_{c}^{t} P(y, z) \frac{(b-y)^{i}}{i!} \frac{(t-z)^{N}}{N!} dz dy = 0, \quad i \in \{0, \dots, M\}; \quad \forall \ t \in [c, d]\] (3.15)

\[\int_{a}^{s} \int_{c}^{t} P(y, z) \frac{(s - y)^{M}}{M!} \frac{(t - z)^{N}}{N!} dz dy \ge 0, \quad \forall \ s \in [a, b]; \quad \forall \ t \in [c, d].\] (3.16)

Remark 22 We may also obtain a corollary of Theorem 20 for \(I^2\) by changing the variables on the right side: \(y \leftrightarrow s\), \(z \leftrightarrow t\).

Corollary 23 Let \(P, f: I^2 \to R\) be functions such that \(f \in C^{(M+1,N+1)}(I^2)\) and P is integrable, where \(i \in \{0,1,2,\cdots,M,M+1\}, j \in \{0,1,2,\dots,N,N+1\}\). Then

\[\begin{split} & \int_{a}^{b} \int_{a}^{b} P(y,z) f(y,z) dz \, dy \\ & = \sum_{i=0}^{M} \sum_{j=0}^{N} \int_{a}^{b} \int_{a}^{b} P(s,t) \frac{(b-s)^{i}}{i!} \frac{(b-t)^{j}}{j!} (-1)^{i+j} f_{(i,j)}(b,b) dt \, ds \\ & + \sum_{j=0}^{N} \int_{a}^{b} \int_{a}^{y} \int_{a}^{b} P(s,t) \frac{(y-s)^{M}}{M!} \frac{(b-t)^{j}}{j!} (-1)^{M+1+j} f_{(M+1,j)}(y,b) dt \, ds \, dy \\ & + \sum_{i=0}^{M} \int_{a}^{b} \int_{a}^{b} \int_{a}^{z} P(s,t) \frac{(b-s)^{i}}{i!} \frac{(z-t)^{N}}{N!} (-1)^{i+N+1} f_{(i,N+1)}(b,z) dt \, ds dz \\ & + \int_{a}^{b} \int_{a}^{b} \int_{a}^{y} \int_{a}^{z} P(s,t) \frac{(y-s)^{M}}{M!} \frac{(z-t)^{N}}{N!} (-1)^{M+N} f_{(M+1,N+1)}(y,z) dt \, ds \, dz \, dy \end{split}\]

holds.

Remark 24 We have obtained the necessary and sufficient conditions by using the results of the previous theorem, where \(\Lambda(f) \ge 0\) holds \(\forall (M+1,N+1) - \nabla\) -convex function (see [1]).

Now we recall a result from [10].

Theorem 25 Let f have continuous partial derivatives \(f_{(i,j)}\) and \(P, f: I \times J \to \mathbb{R}\), which are both integrable functions, where \(i \in \{0,1,2,\cdots,M,M+1\}, j \in\)\(\{0,1,2,\ldots,N,N+1\}\). Then:

\[\int_{c}^{d} \int_{a}^{b} P(y,z) f(y,z) dy dz \qquad (3.17)\] \[= \sum_{i=0}^{M} \sum_{j=0}^{N} \int_{a}^{b} \int_{c}^{d} f_{(i,j)}(a,c) P(s,t) \frac{(s-a)^{i}}{i!} \frac{(t-c)^{j}}{j!} dt ds + \sum_{j=0}^{N} \int_{a}^{b} \int_{y}^{b} \int_{c}^{d} f_{(M+1,j)}(y,c) P(s,t) \frac{(s-y)^{M}}{M!} \frac{(t-c)^{j}}{j!} dt ds dy + \sum_{i=0}^{M} \int_{c}^{d} \int_{a}^{b} \int_{z}^{d} f_{(i,N+1)}(a,z) P(s,t) \frac{(s-a)^{i}}{i!} \frac{(t-z)^{N}}{N!} dt ds dz + \int_{a}^{b} \int_{c}^{d} \int_{y}^{b} \int_{z}^{d} f_{(M+1,N+1)}(y,z) P(s,t) \frac{(s-y)^{M}}{M!} \frac{(t-z)^{N}}{N!} dt ds dz dy.\]

Remark 26 This above result can also be proved by using two-dimensional mathematical induction as given in the proof of Theorem 20.

4 Application to Mean Value Theorem

It is known that the mean value theorem is a valuable tool for obtaining interesting and important results in classical real analysis. In the field of differential calculus, the most demanding theorem is Lagrange and Cauchy's mean value theorem. Here, we provide a generalized mean value theorem of the Lagrange and Cauchy-type.

Theorem 27 Let \(P: I \times J \to \mathbb{R}\), be an integrable function and \(f \in\)\(C^{(M+1,N+1)}(I \times J)\), be a \((M+1,N+1) - \nabla\) -convex function on the interval \(I \times I\). Let \(\Lambda\) be a linear functional as stated in (3.12) and the conditions (3.13), (3.14), (3.15) and (3.16) be true for function P in Theorem 21, then \(\exists (\eta, \zeta) \in I \times J, \exists\)

\[\Lambda(f) = \Lambda(G_0) f_{(M+1,N+1)}(\eta,\zeta), \tag{4.1}\] where

\[G_0(y,z) = (-1)^{M+N} \frac{y^{M+1}}{(M+1)!} \frac{z^{N+1}}{(N+1)!}\]

Proof: Let

\[U = \max_{(y,z)\in I\times J} (-1)^{M+N} f_{(M+1,N+1)}(y,z),\]
\[L = \min_{(y,z)\in I\times J} (-1)^{M+N} f_{(M+1,N+1)}(y,z).\]

\[L = \min_{(y,z) \in I \times J} (-1)^{r_1 r_2 r_3} f_{(M+1,N+1)}(y,z)\]

Then the function

\[G(y,z) = U(-1)^{M+N} \frac{y^{M+1}}{(M+1)!} \frac{z^{N+1}}{(N+1)!} - f(y,z) = UG_0(y,z) - f(y,z),\]

gives us

\[(-1)^{M+N}G_{(M+1,N+1)}(y,z) = U - (-1)^{M+N}f_{(M+1,N+1)}(y,z) \ge 0,\]

i.e. G is a \(\nabla\)-convex function of order (M+1,N+1) on \(I\times J\). Hence, \(\Lambda(G)\geq 0\) using Theorem 21 and we summarize that \(\Lambda(f)\leq U\Lambda(G_0)\). Similarly, \(L\Lambda(G_0)\leq \Lambda(f)\). Now, we can write the above two inequalities as: \(L\Lambda(G_0)\leq \Lambda(f)\leq U\Lambda(G_0)\), which gives the required result (4.1).

5 Conclusion

This article first gave identities for the sums \(\sum_{i=1}^{M} \sum_{j=1}^{N} p_{ij} f(y_i, z_j)\) involving a function and a sequence of two dimensions using \(\nabla\) divided differences and we also obtained a similar result for two \(a_i\) and \(b_j\) sequences by considering \(a_{ij} = a_i b_j\) in the obtained identity. Secondly, we obtained the integral case \(\iint P(y,z) f(y,z) dy dz\) for differentiable functions of higher order with two variables by using two methods: the first method used was two-dimensional mathematical induction and the second was Taylor expansion. Further, we also gave a similar result for the above integral case in the interval \(I^2\). At the end of article, we presented an application to a mean value theorem.