浅谈分块二元Hermite插值

引言

​ 起因是数值分析的一次大作业，在学习了分段三次一元$Hermite$插值并编程实践后，老师布置了如下两个题目：

要求：

​ 自己构造相应的插值函数是困难的，而且不好验证，于是先转向文献调研，再根据调研结果进行分析和编程

最初想法与文献调研

​ 由于题目给出了$f(x,y),\frac{\partial f}{\partial x}(x,y),\frac{\partial f}{\partial y}(x,y),\frac{{\partial }^{2}f}{\partial x\partial y}(x,y)$，并且需要考虑分块插值，一个很自然的想法是考虑$Hermite$插值，当然，$Hermite$插值可以看作是$Lagrange$和$Newton$插值的极限形式，这两种插值也可以考虑，而在文献查阅过程中，考虑$Lagrange$插值一般不考虑导数；而大部分文献似乎并没有将$Newton$插值作为单独的一类插值，而是“$Newton$形式的插值”[1]：
$$p(x,y)=\sum _{i=0}^m\sum _{j=0}^{n}f[x_0.\cdots ,x_i;y_0,\cdots ,y_j]\prod _{h=0}^{i-1}(x-x_h)\prod _{k=0}^{j-1}(x-x_k)$$
​ 但是[1]中提供的关于$Newton$插值的文献大部分是不能访问的，又或者是没有基于导数的。

​ 当然，其实最暴力的方法是线性方程组求导[2]，尽管题目不允许，但这种方法可以帮助我们判断插值方程是否唯一（对于第二题）；样条插值也是一种可以考虑的方法，但是如果$Hermite$插值可行且线性方程组的解唯一，就没必要构造样条插值了。

​ 针对题目，仅考虑二元$Hermite$插值，而事实上，在实际应用中，一般也考虑二元插值，最多考虑三元，二元插值也是研究的最多的[3]。

二元$Hermite$插值的两种“范式”

​ 回顾一元$Hermite$插值的构造过程（方便起见，只考虑一阶导），对于已知$f(x_0),f(x_1),\cdots f(x_n);f^{\prime }(x_0),f^{\prime }(x_1),\cdots f^{\prime }(x_n)$的情形，构造插值函数
$$h(x)=\sum _{i=0}^n{A}_i(x)f(x_i)+\sum _{i=0}^n{B}_i(x)f^{\prime }(x_i)$$
​ 使其满足$h(x_i)=f(x_i),h^{\prime }(x_i)=f^{\prime }(x_i),i=0.1,\cdots ,n$，可以构造
$$\begin{array}{rl}& A_i(x)=[1-2l_i^{\prime }(x_i)(x-x_i)]l_i^2(x)\\ & B_i(x)=(x-x_i)l_i^2(x)\end{array}$$
​ 其中，$l_i(x)$为$Lagrange$插值基函数

​ 以下介绍两种二元$Hermite$插值“范式”（作者自称），方便起见，按作者起名为$Ahlin$范式[4]和$Chawla-Jayarajan$范式[5]。前者的推导思想与一元情形类似，而后者是基于前者的改进

​ 在构造插值函数的过程中，两者都参考了$Spitzbart$[6]在已知任意阶导数的情况下一元$Hermite$插值的构造

$Spitzbart$的构造

​ 假设已知$f^{(k)}(x_j),j=0,1,\cdots ,n;k=0,1,\cdots ,r_j$，设
$$\begin{gathered} p_j(x)=(x-x_0{)}^{r_0+1}\cdots (x-x_n{)}^{r_n+1} \\ {g}_j(x)=[p_j(x){]}^{-1} \end{gathered}$$
​ 可以构造$n+\sum _{j=0}^n{r}_j$次多项式函数$H(x)$，满足
$$H^{(k)}(x_j)=f^{(k)}(x_j)$$
​ $H(x)$如下：
$$H(x)=\sum _{j=0}^n\sum _{k=0}^{r_j}{A}_{jk}{f}^{(k)}(x_j)$$
​ 其中
$$A_{jk}(x)=p_j(x)\frac{(x-x_j{)}^k}{k!}\sum _{t=0}^{r_j-k}\frac{1}{t!}{g}_i^{(t)}(x_j)(x-x_j{)}^t$$
​ 容易验证：
$$A_{j1}=(x-x_j)\frac{p_j(x)}{p_j(x_j)}=(x-x_j)l_j^2(x)$$

