Professor showed this result in the lecture without giving any proof (after proving the existence of the interpolating polynomial in two variables). I've been trying to prove it myself or find a book where is proved but I failed. This is the theorem:
Let
$$ x_0 < x_1 < \cdots < x_n \in [a, b], \quad y_0 < y_1 < \cdots < y_m \in [c, d],$$
$$ M = \{ (x_i, y_j) : 0 \leq i \leq n, 0 \leq j \leq m \}, \quad f \in \mathcal{C}^{m + n + 2}([a,b] \times [c,d]), $$
$$ p \in \Pi_{n, m} : p(x_i, y_j) = f(x_i, y_j) \quad \forall 0 \leq i \leq n, 0 \leq j \leq m. $$
Then, for all $(x, y) \in (x_0, x_n) \times (y_0, y_m)$ exist $\xi, \xi' \in (x_0, x_n), \eta, \eta' \in (y_0, y_m)$ such that
$$ f(x, y) – p(x, y) = \frac{1}{(n + 1)!} \frac{\partial^{n + 1} f(\xi, y)}{\partial x^{n + 1}} \prod_{i = 0}^n (x – x_i) $$
$$ + \frac{1}{(m + 1)!} \frac{\partial^{m + 1} f(x, \eta)}{\partial y^{m + 1}} \prod_{j = 0}^m (y – y_j) $$
$$ – \frac{1}{(n + 1)! (m + 1)!} \frac{\partial^{n + m + 2} f(\xi', \eta')}{\partial x^{n + 1} \partial y^{m + 1}} \prod_{i = 0}^n (x – x_i) \prod_{j = 0}^m (y – y_j) $$
I appreciate any kind of help, even if it is only from where could I start the proof.
Edit 1 (clarification):
$$ \Pi_{n, m} = \{ p(x, y) = \sum_{i = 0}^n \sum_{j = 0}^m a_{i,j} x^i y^j : a_{i, j} \in \mathbb{R} \quad \forall 0 \leq i \leq n, 0 \leq j \leq m \} $$
Best Answer
Let
$$X_{- 1}(x) = 1,\quad X_{i}(x) = \prod_{k = 0}^{i}\left( x - x_{k} \right)\quad\forall 0 \leq i \leq n $$
$$ Y_{- 1}(y) = 1,\quad Y_{j}(y) = \prod_{l = 0}^{j}\left( y - y_{l} \right)\quad\forall 0 \leq j \leq m $$
Considering $f( \cdot ,y) \in C^{n + 1}$ as a univariate function, we'll have by the Newton form for a single variable for interpolation
$$ f(x,y) = \sum_{i = 0}^{n}X_{i - 1}(x)f\left\lbrack x_{0},\ldots,x_{i};y \right\rbrack + X_{n}(x)f\left\lbrack x_{0},\ldots,x_{n},x;y \right\rbrack $$
where
$$ f\left\lbrack x_{0},\ldots,x_{i};y \right\rbrack = \left( f( \cdot ,y) \right)\left\lbrack x_{0},\ldots,x_{i} \right\rbrack $$
Now, $f\left\lbrack x_{0},\ldots,x_{i};y \right\rbrack$ for $0 \leq k \leq n$ are also univariate functions wrt $y$. Because it implicitely depends on (f \in C^{n + m + 2}) by the recursive formula of divided differences, we'll have $f\left\lbrack x_{0},\ldots,x_{i};y \right\rbrack \in C^{m + 1}$ so we can apply Newton formula
$$f\left\lbrack x_{0},\ldots,x_{i};y \right\rbrack = \sum_{j = 0}^{m}Y_{j - 1}(y)f\left\lbrack x_{0},\ldots,x_{i};y_{0},\ldots,y_{j} \right\rbrack + Y_{m}(y)f\left\lbrack x_{0},\ldots,x_{i};y_{0},\ldots,y_{m},y \right\rbrack$$
where
$$f\left\lbrack x_{0},\ldots,x_{i};y_{0},\ldots,y_{j} \right\rbrack = \left( f\left\lbrack x_{0},\ldots,x_{i}; \cdot \right\rbrack \right)\left\lbrack y_{0},\ldots,y_{m} \right\rbrack$$
So if we substitute
$$f(x,y) = \sum_{i = 0}^{n}X_{i - 1}(x)\left( \sum_{j = 0}^{m}Y_{j - 1}(y)f\left\lbrack x_{0},\ldots,x_{i};y_{0},\ldots,y_{j} \right\rbrack + Y_{m}(y)f\left\lbrack x_{0},\ldots,x_{i};y_{0},\ldots,y_{m},y \right\rbrack \right)$$ $$+ X_{n}(x)f\left\lbrack x_{0},\ldots,x_{n},x;y \right\rbrack$$ $$= \sum_{i = 0}^{n}\sum_{j = 0}^{m}X_{i - 1}(x)Y_{j - 1}(y)f\left\lbrack x_{0},\ldots,x_{i};y_{0},\ldots,y_{j} \right\rbrack$$ $$+ Y_{m}(y)\sum_{i = 0}^{n}X_{i - 1}(x)f\left\lbrack x_{0},\ldots,x_{i};y_{0},\ldots,y_{m},y \right\rbrack + X_{n}(x)f\left\lbrack x_{0},\ldots,x_{n},x;y \right\rbrack$$
Let
$$p^\ast (x,y) = \sum_{i = 0}^{n}\sum_{j = 0}^{m}X_{i - 1}(x)Y_{j - 1}(y)f\left\lbrack x_{0},\ldots,x_{i};y_{0},\ldots,y_{j} \right\rbrack \in \Pi_{n,m}$$
$$r(x,y) = Y_{m}(y)\sum_{i = 0}^{n}X_{i - 1}(x)f\left\lbrack x_{0},\ldots,x_{i};y_{0},\ldots,y_{m},y \right\rbrack + X_{n}(x)f\left\lbrack x_{0},\ldots,x_{n},x;y \right\rbrack$$
So we have $f(x,y) = p^\ast (x,y) + r(x,y)$. And it is obvious that
$$X_{n}\left( x_{i} \right) = 0,\ Y_{m}\left( y_{j} \right) = 0\quad\forall 0 \leq i \leq n,\ 0 \leq j \leq m$$
So we have $r\left( x_{i},y_{j} \right) = 0\ \forall i,j$. I. e., $p^{\ast \left( x_{i},y_{j} \right)} = f\left( x_{i},y_{j} \right)\ \forall i,j$. By the uniqueness of the interpolator polynomial in $\Pi_{n,m}$ we have that $p(x,y) = p^\ast (x,y)$. So $f(x,y) - p(x,y) = r(x,y)$.
Now, if we consider $f\left\lbrack x;y_{0},\ldots,y_{m} \right\rbrack$, we know that
$$f\left\lbrack x;y_{0},\ldots,y_{m} \right\rbrack = \left( f\lbrack x; \cdot \rbrack \right)\left\lbrack y_{0},\ldots,y_{m},y \right\rbrack$$
Now, by the definition of divided differences
$$\left( f\left\lbrack \cdot ;y_{0},\ldots,y_{m},y \right\rbrack \right)\lbrack x\rbrack = \left( f\left\lbrack \cdot ;y_{0},\ldots,y_{m},y \right\rbrack \right)(x) = f\left\lbrack x;y_{0},\ldots,y_{m} \right\rbrack$$
So we can consider it as a function of $x$. Because it implicitely depends on $f \in C^{n + m + 2}$ it will be at least $C^{n + 1}$ so we can apply Newton formula
$$f\left\lbrack x;y_{0},\ldots,y_{m},y \right\rbrack = \sum_{i = 0}^{n}X_{i - 1}(x)f\left\lbrack x_{0},\ldots,x_{i};y_{0},\ldots,y_{m},y \right\rbrack + X_{n}(x)f\left\lbrack x_{0},\ldots,x_{n},x;y_{0},\ldots,y_{m},y \right\rbrack$$
Manipulating,
$$Y_{m}(y)f\left\lbrack x;y_{0},\ldots,y_{m},y \right\rbrack =$$ $$Y_{m}(y)\sum_{i = 0}^{n}X_{i - 1}(x)f\left\lbrack x_{0},\ldots,x_{i};y_{0},\ldots,y_{m},y \right\rbrack + Y_{m}(y)X_{n}(x)f\left\lbrack x_{0},\ldots,x_{n},x;y_{0},\ldots,y_{m},y \right\rbrack$$
$$ Y_{m}(y)\sum_{i = 0}^{n}X_{i - 1}(x)f\left\lbrack x_{0},\ldots,x_{i};y_{0},\ldots,y_{m},y \right\rbrack = $$ $$Y_{m}(y)f\left\lbrack x;y_{0},\ldots,y_{m},y \right\rbrack - Y_{m}(y)X_{n}(x)f\left\lbrack x_{0},\ldots,x_{n},x;y_{0},\ldots,y_{m},y \right\rbrack$$
We remember that
$$r(x,y) = Y_{m}(y)\sum_{i = 0}^{n}X_{i - 1}(x)f\left\lbrack x_{0},\ldots,x_{i};y_{0},\ldots,y_{m},y \right\rbrack + X_{n}(x)f\left\lbrack x_{0},\ldots,x_{n},x;y \right\rbrack$$
So if we substitute
$$r(x,y) = X_{n}(x)f\left\lbrack x_{0},\ldots,x_{n},x;y \right\rbrack + Y_{m}(y)f\left\lbrack x;y_{0},\ldots,y_{m},y \right\rbrack - X_{n}(x)Y_{m}(y)f\left\lbrack x_{0},\ldots,x_{n},x;y_{0},\ldots,y_{m},y \right\rbrack$$
Now,
$$f\left\lbrack x_{0},\ldots,x_{n},x;y \right\rbrack = \left( f\lbrack \cdot ,y\rbrack \right)\left\lbrack x_{0},\ldots,x_{n},x \right\rbrack,\quad f\left\lbrack x;y_{0},\ldots,y_{m},y \right\rbrack = \left( f(x, \cdot ) \right)\left\lbrack y_{0},\ldots,y_{m},y \right\rbrack$$
Applying mean value theorem for univariate divided differences, we have that exist $\xi \in \left( x_{0},x_{n} \right),\ \eta \in \left( y_{0},y_{m} \right)$ such that
$$f\left\lbrack x_{0},\ldots,x_{n},x;y \right\rbrack = \frac{1}{(n + 1)!} \cdot \frac{\partial^{n + 1}f}{\partial x^{n + 1}}(\xi,y)$$
$$f\left\lbrack x;y_{0},\ldots,y_{m},y \right\rbrack = \frac{1}{(m + 1)!} \cdot \frac{\partial^{m + 1}f}{\partial y^{m + 1}}(x,\eta)$$
Finally we consider $f\left\lbrack x_{0},\ldots,x_{n},x;y_{0},\ldots,y_{m},y \right\rbrack$.
On the one hand,
$$f\left\lbrack x_{0},\ldots,x_{n},x;y_{0},\ldots,y_{m},y \right\rbrack = \left( f\left\lbrack \cdot ;y_{0},\ldots,y_{m},y \right\rbrack \right)\left\lbrack x_{0},\ldots,x_{n},x \right\rbrack$$
Applying mean value theorem to $f\left\lbrack \cdot ;y_{0},\ldots,y_{m},y \right\rbrack$, which we know it is $C^{n + 1}$ because it implicitely depends on $f$, we'll have that exists $\xi\prime \in \left( x_{0},x_{n} \right)$ such that
$$\left( f\left\lbrack \cdot ;y_{0},\ldots,y_{m},y \right\rbrack \right)\left\lbrack x_{0},\ldots,x_{n},x \right\rbrack = \frac{1}{(n + 1)!} \cdot \frac{\partial^{n + 1}f\left\lbrack \cdot ;y_{0},\ldots,y_{m},y \right\rbrack}{\partial x^{n + 1}}(\xi\prime)$$
Now,
$$\frac{\partial^{n + 1}f\left\lbrack \cdot ;y_{0},\ldots,y_{m},y \right\rbrack}{\partial x^{n + 1}}(\xi\prime) = \frac{\partial^{n + 1}f\left\lbrack \xi\prime;y_{0},\ldots,y_{m},y \right\rbrack}{\partial x^{n + 1}} = \left( \frac{\partial^{n + 1}f\lbrack\xi\prime; \cdot \rbrack}{\partial x^{n + 1}} \right)\left\lbrack y_{0},\ldots,y_{m},y \right\rbrack$$
And $\frac{\partial^{n + 1}f\lbrack\xi\prime; \cdot \rbrack}{\partial x^{n + 1}}$ is an expression obtained derivating $n + 1$ times wrt $x$ the function $f \in C^{n + m + 2}$. So this function will be at least $C^{n + m + 2 - (n + 1)} = C^{m + 1}$ and we can apply the mean value theorem, so it exists $\eta\prime \in \left( y_{0},y_{m} \right)$ such that
$$\left( \frac{\partial^{n + 1}f\lbrack\xi\prime; \cdot \rbrack}{\partial x^{n + 1}} \right)\left\lbrack y_{0},\ldots,y_{m},y \right\rbrack = \frac{1}{(m + 1)!} \cdot \frac{\partial^{n + 1}\frac{\partial^{m + 1}f\lbrack\xi\prime;\eta\prime\rbrack}{\partial y^{m + 1}}}{\partial x^{n + 1}}$$
Now,
$$f\lbrack\xi\prime;\eta\prime\rbrack = \left( f( \cdot ,\eta\prime) \right)\lbrack\xi\prime\rbrack = \left( f( \cdot ,\eta\prime) \right)(\xi\prime) = f(\xi\prime,\eta\prime)$$
So
$$\frac{\partial^{n + 1}\frac{\partial^{m + 1}f\lbrack\xi\prime;\eta\prime\rbrack}{\partial y^{m + 1}}}{\partial x^{n + 1}} = \frac{\partial^{n + 1}\frac{\partial^{m + 1}f(\xi\prime,\eta\prime)}{\partial y^{m + 1}}}{\partial x^{n + 1}} = \frac{\partial^{n + m + 2}f}{\partial x^{n + 1}y^{m + 1}}(\xi\prime,\eta\prime)$$
So we have
$$f\left\lbrack x_{0},\ldots,x_{n},x;y_{0},\ldots,y_{m},y \right\rbrack = \frac{1}{(n + 1)!(m + 1)!} \cdot \frac{\partial^{n + m + 2}f}{\partial x^{n + 1}y^{m + 1}}(\xi\prime,\eta\prime)$$
Substituting in $r(x,y)$ we'll have
$$r(x,y) = \frac{1}{(n + 1)!} \cdot \frac{\partial^{n + 1}f}{\partial x^{n + 1}}(\xi,y) \cdot X_{n}(x) + \frac{1}{(m + 1)!} \cdot \frac{\partial^{m + 1}f}{\partial y^{m + 1}}(x,\eta) \cdot Y_{m}(y)$$ $$- \frac{1}{(n + 1)!(m + 1)!} \cdot \frac{\partial^{n + m + 2}f}{\partial x^{n + 1}y^{m + 1}}(\xi\prime,\eta\prime)X_{n}(x)Y_{m}(y)$$
As we wanted to show.