$$\begin{array}{rl}{A}_{j0}& =p_j(x)[g_j(x_j)+g_j^{\prime }(x_j)(x-x_j)]\\ & =p_j(x)\left[\frac{1}{p_j(x_j)}+\frac{p_j^{\prime }(x_j)}{[p_j(x_j){]}^2}(x-x_j)\right]\\ & =\left[1-\frac{p_j^{\prime }(x_j)}{p_j(x_j)}(x-x_j)\right]l_j^2(x)\\ & =[1-2l_i^{\prime }(x_i)(x-x_i)]l_i^2(x)\end{array}$$

​ 这与我们的构造结果是一致的

$Ahlin$范式

​ 已知$f^{(p,q)}(x_r,y_i),p,q=0,\cdots ,k-1;r,i=1,\cdots ,n$，假设插值函数唯一，记为$p(x,y)$，则
$$\frac{{\partial }^{p+q}}{\partial x^p\partial y^q}p(x_r,y_i)=f^{(p,q)}(x_r,y_i)$$
​ 具体的推导过程比较复杂，可以参照原文（因为我没完全看懂），作者先是在复平面上考虑插值余项，进而构造插值函数，最后考虑$x,y$一阶偏导和混合偏导的特殊情形，从一个非常复杂的式子得到一个漂亮的结果：
$$\begin{array}{rl}f(x,y)\approx p(x,y)=& \sum _{i=1}^n\sum _{r=1}^n{h}_r(x)g_i(y)f(x_r,y_i)+\sum _{i=1}^n\sum _{r=1}^n{h}_r(x){\tilde{g}}_i(y)\frac{\partial }{\partial y}f(x_r,y_i)\\ & +\sum _{i=1}^n\sum _{r=1}^n{\tilde{h}}_r(x)g_i(y)\frac{\partial }{\partial x}f(x_r,y_i)+\sum _{i=1}^n\sum _{r=1}^n{\tilde{h}}_r(x){\tilde{g}}_i(y)\frac{{\partial }^2}{\partial x\partial y}f(x_r,y_i)\end{array}$$
​ 其中，
$$\begin{array}{rl}& h_i(x)=[1-2l_r^{\prime }(x_r)(x-x_r)]l_r^2(x)\\ & {\tilde{h}}_i(x)=(x-x_r)l_r^2(x)\\ & g_i(y)=[1-2m_r^{\prime }(y_i)(y-y_i)]m_i^2(y)\\ & {\tilde{h}}_i(y)=(y-y_i)m_i^2(y)\end{array}$$
​ $l_r(x),m_i(y)$分别为关于$x,y$的$Lagrange$插值基函数

​ 上式与一元的情形是相似的，且容易验证，满足$\frac{{\partial }^{p+q}}{\partial x^p\partial y^q}p(x_r,y_i)=f^{(p,q)}(x_r,y_i),p,q=0,1$

$Chawla-Jayarajan$范式

​ $Chawla-Jayarajan$范式的推导过程与一元情形类似，对于$Ahlin$范式改进的地方在于，考虑了对于$x,y$已知条件不对称的情形

​ 已知$f^{(k,l)}(x_i,y_j),i=0,\cdots ,m,j=0,\cdots ,n;k=1,\cdots ,r_i,l=1,\cdots ,r_j$，假设插值函数唯一，记为$H_{M,N}(x,y)$，其中$M=m+\sum _{i=0}^m{r}_i,N=n+\sum _{j=0}^n{r}_j$，构造
$$H_{m,n}(x,y)=\sum _{i=0}^m\sum _{j=0}^n\sum _{k=0}^{r_i}\sum _{l=0}^{s_j}{A}_{ik}(x)B_{jl}(y)f^{(k,l)}(x_i,y_j)$$
​ 其中$A_{ik}(x),B_{jl}(y)$作者参照了$Spitzbart$的研究，这里不再赘述

​ 对于已知两对$x,y$的函数值、偏导和混合偏导的特殊情形，即$\frac{{\partial }^{k+l}}{\partial x^k\partial y^l}{H}_{2m-1,2n-1}(x_i,y_j)=f^{(k,l)}(x_i,y_j),i,j=0,1;k=0,1,\cdots ,m-1,l=0,1,\cdots ,n-1$，作者给出了以下公式：
$$\begin{array}{rl}{H}_{2m-1,2n-1}& =(x-x_0{)}^m(y-y_0{)}^n\sum _{k=0}^{m-1}\sum _{l=0}^{n-1}{A}_{kl}\frac{(x-x_1{)}^k}{k!}\frac{(y-y_1{)}^l}{l!}\\ & +(x-x_0{)}^m(y-y_1{)}^n\sum _{k=0}^{m-1}\sum _{l=0}^{n-1}{B}_{kl}\frac{(x-x_1{)}^k}{k!}\frac{(y-y_0{)}^l}{l!}\\ & \\ & +(x-x_1{)}^m(y-y_0{)}^n\sum _{k=0}^{m-1}\sum _{l=0}^{n-1}{C}_{kl}\frac{(x-x_0{)}^k}{k!}\frac{(y-y_1{)}^l}{l!}\\ & +(x-x_1{)}^m(y-y_1{)}^n\sum _{k=0}^{m-1}\sum _{l=0}^{n-1}{D}_{kl}\frac{(x-x_0{)}^k}{k!}\frac{(y-y_0{)}^l}{l!}\end{array}$$
​ 其中
$$A_{kl}={\left[\frac{{\partial }^{k+l}}{\partial x^k\partial y^l}\frac{f(x,y)}{(x-x_0{)}^m}(y-y_0{)}^n\right]}_{x-x_1,y=y_1}$$
​ $B_{kl},C_{kl},D_{kl}$类似

解题与编程

范式的选择

​ 对于上述两种范式，需要根据需求选择一个。从编程的难易程度而言，$Ahlin$范式是优于$Chawla-Jayarajan$范式的，因为$Ahlin$范式与一元的$Hermite$插值公式相似，并且方便定义函数和构造循环体；而对于$Chawla-Jayarajan$范式，比较头疼的是$A_{kl}$的计算，尽管对于特定的$k,l$，$A_{kl},B_{kl},C_{kl},D_{kl}$可以一起实现，但其中偏导的部分在程序中是难以实现的，对于每阶偏导需要单独手算然后根据手算的结果代入数据，就非常繁琐。而且，当插值函数唯一时，这两种范式的表达式理论上是统一的。

​ 于是我们在解题时选择$Ahlin$范式。

分块的实现

​ 回顾一元分段三次$Hermite$插值，每个子区间$[x_i,x_{i+1}]$上的$p_3(x)$可以写成
$$\begin{array}{rl}{p}_3(x)=& \left(1-2\frac{x-x_i}{x_i-x_{i+1}}\right){\left(\frac{x-x_{i+1}}{x_i-x_{i+1}}\right)}^{2}f(x_i)+\left(1-2\frac{x_{i+1}-x}{x_i-x_{i+1}}\right){\left(\frac{x-x_i}{x_i-x_{i+1}}\right)}^{2}f(x_{i+1})\\ & +(x-x_i){\left(\frac{x-x_i}{x_i-x_{i+1}}\right)}^2{f}^{\prime }(x_i)+(x-x_{i+1}){\left(\frac{x-x_i}{x_{i+1}-x_i}\right)}^2{f}^{\prime }(x_{i+1})\end{array}$$
​ 于是，类似地，可以构造二元分块$3\times 3$次$Hermite$插值，为方便起见，对于$x$引入插值基函数：
$$\begin{gathered} {\varphi }_0(x)=\{\begin{array}{ll}\left(1-2\frac{x-x_0}{x_0-x_1}\right){\left(\frac{x-x_1}{x_0-x_1}\right)}^2,& x\in [x_0,x_1]\\ 0,& else\end{array} \\ {\varphi }_i(x)=\{\begin{array}{ll}\left(1-2\frac{x-x_i}{x_i-x_{i+1}}\right){\left(\frac{x-x_{i+1}}{x_i-x_{i+1}}\right)}^2,& x\in [x_{i-1},x_i)\\ \left(1-2\frac{x_{i+1}-x}{x_i-x_{i+1}}\right){\left(\frac{x-x_i}{x_i-x_{i+1}}\right)}^2,& x\in [x_i,x_{i+1}]\\ 0,& else\end{array} \\ {\varphi }_n(x)=\{\begin{array}{ll}\left(1-2\frac{x_n-x}{x_{n-1}-x_n}\right){\left(\frac{x-x_{n-1}}{x_{n-1}-x_n}\right)}^2,& x\in [x_{n-1},x_n]\\ 0,& else\end{array} \\ \widetilde{{\varphi }_0}(x)=\{\begin{array}{ll}(x-x_0){\left(\frac{x-x_0}{x_i-x_1}\right)}^2,& x\in [x_0,x_1]\\ 0,& else\end{array} \\ \widetilde{{\varphi }_i}(x)=\{\begin{array}{ll}(x-x_i){\left(\frac{x-x_i}{x_i-x_{i+1}}\right)}^2,& x\in [x_{i-1},x_i)\\ (x-x_{i+1}){\left(\frac{x-x_i}{x_{i+1}-x_i}\right)}^{2}2,& x\in [x_i,x_{i+1}]\\ 0,& else\end{array} \\ \widetilde{{\varphi }_n}(x)=\{\begin{array}{ll}(x-x_n){\left(\frac{x-x_{n-1}}{x_n-x_{n-1}}\right)}^2,& x\in [x_{n-1},x_n]\\ 0,& else\end{array} \end{gathered}$$
​ 对于$y$：
$$\begin{gathered} {\psi }_0(y)=\{\begin{array}{ll}\left(1-2\frac{y-y_0}{y_0-y_1}\right){\left(\frac{y-y_1}{y_0-y_1}\right)}^2,& y\in [y_0,y_1]\\ 0,& else\end{array} \\ {\psi }_i(y)=\{\begin{array}{ll}\left(1-2\frac{y-y_i}{y_i-y_{i+1}}\right){\left(\frac{y-y_{i+1}}{y_i-y_{i+1}}\right)}^2,& y\in [y_{i-1},y_i)\\ \left(1-2\frac{y_{i+1}-y}{y_i-y_{i+1}}\right){\left(\frac{y-y_i}{y_i-y_{i+1}}\right)}^2,& y\in [y_i,y_{i+1}]\\ 0,& else\end{array} \\ {\psi }_n(y)=\{\begin{array}{ll}\left(1-2\frac{y_n-y}{y_{n-1}-y_n}\right){\left(\frac{y-y_{n-1}}{y_{n-1}-y_n}\right)}^2,& y\in [y_{n-1},y_n]\\ 0,& else\end{array} \\ \widetilde{{\psi }_0}(y)=\{\begin{array}{ll}(y-y_0){\left(\frac{y-y_0}{y_i-y_1}\right)}^2,& y\in [y_0,y_1]\\ 0,& else\end{array} \\ \widetilde{{\psi }_i}(y)=\{\begin{array}{ll}(y-y_i){\left(\frac{y-y_i}{y_i-y_{i+1}}\right)}^2,& y\in [y_{i-1},y_i)\\ (y-y_{i+1}){\left(\frac{y-y_i}{y_{i+1}-y_i}\right)}^{2}2,& y\in [y_i,y_{i+1}]\\ 0,& else\end{array} \\ \widetilde{{\psi }_n}(y)=\{\begin{array}{ll}(y-y_n){\left(\frac{y-y_{n-1}}{y_n-y_{n-1}}\right)}^2,& y\in [y_{n-1},y_n]\\ 0,& else\end{array} \end{gathered}$$
​ 于是，每个子块$[x_i,x_{i+1}]\times [y_j,y_{j+1}]$上的$p_{3,3}(x,y)$可以写成
$$\begin{array}{rl}{p}_{3.3}(x,y)=& \sum _{j=1}^n\sum _{i=1}^n{\varphi }_i(x)g_j(y)f(x_i,y_j)+\sum _{j=1}^n\sum _{i=1}^n{\varphi }_i(x){\tilde{\psi }}_j(y)\frac{\partial }{\partial y}f(x_i,y_j)\\ & +\sum _{j=1}^n\sum _{i=1}^n{\tilde{\varphi }}_i(x)g_j(y)\frac{\partial }{\partial x}f(x_i,y_j)+\sum _{j=1}^n\sum _{i=1}^n{\tilde{\varphi }}_i(x){\tilde{\psi }}_j(y)\frac{{\partial }^2}{\partial x\partial y}f(x_i,y_j)\end{array}\cdots (\ast )$$

问题1的解决

​ 根据上述分析，算法如下：

	用$C/C++$进行编程实现，顺便对于区间的选择和分块插值次数做了泛化，代码如下：

#include <stdio.h>
#include <math.h>
#include <stdlib.h>

void data(double *x, double *y, double *f, double *dx, double *dy, double *dxy, int n) {
    // 生成sin(x^2y+1)的dx,dy,dxy
    int i, j;

    for (i = 0; i <= n; i++) {
        for (j = 0; j <= n; j++) {
            f[i * n + j] = sin(x[i] * x[i] * y[j] + 1);
            dx[i * n + j] = 2 * x[i] * y[j] * cos(x[i] * x[i] * y[j] + 1);
            dy[i * n + j] = x[i] * x[i] * cos(x[i] * x[i] * y[j] + 1);
            dxy[i * n + j] = 2 * x[i] * cos(x[i] * x[i] * y[j] + 1) - 2 * pow(x[i], 3) * y[i] * sin(x[i] * x[i] * y[j] + 1);
        }
    }
}

double lagrange(double *z, double zs, int n, int i) {
    // 求lagrange基函数
    int j;
    double t = 1;

    for (j = 0; j < n; j++) {
        if (i != j) {
            t *= ((zs - z[j]) / (z[i] - z[j]));
        }
    }
    return t;
}

double lagrange_(double *z, int i, int n) {
    // 求lagrange基函数导数
    int j;
    double r = 0;

    for (j = 0; j < n; j++) {
        if (i != j) {
            r += (1 / (z[i] - z[j]));
        }
    }
    return r;
}

double Hermite(double *x, double *y, double *f, double *dx, double *dy, double *dxy, int m, double xs, double ys, int s) {
    int i, j, k, t;
    double r = 0, h, h_, g, g_;
    double *x_, *y_;


    x_ = (double *) malloc(sizeof(double) * (s - 1));
    y_ = (double *) malloc(sizeof(double) * (s - 1));

    for (i = 0; i < m; i++) {
        for (j = 0; j < m; j++) {
            if (x[i] <= xs && xs <= x[i + s - 2] && y[j] <= ys && ys <= y[j + s - 2] && i % (s - 2) == 0 && j % (s - 2) == 0 && i + s - 2 <= m && j + s - 2 <= m) {
                for (k = 0; k < s - 1; k++) {
                    x_[k] = x[i + k];
                    y_[k] = y[j + k];
                }

                for (k = 0; k <= s - 2; k++) {
                    for (t = 0; t <= s - 2; t++) {
                        h = (1 - 2 * lagrange_(x_, k, s - 1) * (xs - x_[k])) * lagrange(x_, xs, s - 1, k) * lagrange(x_, xs, s - 1, k);
                        h_ = (xs - x_[k]) * lagrange(x_, xs, s - 1, k) * lagrange(x_, xs, s - 1, k);
                        g = (1 - 2 * lagrange_(y_, t, s - 1) * (ys - y_[t])) * lagrange(y_, ys, s - 1, t) * lagrange(y_, ys, s - 1, t);
                        g_ = (ys - y_[t]) * lagrange(y_, ys, s - 1, t) * lagrange(y_, ys, s - 1, t);
                        r += (h * g * f[(i + k) * m + (j + t)] + h * g_ * dy[(i + k) * m + (j + t)] + h_ * g * dx[(i + k) * m + (j + t)] + h_ * g_ * dxy[(i + k) * m + (j + t)]);
                    }
                }
            }
        }
    }

    free(x_);
    free(y_);

    return r;
}

int main() {
    int n = 5; // 区间n等分
    int s = 3; // s次插值
    double *x, *y, *dx, *dy, *dxy, *f;
    double xs, ys;
    double left = 0, right = 1;
    int i, j;
    double value, result;

    xs = 1.0 / 3.0;
    ys = 2.0 / 3.0;
    x = (double *) malloc(sizeof(double) * (n + 1));
    dx = (double *) malloc(sizeof(double) * (n + 1) * (n + 1));
    y = (double *) malloc(sizeof(double) * (n + 1));
    dy = (double *) malloc(sizeof(double) * (n + 1) * (n + 1));
    dxy = (double *) malloc(sizeof(double) * (n + 1) * (n + 1));
    f = (double *) malloc(sizeof(double) * (n + 1) * (n + 1));

    // 等分区间
    for (i = 0; i <= n; i++) {
        x[i] = ((right - left) / n) * i + left;
        y[i] = ((right - left) / n) * i + left;
    }
    data(x, y, f, dx, dy, dxy, n);

    value = sin(xs * xs * ys + 1);
    printf("Actual Value: %.15lf\n", value);
    result = Hermite(x, y, f, dx, dy, dxy, n, xs, ys, s);
    printf("Hermite Value: %.15lf\n", result);
    printf("error: %.15lf", abs(result - value));

    free(x);
    free(dx);
    free(y);
    free(dy);
    free(f);
    free(dxy);

    return 0;
}

​ 结果为：

问题2的解决

计算插值多项式在特定点的值

​ 如果要满足题目所需条件，即满足$p(x,y)=f(x,y),\frac{\partial }{\partial x}p(x,y)=\frac{\partial }{\partial x}f(x,y),\frac{\partial }{\partial y}p(x,y)=\frac{\partial }{\partial y}f(x,y)$，显然，$(\ast )$式是满足条件的，因此只需将上述代码中Hermite函数中的

r += (h * g * f[(i + k) * m + (j + t)] + h * g_ * dy[(i + k) * m + (j + t)] + h_ * g * dx[(i + k) * m + (j + t)] + h_ * g_ * dxy[(i + k) * m + (j + t)]);

改成

r += (h * g * f[(i + k) * m + (j + t)] + h * g_ * dy[(i + k) * m + (j + t)] + h_ * g * dx[(i + k) * m + (j + t)]);

即可

​ 结果为：

连续性和可微性的讨论

​ 对于连续性，在$(x_i,x_{i+1}),i=0,1,\cdots ,n-1$，由于$p_3(x)$为多项式函数，显然连续；而在点$x=x_i,i=0,1,\cdots ,n$处，$p_3(x)=0$，显然连续。

​ 综上，$p_3(x)$在$[x_0,x_n]$上连续。

​ 对于可微性，由于$\frac{{\partial }^2}{\partial x\partial y}f(x,y)$没有数据（实际上是按0处理了），则$\forall x\in [x_0,x_n]$，$\frac{{\partial }^2}{\partial x\partial y}{p}_3(x,y)\equiv 0$；由于$p_3(x)$为多项式函数，显然其关于$x,y$的偏导数存在，故关于$x,y$可微。

​ 综上，$p_3(x)$在$[x_0,x_n]$上可微。

对于问题2的疑惑

​ 从上述对于问题2的讨论中可以发现，由于$\frac{{\partial }^2}{\partial x\partial y}f(x,y)$没有数据，我们的处理是直接将$\frac{{\partial }^2}{\partial x\partial y}f(x,y)$设置成0，这就导致了我们的误差在问题1结果的基础上又加上了
$$\sum _{j=1}^n\sum _{i=1}^n{\tilde{\varphi }}_i(x){\tilde{\psi }}_j(y)\frac{{\partial }^2}{\partial x\partial y}f(x_i,y_j)$$
​ 这直接导致误差大了十倍甚至九倍，是否合理呢？

​ 首先，我们的已知条件由16个缩减到了12个，于是我们只能构造12个方程，因此我们无法构造$p(x,y)=\sum _{i=0}^3\sum _{j=0}^3{c}_{ij}{x}^i{y}^j$，因为有16个未知数。我们只能构造$p(x)=\sum _{i=0}^2\sum _{j=0}^2{c}_{ij}{x}^i{y}^j+q(x,y)$，其中$q(x,y)$包含$x^k{y}^{3-k},k=0,1,2,3$中的三项并乘一个常数，因此插值函数是不唯一的

​ 而在Tomas Sauer和Yuan Xu的研究中[7]，对于多元$Hermite$插值的不同插值条件，即已知不同方向导数的情形，做了进一步的讨论，其引入了适定性(posedness)的概念（事实上，现代对于插值的大部分研究都会引入并考虑这一概念），对于类似问题2的情形做了进一步的阐释

​ 总之，由于数据的缺失，误差是无法避免的

误差分析

​ 回顾一元分段三次$Hermite$插值的误差分析，其余项为：
$$R(x)=\frac{f^{(4)}({\xi }_i)}{4!}(x-x_i{)}^2(x-x_{i+1}{)}^2,x_i<{\xi }_i<x_{i+1}$$
​ 则其最大误差限为：
$$|R(x)|=\left|\frac{f^{(4)}({\xi }_i)}{4!}(x-x_i{)}^2(x-x_{i+1}{)}^2\right|\le \frac{h^4}{384}\underset{x_i<x<x_{i+1}}{\max }|f^{(4)}(x)|,x_i<{\xi }_i<x_{i+1}$$
​ 对于$Ahlin$范式和$Chawla-Jayarajan$范式，两个范式的作者分别给出了各自的误差分析

​ $Ahlin$是在复平面中分析的，考虑包含$x,x_0,x_1,\cdots ,x_{n-1}$的$Jordan$曲线$C_1$和$y,y_0,y_1,\cdots ,y_{n-1}$的$Jordan$曲线，且$k-1$为已知的最高次偏导，则插值余项为：
$$\begin{array}{rl}f(x,y)-p(x,y)=& \frac{[\lambda (x){]}^k}{2\pi i}{\int }_{C_1}\frac{f(s,y)}{[\lambda (s){]}^k(s-x)}ds+\frac{[\mu (y){]}^k}{2\pi i}{\int }_{C_2}\frac{f(x,t)}{[\mu (s){]}^k(t-y)}dt\\ & -\frac{[\lambda (x)\mu (y){]}^k}{2\pi i}{\int }_{C_1}{\int }_{C_2}\frac{f(s,t)}{[\lambda (s)\mu (t){]}^k(s-x)(t-y)}dsdt\end{array}$$
​ 其中
$$\begin{gathered} \lambda (x)=\prod _{i=0}^{n-1}(x-x_i) \\ \mu (y)=\prod _{i=0}^{n-1}(y-y_i) \end{gathered}$$
​ $Chawla-Jayarajan$范式的余项为：
$$\begin{array}{rl}{R}_{M,N}(x,y)=& \frac{\alpha (x)}{(M+1)!}\frac{{\partial }^{M+1}}{\partial x^{M+1}}f(\xi ,y)+\frac{\beta (y)}{(N+1)!}\frac{{\partial }^{N+1}}{\partial y^{N+1}}f(x,\eta )\\ & -\frac{\alpha (x)\beta (y)}{(M+1)!(N+1)!}\frac{{\partial }^{M+N+2}}{\partial x^{M+1}\partial y^{N+1}}f(\xi ,\eta ),\\ & x_i<\xi <x_{i+1},y_j<\eta <y_{j+1}\end{array}$$
​ 其中
$$\begin{gathered} \alpha (x)=(x-x_0{)}^{r_0+1}\cdots (x-x_m{)}^{r_m+1} \\ \beta (y)=(y-y_0{)}^{s_0+1}\cdots (y-y_m{)}^{s_n+1} \end{gathered}$$

​ 由于$Ahlin$范式的误差分析是在复平面上分析，需要用到留数定理进行误差分析；而$Chawla-Jayarajan$范式的余项相对简洁的多，而且和一元$Hermite$插值的余项更加相似。

​ 我们只考虑插值公式唯一的情形（问题1），对于不唯一的情形（问题2）可以在插值公式唯一的情形加上或减去缺失的项。

​ 由于插值公式唯一，所以插值余项理论上也是唯一的，于是我们选择更为简洁的$Chawla-Jayarajan$范式的余项进行误差分析。

​ 则问题1的余项为：
$$\begin{array}{rl}{R}_{3,3}(x,y)=& \frac{(x-0.2{)}^2(x-0.4{)}^2}{4!}\frac{{\partial }^4}{\partial x^4}f(\xi ,y)+\frac{(y-0.6{)}^2(y-0.8{)}^2}{4!}\frac{{\partial }^4}{\partial y^4}f(x,\eta )\\ & -\frac{(x-0.2{)}^2(x-0.4{)}^2(y-0.6{)}^2(y-0.8{)}^2}{4!4!}\frac{{\partial }^8}{\partial x^4\partial y^4}f(\xi ,\eta )\end{array}$$

​ 则最大误差限为：
$$\begin{array}{rl}|R_{3,3}(x,y)|& \le \left|\frac{0.2^4}{384}\frac{{\partial }^4}{\partial x^4}f(\xi ,y)\right|+\left|\frac{0.2^4}{384}\frac{{\partial }^4}{\partial x^4}f(x,\eta )\right|+\left|\frac{0.2^8}{384^2}\frac{{\partial }^8}{\partial x^4\partial y^4}f(\xi ,\eta )\right|\end{array}$$
​ 对于$f(x,y)=sin(x^{2}y+1),0.2<x<0.4,0.6<y<0.8$，利用$Sagemath$求导可得：

sage: x, y = var('x,y')
sage: f = sin(x^2*y+1)
sage: f.diff(x, 4)
16*x^4*y^4*sin(x^2*y + 1) - 48*x^2*y^3*cos(x^2*y + 1) - 12*y^2*sin(x^2*y + 1)
sage: f.diff(y, 4)
x^8*sin(x^2*y + 1)
sage: f.diff(x, 4, y, 4)
16*x^12*y^4*sin(x^2*y + 1) - 304*x^10*y^3*cos(x^2*y + 1) - 1740*x^8*y^2*sin(x^2*y + 1) + 3360*x^6*y*cos(x^2*y + 1) + 1680*x^4*sin(x^2*y + 1)

​ 代入得：
$$\begin{array}{rl}|R_{3,3}(x,y)|& \le \left|\frac{0.2^4}{384}[16{\xi }^4{y}^{4}sin({\xi }^{2}y+1)-48{\xi }^2{y}^{3}cos({\xi }^{2}y+1)-12y^{2}sin({\xi }^{2}y+1)]\right|\\ & +\left|\frac{0.2^4}{384}{x}^{8}sin(x^2\eta +1)\right|+|\frac{0.2^8}{384^2}[16{\xi }^{12}{\eta }^{4}sin({\xi }^2\eta +1)-304{\xi }^{10}{\eta }^{3}cos({\xi }^2\eta +1)\\ & -1740{\xi }^8{\eta }^{2}sin({\xi }^2\eta +1)+3360{\xi }^6\eta cos({\xi }^2\eta +1)+1680{\xi }^{4}sin({\xi }^2\eta +1)]|\\ & \le 3.806762660674352\times 10^{-5}\end{array}$$

参考文献

[1] M. Gasca and T. Sauer, “On the history of multivariate polynomial interpolation,” J. Comput. Appl. Math., vol. 122, no. 1–2, pp. 23–35, 2000.

[2] Wikipedia contributors. (2022, September 21). Bicubic interpolation. In Wikipedia, The Free Encyclopedia. Retrieved 08:11, October 9, 2022, from https://en.wikipedia.org/w/index.php?title=Bicubic_interpolation&oldid=1111555946.

[3] Wikipedia contributors. (2022, August 14). Multivariate interpolation. In Wikipedia, The Free Encyclopedia. Retrieved 08:20, October 9, 2022, from https://en.wikipedia.org/w/index.php?title=Multivariate_interpolation&oldid=1104415950.

[4] A. C. Ahlin, ”A bivariate generalization of Hermite’s interpolation formula”, Math. Comp.
18 (1964), 264-273.

[5] M. M. Chawla and N. Jayarajan, “A generalization of Hermite’s interpolation formula in two variables”, J. Austral. Math. Soc., vol. 18, no. 4, pp. 402-410, Dec. 1974.

[6] A. Spitzbart, “A generalization of Hermite’s interpolation formula”, Amer. Math. Monthly 67 (1960), 42-6.

[7] T. Sauer and Yuan Xu, “On multivariate Lagrange interpolation”, Math. Comp. 64 (1995) 1147–1170.