计算几何模型 — 从多项式插值到散乱数据拟合

主讲人：复旦大学数学科学学院 曹沅

---

插值问题的提出

问题 1（型值点插值）. 已知平面上的有序点列

$$P_i(x_i,y_i),i=0,1,\cdots ,n$$

假设

$$x_0<x_1<\cdots <x_{n-1}<x_n$$

求平面上的连续曲线 $C：y=f(x)$，要求 $C$ 通过 $P_i$，即 $y_i=f(x_i)$。

这是函数插值问题：已知的信息只是有限个点，要求补全为一条整体的连续曲线。

问题 2（函数逼近意义下的插值收敛性）. 已知平面上的曲线函数 $f(x)$，$x\in [a,b]$，在 $f(x)$ 上取有序点列 $\{P_i(x_i,y_i){\}}_{i=0}^n$，求平面上的连续曲线函数 $f_n(x)$，使 $f_n(x)$ 按序通过点 $P_i(x_i,y_i)$，且满足

$$\underset{n\to \infty }{\lim }\Vert f(x)-f_n(x){\Vert }_{\infty }=0$$

也就是说，我们不仅要求插值曲线通过给定点，还希望当采样点越来越密时，插值曲线能够一致收敛到原始曲线 $f(x)$。这是衡量一个插值方法"好坏"的核心标准。

---

多项式插值

直接待定系数法

最自然的想法是用一个 $n$ 次多项式函数去拟合 $n+1$ 个点：

$$y=a_0+a_{1}x+\cdots +a_{n-1}{x}^{n-1}+a_n{x}^n$$

由插值条件

$$y_i=f(x_i),i=0,1,\cdots ,n$$

得到关于系数 $a_0,a_1,\cdots ,a_n$ 的线性方程组：

$$\{\begin{array}{l}{y}_0=a_0+a_1{x}_0+\cdots +a_{n-1}{x}_0^{n-1}+a_n{x}_0^n\\ ⋮\\ y_n=a_0+a_1{x}_n+\cdots +a_{n-1}{x}_n^{n-1}+a_n{x}_n^n\end{array}$$

写成矩阵形式：

$$\left(\begin{array}{ccccc}1& x_0& \cdots & x_0^{n-1}& x_0^n\\ ⋮& & & & ⋮\\ 1& x_n& \cdots & x_n^{n-1}& x_n^n\end{array}\right)\left(\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_{n-1}\\ a_n\end{array}\right)=\left(\begin{array}{c}{y}_0\\ y_1\\ ⋮\\ y_{n-1}\\ y_n\end{array}\right)$$

系数矩阵是 Vandermonde 矩阵，其行列式为

$$\det =\prod _{0\le i<j\le n}(x_j-x_i)$$

当 $x_0<x_1<\cdots <x_n$（即节点两两不同）时，行列式非零，因而线性方程组有唯一解。解出 $a_i,i=0,1,\cdots ,n$，便得到所要求的 $n$ 次插值多项式 $y=f(x)$。

这种方法概念简单，但需要求解一个 $(n+1)\times (n+1)$ 的线性方程组，当 $n$ 较大时计算量大、数值不稳定（Vandermonde 矩阵条件数很差）。因此人们更常用下面的 Lagrange 显式表达式。

Lagrange 插值公式

定理（Lagrange 插值）. 通过点 $\{(x_i,y_i){\}}_{i=0}^n$（节点两两不同）的 $n$ 次插值多项式可以直接写为

$$y=\sum _{i=0}^n{y}_i{l}_i(x)$$

其中 $l_i(x)$ 称为 Lagrange 基函数：

$$l_i(x)=\prod _{\begin{array}{c}j=0\\ j\ne i\end{array}}^n\frac{x-x_j}{x_i-x_j}=\frac{(x-x_0)\cdots (x-x_{i-1})(x-x_{i+1})\cdots (x-x_n)}{(x_i-x_0)\cdots (x_i-x_{i-1})(x_i-x_{i+1})\cdots (x_i-x_n)}$$

验证插值性质： 将 $x_j$ 代入 $l_i(x)$，可得

$$l_i(x_j)={\delta }_{ij}=\{\begin{array}{ll}1,& i=j\\ 0,& i\ne j\end{array}$$

这是因为当 $j\ne i$ 时，$l_i(x)$ 的分子中含有因子 $(x-x_j)$，代入 $x=x_j$ 后该因子为零；而当 $j=i$ 时，分子分母完全相同，故比值为 1。由此立即验证

$$y(x_i)=\sum _{k=0}^n{y}_k{l}_k(x_i)=\sum _{k=0}^n{y}_k{\delta }_{ki}=y_i$$

对任意 $x\in [x_0,x_n]$，插值曲线 $y=f(x)$ 均可用此公式直接计算，无需求解线性方程组。

Lagrange 插值的缺陷：Runge 现象

Lagrange 插值虽然表达式简洁，但存在明显缺点：

计算不稳定： 节点数增多时多项式次数随之升高，插值曲线容易发生剧烈振荡；这就是所谓的 Runge 现象。

Runge 在研究多项式插值时发现，节点取得越多，插值多项式未必越精确。著名的反例是

$$f(x)=\frac{1}{1+25x^2},x\in [-1,1]$$

当采用等距节点做高次多项式插值时，插值函数在区间两端会发生剧烈的振荡，误差随次数增高反而增大。其根本原因在于高次多项式的系数对数据的微小扰动（包括舍入误差）极其敏感，且等距节点的 Lebesgue 常数随 $n$ 指数增长。

插值收敛性的数学提法

为了得到稳定、可控的插值方法，我们把"好的插值方法"用如下数学语言精确表述：

设 $f(x)\in C[a,b]$ 是任一被插值函数。取节点 $x_i\in [a,b],i=1,\cdots ,n$。求 $s(x)\in C[a,b]$，使得

$$s(x_i)=f(x_i),i=1,\cdots ,n$$

且当 $n\to \infty$（节点越取越多、越取越密）时，

$$\Vert f(x)-s(x){\Vert }_{\infty }\to 0$$

进一步地，如果 $f(x)\in C^l[a,b]$（$l=0,1,\cdots ,d$，即 $f$ 具有到 $d$ 阶的连续导数），我们希望 $s(x)\in C^r[a,b]$（$r=0,1,\cdots ,h$，$h<d$），并且各阶导数也同时一致收敛：

$$\Vert f^{(r)}(x)-s^{(r)}(x){\Vert }_{\infty }\to 0,r=0,1,\cdots ,h$$

样条函数的产生动机： 高次多项式插值虽然连续性好（$C^{\infty }$），但数值不稳定；要消除 Runge 现象、获得数值稳定性，一个自然的想法是适当降低整体的光滑性（连续阶数），换取数值上的稳定性——即不用单一的高次多项式，而是用分段低次多项式并在分段点处满足有限阶的连续性条件。这就是样条函数的基本思想。

---

三次样条插值

物理背景

样条（spline）一词来自制图工艺中使用的富有弹性的木条或钢条（"压铁样条"）。将其固定在若干型值点上任其自然弯曲，得到的曲线在工程制图中被认为是最"光顺"的曲线。根据对木样条的力学分析（细梁的小挠度弯曲方程），在忽略高阶小量的假设下，木样条在两个约束点之间近似满足三次多项式关系。这就是三次样条函数的物理来源。

记号与问题的提出

设分割

$$x_0<x_1<\cdots <x_{n-1}<x_n$$

• 节点： $\{x_i{\}}_{i=0}^n$；
• 内节点： $\{x_i{\}}_{i=1}^{n-1}$；
• $n$ 个小区间： $[x_{i-1},x_i],i=1,2,\cdots ,n$，记 $h_i=x_i-x_{i-1}$。

目标： 寻找 $[x_0,x_n]$ 上的函数 $y=s(x)$，使得

1. $s(x)$ 在每个小区间 $[x_{i-1},x_i]$ 上都是次数不超过 3 的多项式；
2. $s(x)$ 在所有内节点处达到 $C^2$ 连续（即 $s,s^{\prime },s^{\prime \prime }$ 均连续）；
3. $s(x_i)=y_i,i=0,1,\cdots ,n$（插值条件）。

这样的 $s(x)$ 就称为三次样条插值函数。

用二阶导数 $M_i$ 表示 $s(x)$

设 $s(x)$ 在节点 $\{x_i{\}}_{i=0}^n$ 处的二阶导数为

$$s^{\prime \prime }(x_i)=M_i,i=0,1,\cdots ,n$$

在小区间 $[x_{i-1},x_i]$ 上，已知端点数据：

$$x_{i-1},y_{i-1},M_{i-1};x_i,y_i,M_i$$

共 4 个条件（两个函数值、两个二阶导数值），恰好可以唯一确定一个 3 次多项式。

第一步：由 $M_{i-1},M_i$ 得到 $s^{\prime \prime }(x)$。

因为 $s(x)$ 在 $[x_{i-1},x_i]$ 上是三次多项式，所以 $s^{\prime \prime }(x)$ 在该区间上是关于 $x$ 的线性函数。利用两点式，满足 $s^{\prime \prime }(x_{i-1})=M_{i-1}$，$s^{\prime \prime }(x_i)=M_i$ 的线性函数为

$$s^{\prime \prime }(x)=M_{i-1}\frac{x_i-x}{h_i}+M_i\frac{x-x_{i-1}}{h_i},x\in [x_{i-1},x_i]$$

第二步：积分两次得到 $s(x)$。

对上式积分一次：

$$s^{\prime }(x)=-M_{i-1}\frac{(x_i-x{)}^2}{2h_i}+M_i\frac{(x-x_{i-1}{)}^2}{2h_i}+c$$

再积分一次，并把积分常数写成有对称结构的形式（即用两个待定常数代替一个）：

$$s(x)=M_{i-1}\frac{(x_i-x{)}^3}{6h_i}+M_i\frac{(x-x_{i-1}{)}^3}{6h_i}+c_{i-1}(x_i-x)+c_i(x-x_{i-1})$$

这里之所以能够用 $c_{i-1}(x_i-x)+c_i(x-x_{i-1})$ 代替一次多项式 $cx+d$，是因为 $\{x_i-x,x-x_{i-1}\}$ 与 $\{1,x\}$ 张成同一个一次多项式空间，这样写是为了后面代入端点条件时形式对称、简洁。

第三步：由插值条件确定 $c_{i-1},c_i$。

代入插值条件 $s(x_{i-1})=y_{i-1}$，$s(x_i)=y_i$：

当 $x=x_{i-1}$ 时，$(x_i-x_{i-1})=h_i$，$(x-x_{i-1})=0$，故

$$y_{i-1}=M_{i-1}\frac{h_i^3}{6h_i}+c_{i-1}{h}_i=\frac{M_{i-1}{h}_i^2}{6}+c_{i-1}{h}_i⟹c_{i-1}=\frac{y_{i-1}}{h_i}-\frac{h_i{M}_{i-1}}{6}$$

同理，当 $x=x_i$ 时，

$$c_i=\frac{y_i}{h_i}-\frac{h_i{M}_i}{6}$$

第四步：$s(x)$ 的显式表达式。

$$s(x)=M_{i-1}\frac{(x_i-x{)}^3}{6h_i}+M_i\frac{(x-x_{i-1}{)}^3}{6h_i}+\left(\frac{y_{i-1}}{h_i}-\frac{h_i{M}_{i-1}}{6}\right)(x_i-x)+\left(\frac{y_i}{h_i}-\frac{h_i{M}_i}{6}\right)(x-x_{i-1}),$$

$$x\in [x_{i-1},x_i]$$

对其求导，得到一阶导数：

$$s^{\prime }(x)=-M_{i-1}\frac{(x_i-x{)}^2}{2h_i}+M_i\frac{(x-x_{i-1}{)}^2}{2h_i}+\frac{y_i-y_{i-1}}{h_i}-\frac{h_i(M_i-M_{i-1})}{6}$$

矩阵形式： 上述表达式可以整理为幂基下的矩阵乘积形式，便于编程实现：

$$s(x)=\left(\begin{array}{cccc}1,& (x-x_{i-1}),& (x-x_{i-1}{)}^2,& (x-x_{i-1}{)}^3\end{array}\right)\left(\begin{array}{cccc}1& 0& 0& 0\\ \frac{-1}{h_i}& \frac{1}{h_i}& \frac{-h_i}{3}& \frac{-h_i}{6}\\ 0& 0& \frac{1}{2}& 0\\ 0& 0& \frac{-1}{6h_i}& \frac{1}{6h_i}\end{array}\right)\left(\begin{array}{c}{y}_{i-1}\\ y_i\\ M_{i-1}\\ M_i\end{array}\right),x\in [x_{i-1},x_i]$$

到这里，$s(x)$ 的表达式已经完全由 $\{y_i\}$（已知）和 $\{M_i\}$（未知）确定。现在问题归结为如何求得 $M_i,i=0,1,\cdots ,n$。

连续性方程

以上构造保证了 $s(x)$ 在每个小区间上是三次多项式，且 $s,s^{\prime \prime }$ 在内节点处自动连续（因为两侧的表达式都以 $M_i,y_i$ 为公共端点值）。尚未使用的条件是 $s^{\prime }(x)$ 在内节点处的连续性，这将给出确定 $M_i$ 的方程。

要求：

$$s^{\prime }(x_i-0)=s^{\prime }(x_i+0),i=1,2,\cdots ,n-1$$

在 $[x_{i-1},x_i]$ 上代入 $x=x_i$（即 $x_i-x=0$，$x-x_{i-1}=h_i$）：

$$s^{\prime }(x_i-0)=\frac{h_i}{6}{M}_{i-1}+\frac{h_i}{3}{M}_i+\frac{y_i-y_{i-1}}{h_i}$$

在 $[x_i,x_{i+1}]$ 上代入 $x=x_i$（即 $x_{i+1}-x=h_{i+1}$，$x-x_i=0$）：

$$s^{\prime }(x_i+0)=-\frac{h_{i+1}}{3}{M}_i-\frac{h_{i+1}}{6}{M}_{i+1}+\frac{y_{i+1}-y_i}{h_{i+1}}$$

令两者相等：

$$\frac{h_i}{6}{M}_{i-1}+\frac{h_i}{3}{M}_i+\frac{y_i-y_{i-1}}{h_i}=-\frac{h_{i+1}}{3}{M}_i-\frac{h_{i+1}}{6}{M}_{i+1}+\frac{y_{i+1}-y_i}{h_{i+1}}$$

移项整理：

$$\frac{h_i}{6}{M}_{i-1}+\frac{h_i+h_{i+1}}{3}{M}_i+\frac{h_{i+1}}{6}{M}_{i+1}=\frac{y_{i+1}-y_i}{h_{i+1}}-\frac{y_i-y_{i-1}}{h_i}$$

归一化： 两边同乘 $\frac{6}{h_i+h_{i+1}}$，并记

$${\lambda }_i=\frac{h_{i+1}}{h_i+h_{i+1}},{\mu }_i=\frac{h_i}{h_i+h_{i+1}}({\mu }_i+{\lambda }_i=1),d_i=\frac{6}{h_i+h_{i+1}}\left(\frac{y_{i+1}-y_i}{h_{i+1}}-\frac{y_i-y_{i-1}}{h_i}\right)$$

上式化为标准的三对角形式：

$${\mu }_i{M}_{i-1}+2M_i+{\lambda }_i{M}_{i+1}=d_i,i=1,2,\cdots ,n-1$$

这给出 $n-1$ 个方程，但未知数是 $M_i,i=0,1,\cdots ,n$，共 $n+1$ 个，方程数比未知数少 2 个。缺少的两个条件需要由边界条件补充。

边界条件

在两端点 $x_0,x_n$ 处，需要额外给定一阶或二阶导数信息，常见有三种类型：

I 型边界条件（给定一阶导数，即"固支"边界）：

$$s^{\prime }(x_0)=y_0^{\prime },s^{\prime }(x_n)=y_n^{\prime }$$

将 $x=x_0$ 代入区间 $[x_0,x_1]$ 上 $s^{\prime }$ 的表达式：

$$s^{\prime }(x_0+0)=-\frac{h_1}{3}{M}_0-\frac{h_1}{6}{M}_1+\frac{y_1-y_0}{h_1}=y_0^{\prime }$$

同理在 $[x_{n-1},x_n]$ 上代入 $x=x_n$：

$$s^{\prime }(x_n-0)=\frac{h_n}{6}{M}_{n-1}+\frac{h_n}{3}{M}_n+\frac{y_n-y_{n-1}}{h_n}=y_n^{\prime }$$

整理后即可得到关于 $M_0,M_1$ 与 $M_{n-1},M_n$ 的两个附加线性方程。

II 型边界条件（直接给定二阶导数）：

$$M_0=y_0^{\prime \prime },M_n=y_n^{\prime \prime }$$

特别地，当 $y_0^{\prime \prime }=y_n^{\prime \prime }=0$ 时，称为自然边界条件，所得样条称为自然三次样条函数——此时样条曲线在两端的曲率为零，趋于"平直"，这也是物理样条自由端的真实状态（木条两端不受外力弯矩）。

III 型边界条件（一般/周期型边界条件）： 直接给定形如内部方程结构的两个关系式：

$$\{\begin{array}{l}2M_0+{\lambda }_0{M}_1=d_0\\ {\mu }_n{M}_{n-1}+2M_n=d_n\end{array}(0\le {\lambda }_0,{\mu }_n\le 1)$$

这种写法也称为一般边界条件，因为适当选取 ${\lambda }_0,{\mu }_n,d_0,d_n$ 可以把 I 型、II 型边界条件都统一表示成这个形式，从而在理论分析与程序实现中可以只处理一种统一的线性方程组结构。

方程组的求解

以 III 型边界条件为例，把边界方程与内部的连续性方程放在一起，得到关于 $M_0,M_1,\cdots ,M_n$ 的完整线性方程组：

$$\{\begin{array}{ll}2M_0+{\lambda }_0{M}_1=d_0& \\ {\mu }_i{M}_{i-1}+2M_i+{\lambda }_i{M}_{i+1}=d_i,& i=1,2,\cdots ,n-1\\ {\mu }_n{M}_{n-1}+2M_n=d_n& \end{array}$$

写成矩阵形式：

$$\left(\begin{array}{ccccccc}2& {\lambda }_0& & & & & \\ {\mu }_1& 2& {\lambda }_1& & & & \\ & {\mu }_2& 2& {\lambda }_2& & & \\ & & \ddots & \ddots & \ddots & & \\ & & & {\mu }_{n-2}& 2& {\lambda }_{n-2}& \\ & & & & {\mu }_{n-1}& 2& {\lambda }_{n-1}\\ & & & & & {\mu }_n& 2\end{array}\right)\left(\begin{array}{c}{M}_0\\ M_1\\ M_2\\ ⋮\\ M_{n-2}\\ M_{n-1}\\ M_n\end{array}\right)=\left(\begin{array}{c}{d}_0\\ d_1\\ d_2\\ ⋮\\ d_{n-2}\\ d_{n-1}\\ d_n\end{array}\right)$$

可解性： 由 ${\mu }_i+{\lambda }_i=1$（$1\le i\le n-1$）以及对角元恒为 2，可知系数矩阵每一行都满足

$$|2|>{\mu }_i+{\lambda }_i=1$$

即系数矩阵是严格对角占优的三对角矩阵。严格对角占优矩阵必非奇异（这是 Gershgorin 圆盘定理的直接推论），因此方程组的解存在且唯一。

数值解法： 由于系数矩阵是三对角矩阵，可以用高斯消元法的特殊形式——**追赶法（Thomas 算法）**求解，计算量仅为 $O(n)$，远优于一般线性方程组的 $O(n^3)$。

三次样条函数的基表示

正如 Lagrange 插值把插值多项式写成基函数的线性组合一样，三次样条函数也可以写成基样条函数的线性组合，从而避免每次都要重新解三对角方程组。

对于 I 型边界条件，样条函数可写为：

$$s(x)=\sum _{i=0}^n{y}_i{\phi }_i(x)+y_0^{\prime }{\phi }_{n+1}(x)+y_n^{\prime }{\phi }_{n+2}(x)$$

对于 II 型边界条件，样条函数可写为：

$$s(x)=\sum _{i=0}^n{y}_i{\tilde{\phi }}_i(x)+y_0^{\prime \prime }{\tilde{\phi }}_{n+1}(x)+y_n^{\prime \prime }{\tilde{\phi }}_{n+2}(x)$$

其中 ${\phi }_i(x)$（及 ${\tilde{\phi }}_i(x)$）是仅依赖于节点分割 $\{x_i\}$ 而与具体数据 $\{y_i\}$ 无关的基样条函数，它们满足类似 Lagrange 基函数的插值性质。这种表示称为三次样条函数的基样条表示。它的优点是：一旦节点分割给定，基函数只需计算一次；对不同的数据 $\{y_i\}$，样条函数就是这些基函数的线性组合，不必重新求解线性方程组。

---

B 样条曲线

在上面的基样条表示的启发下，可以进一步发展出更一般、局部性更好的样条表示方式——B 样条（B-spline）。

已知 $n+1$ 个控制点：

$${\vec{P}}_0,{\vec{P}}_1,\cdots ,{\vec{P}}_n$$

依次用直线段连接相邻控制点，所得折线称为控制多边形。

给定次数 $p$，以及非降的节点序列（节点向量）：

$$t_0\le t_1\le \cdots \le t_{n+p+1}$$

则 $p$ 次 B 样条曲线定义为：

$$C(t)=\sum _{i=0}^n{N}_{i,p}(t){\vec{P}}_i$$

注意： 这里是控制点 ${\vec{P}}_i$（可以是空间中的向量），而不是型值点上的函数值 $y_i$——B 样条曲线一般并不通过控制点，而是被控制多边形"逼近/引导"，这一点与前面三次样条插值有本质区别。

其中 $N_{i,p}(t)$ 是 $p$ 次 B 样条基函数，由如下 Cox–de Boor 递推公式定义：

$$N_{i,0}(t)=\{\begin{array}{ll}1,& t_i\le t<t_{i+1}\\ 0,& \text{其它}\end{array}$$

$$N_{i,p}(t)=\frac{t-t_i}{t_{i+p}-t_i}{N}_{i,p-1}(t)+\frac{t_{i+p+1}-t}{t_{i+p+1}-t_{i+1}}{N}_{i+1,p-1}(t)(\text{约定 }\frac{0}{0}=0)$$

B 样条基函数具有以下重要性质，使其在曲线设计中远比全局的 Bézier/多项式方法灵活：

• 局部支集： $N_{i,p}(t)$ 只在 $[t_i,t_{i+p+1}]$ 上非零，因此移动一个控制点只会影响曲线的局部形状；
• 非负性： $N_{i,p}(t)\ge 0$；
• 单位分解性： $\sum _{i=0}^n{N}_{i,p}(t)=1$（在有效参数区间内），从而曲线 $C(t)$ 落在控制多边形的凸包内；
• 光滑性： 若节点无重复，$N_{i,p}(t)$ 具有 $C^{p-1}$ 连续性；节点处的重数越高，光滑阶数越低。

三次样条插值可以看作是节点处于插值型值点、且带有插值约束的特殊 B 样条情形，二者思想一脉相承。

---

三次样条函数的性质

以下三个定理说明了三次样条插值（尤其是自然样条）在理论上的最优性，回答了"为什么样条比高次多项式更好"的问题。

定理 1（极小模性质）. 设 $\Delta$ 是 $[a,b]$ 的一个分割，$f(x)\in C^2[a,b]$ 是任一被插值函数（即 $f(x_i)=y_i$）。$s(x)$ 是分割 $\Delta$ 上的自然插值三次样条函数（即 $s^{\prime \prime }(x_0)=0,s^{\prime \prime }(x_n)=0$，且 $s(x_i)=y_i$），则

$${\int }_a^b[s^{\prime \prime }(x){]}^{2}dx\le {\int }_a^b[f^{\prime \prime }(x){]}^{2}dx$$

等号成立当且仅当 $s(x)\equiv f(x)$。

几何意义： 曲线 $f(x)$ 的曲率为

$$k_1=\frac{f^{\prime \prime }(x)}{(1+(f^{\prime }(x){)}^2{)}^{3/2}}$$

当 $|f^{\prime }(x)|\ll 1$（曲线接近水平、斜率很小）时，有近似

$$k_1\approx f^{\prime \prime }(x)$$

类似地 $s(x)$ 的曲率 $k_2\approx s^{\prime \prime }(x)$。因此定理 1 的不等式，在积分（能量）的意义下，表明 $s(x)$ 的曲率不超过 $f(x)$ 的曲率，即 $s(x)$ 的曲率半径不小于 $f(x)$ 的曲率半径——也就是说，在所有通过给定型值点的插值曲线中，自然三次样条是"弯曲能量"（即弯曲势能 $\int (s^{\prime \prime }{)}^{2}dx$，正比于细梁弯曲能量）最小的曲线，最符合物理样条"最省力、最平顺"的自然状态，这正是"样条"一词的由来。

定理 2（最佳逼近性质）. 设 $\Delta$ 是 $[a,b]$ 的一个分割，$f(x)\in C^2[a,b]$。设 $s_f(x)$ 是 $\Delta$ 上带有斜率边界条件（I 型边界条件，且斜率取自 $f$ 本身）的插值三次样条函数，即

$$s_f(x_i)=f(x_i)(i=0,\cdots ,n),s_f^{\prime }(x_0)=f^{\prime }(x_0),s_f^{\prime }(x_n)=f^{\prime }(x_n)$$

设 $s(x)$ 是分割 $\Delta$ 上任意满足 $s(x_i)=f(x_i)$ 的三次样条函数，则

$${\int }_a^b[f^{\prime \prime }(x)-s_f^{\prime \prime }(x){]}^{2}dx\le {\int }_a^b[f^{\prime \prime }(x)-s^{\prime \prime }(x){]}^{2}dx$$

即在所有满足同样插值条件的三次样条函数中，$s_f$（用真实斜率作为边界条件）在二阶导数（曲率）的 $L^2$ 意义下是对 $f$ 的最佳逼近。

定理 3（误差估计）. 设 $\Delta$ 是 $[a,b]$ 的一个分割，$f(x)\in C^4[a,b]$。$s(x)$ 是关于 $f(x)$ 的、带 I 型或 II 型边界条件的插值三次样条函数，则对 $r=0,1,2,3$，有

$$\Vert (f(x)-s(x){)}^{(r)}{\Vert }_{\infty }\le C_r\Vert f^{(4)}{\Vert }_{\infty }{h}^{4-r}$$

其中

$$C_0=\frac{5}{384},C_1=\frac{1}{24},C_2=\frac{3}{8},C_3=\frac{\beta +{\beta }^{-1}}{2}$$

$$h=\underset{i}{\max }{h}_i,\beta =\frac{{\max }_i{h}_i}{{\min }_i{h}_i}$$

（此即经典的 Hall–Meyer 误差估计。注意函数值误差是 $O(h^4)$ 阶的，这正是三次样条相比分段线性插值（$O(h^2)$）精度大幅提升的原因；同时误差与 $f^{(4)}$ 的大小成正比，与高次多项式插值不同，这里的界不依赖节点个数 $n$ 的高次幂，因而不会出现 Runge 现象那样的振荡发散。）

---

参数曲线样条

前面第一至五节讨论的都是"函数型"插值，即假设 $x_0<x_1<\cdots <x_n$，每个 $x$ 对应唯一的 $y$。但很多平面曲线（例如"打圈"的曲线、自相交曲线）无法写成 $y=f(x)$ 的形式。为处理这种一般情形，需要引入参数曲线样条。

问题： 已知平面上的有序点列 $P_i(x_i,y_i),i=0,1,\cdots ,n$，不再假设 $x_0<x_1<\cdots <x_n$（点列可以任意顺序分布，曲线可以自交、打圈）。求平面上的连续曲线 $C$，要求 $C$ 依次通过 $P_i$。

解决思路：引入参数 $t$，把 $x,y$ 都看成 $t$ 的函数分别插值。

常用的做法是累加弦长参数化：用相邻两点间的直线段长度（弦长）近似曲线弧长，累加得到每个型值点对应的参数值。

记相邻两点间的弦长：

$$l_i=\sqrt{(x_i-x_{i-1}{)}^2+(y_i-y_{i-1}{)}^2},i=1,2,\cdots ,n$$

弦长总和：

$$L=\sum _{j=1}^n{l}_j$$

定义各型值点处的参数值：

$$t_0=0,t_i=\sum _{j=1}^i{l}_j,i=1,2,\cdots ,n$$

参数范围为 $[0,L]$，且 $t_0=0<t_1<\cdots <t_n=L$（严格递增，这正是前面函数型插值所需要的条件！）。

（也可以将参数单位化到 $[0,1]$：$t_0=0,t_i=\frac{{\sum }_{j=1}^i{l}_j}{L}$，则 $t_0=0<t_1<\cdots <t_n=1$。）

于是每个型值点变为三元组 $P_i(x_i,y_i,t_i),i=0,1,\cdots ,n$，其中 $t_i$ 严格递增。

分别对 $(t_i,x_i)$ 与 $(t_i,y_i)$ 做前述的三次样条插值，得到两条以 $t$ 为自变量的样条函数

$$s_x(t),s_y(t)$$

满足

$$s_x(t_i)=x_i,s_y(t_i)=y_i,i=0,1,\cdots ,n$$

于是参数样条曲线

$$(s_x(t),s_y(t)),t\in [0,L]$$

就是依次通过所有 $P_i$ 的连续曲线，且不再要求 $x$ 坐标单调，从而可以描述打圈、自交等复杂曲线。

空间曲线（$P_i(x_i,y_i,z_i)$）的插值问题可以完全类似地处理：对 $(t_i,x_i)$、$(t_i,y_i)$、$(t_i,z_i)$ 分别做三次样条插值。但由于弦长参数化在三维空间中对曲率、挠率变化剧烈的曲线近似误差更大，所以处理效果一般不如平面曲线理想。

---

Hermite 插值

三次样条插值中，$M_i=s^{\prime \prime }(x_i)$ 是通过整体线性方程组间接求出的（依赖所有型值点），这带来了"全局性"——修改一个型值点会影响整条曲线。Hermite 插值则换一个角度：直接在每个型值点处指定一阶导数值，从而可以做到分段独立构造（每一段只依赖该段两端的信息，具有更好的局部性）。

问题： 已知平面上的有序点列 $P_i(x_i,y_i),i=0,1,\cdots ,n$，假设

$$x_0<x_1<\cdots <x_{n-1}<x_n$$

并已知在 $P_i(x_i,y_i)$ 处的一阶导数值 $\{m_i{\}}_{i=0}^n$。求平面上的连续曲线 $C:y=f(x)$，满足

$$f(x_i)=y_i,f^{\prime }(x_i)=m_i,i=0,1,\cdots ,n$$

构造思路： 在每个小区间 $[x_{i-1},x_i]$ 上，共有 4 个条件——$y_{i-1},y_i,m_{i-1},m_i$——恰好唯一确定一个三次多项式（与三次样条插值中用 $y_{i-1},y_i,M_{i-1},M_i$ 确定三次多项式的思路完全类似，只是这里用一阶导数代替二阶导数作为附加条件）。

利用三次 Hermite 基函数直接构造，可以验证：在 $[x_{i-1},x_i]$ 上，

$$\begin{array}{rl}f(x)& =m_{i-1}\frac{(x_i-x{)}^2(x-x_{i-1})}{h_i^2}-m_i\frac{(x-x_{i-1}{)}^2(x_i-x)}{h_i^2}\\ & +y_{i-1}\frac{(x_i-x{)}^2[2(x-x_{i-1})+h_i]}{h_i^3}+y_i\frac{(x-x_{i-1}{)}^2[2(x_i-x)+h_i]}{h_i^3}\end{array}$$

其中 $h_i=x_i-x_{i-1}$，$x\in [x_{i-1},x_i]$。

验证要点：

• 代入 $x=x_{i-1}$：第一、二项含因子 $(x-x_{i-1})$ 或 $(x_i-x{)}^2(x-x_{i-1})$ 中的 $(x-x_{i-1})$ 均为零，只剩第三项 $y_{i-1}\cdot \frac{h_i^2[0+h_i]}{h_i^3}=y_{i-1}$，第四项含 $(x-x_{i-1}{)}^2=0$；故 $f(x_{i-1})=y_{i-1}$，同理 $f(x_i)=y_i$。
• 对 $f(x)$ 求导并代入端点，可验证 $f^{\prime }(x_{i-1})=m_{i-1}$，$f^{\prime }(x_i)=m_i$（读者可自行展开验证）。

由于每一段的三次多项式只依赖该段端点处的 $(y,m)$ 值，Hermite 插值天然是 $C^1$ 连续的分段构造，比三次样条更"局部"，但代价是需要额外提供导数信息 $\{m_i\}$（工程中常用相邻弦的斜率或有限差分近似给出）。

---

曲线的最小二乘逼近

前面讨论的都是插值问题：要求曲线严格通过所有数据点。但当数据点很多、且带有测量噪声时，严格插值反而会拟合噪声、产生振荡，这时更合理的做法是逼近（拟合）：只要求曲线在某种整体误差度量下"最接近"数据点，不必精确通过每一个点。

多项式最小二乘逼近

设逼近曲线为 $N$ 次多项式（$N\ll M$，见下）：

$$y=a_0+a_{1}x+\cdots +a_{N-1}{x}^{N-1}+a_N{x}^N$$

取平面上 $M$ 个数据点 $P_i(x_i,y_i),i=1,\cdots ,M$，其中数据点个数 $M$ 远大于多项式的自由度 $N+1$（即 $N\ll M$），因此一般无法要求曲线精确通过所有点。

定义残差平方和：

$$F(a_0,a_1,\cdots ,a_N)=\sum _{i=1}^M(y(x_i)-y_i{)}^2$$

我们要求

$$\underset{a_0,\cdots ,a_N}{\min }F(a_0,a_1,\cdots ,a_N)$$

最小二乘法： 令 $F$ 关于每个 $a_j$ 的偏导数为零（这是 $F$ 取极小值的必要条件，且 $F$ 是关于 $\{a_j\}$ 的凸二次函数，故驻点即为全局最小点）：

$$\frac{\partial F}{\partial a_j}=2\sum _{i=1}^M(y(x_i)-y_i)x_i^j=0,j=0,1,\cdots ,N$$

将 $y(x_i)={\sum }_{k=0}^N{a}_k{x}_i^k$ 代入并展开，得到

$$a_0\sum _{i=1}^M{x}_i^j+a_1\sum _{i=1}^M{x}_i^{j+1}+\cdots +a_N\sum _{i=1}^M{x}_i^{j+N}=\sum _{i=1}^M{x}_i^j{y}_i,j=0,1,\cdots ,N$$

写成矩阵形式：

$$\left(\begin{array}{cccc}M& {\sum }_{i=1}^M{x}_i& \cdots & {\sum }_{i=1}^M{x}_i^N\\ {\sum }_{i=1}^M{x}_i& {\sum }_{i=1}^M{x}_i^2& \cdots & {\sum }_{i=1}^M{x}_i^{N+1}\\ ⋮& ⋮& \ddots & ⋮\\ {\sum }_{i=1}^M{x}_i^N& {\sum }_{i=1}^M{x}_i^{N+1}& \cdots & {\sum }_{i=1}^M{x}_i^{2N}\end{array}\right)\left(\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_N\end{array}\right)=\left(\begin{array}{c}{\sum }_{i=1}^M{y}_i\\ {\sum }_{i=1}^M{x}_i{y}_i\\ ⋮\\ {\sum }_{i=1}^M{x}_i^N{y}_i\end{array}\right)$$

这一线性方程组称为法方程（正规方程，normal equations）。解得 $a_i,i=0,1,\cdots ,N$，即得所求的最小二乘逼近多项式。

一般基函数系下的最小二乘逼近

多项式基 $\{1,x,x^2,\cdots ,x^N\}$ 只是众多可能的基函数系之一，一般情形下可以用任意一组线性无关函数系来逼近。

设 $\{{\phi }_j(x){\}}_{j=0}^N$ 为 $[a,b]$ 上一组线性无关的函数，要求用它们的线性组合

$$y=a_0{\phi }_0(x)+a_1{\phi }_1(x)+\cdots +a_N{\phi }_N(x)=\sum _{j=0}^N{a}_j{\phi }_j(x)$$

来逼近平面上给定的数据点 $P_i(x_i,y_i),i=1,\cdots ,M$（$x_i\in [a,b]$，$N\ll M$）。

同样定义残差平方和：

$$F(a_0,a_1,\cdots ,a_N)=\sum _{i=1}^M(y(x_i)-y_i{)}^2=\sum _{i=1}^M{\left(\sum _{j=0}^N{a}_j{\phi }_j(x_i)-y_i\right)}^2$$

对 $a_k$ 求偏导数并令其为零：

$$\frac{\partial F}{\partial a_k}=2\sum _{i=1}^M\left(\sum _{j=0}^N{a}_j{\phi }_j(x_i)-y_i\right){\phi }_k(x_i)=0$$

即

$$\sum _{i=1}^M\sum _{j=0}^N{a}_j{\phi }_j(x_i){\phi }_k(x_i)=\sum _{i=1}^M{y}_i{\phi }_k(x_i),k=0,1,\cdots ,N$$

记

$$S_{j,k}=\sum _{i=1}^M{\phi }_j(x_i){\phi }_k(x_i),T_k=\sum _{i=1}^M{y}_i{\phi }_k(x_i),k=0,1,\cdots ,N$$

则法方程可写为

$$\sum _{j=0}^N{a}_j{S}_{j,k}=T_k,k=0,1,\cdots ,N$$

矩阵形式：

$$\left(\begin{array}{cccc}{S}_{0,0}& S_{1,0}& \cdots & S_{N,0}\\ S_{0,1}& S_{1,1}& \cdots & S_{N,1}\\ ⋮& ⋮& \ddots & ⋮\\ S_{0,N}& S_{1,N}& \cdots & S_{N,N}\end{array}\right)\left(\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_N\end{array}\right)=\left(\begin{array}{c}{T}_0\\ T_1\\ ⋮\\ T_N\end{array}\right)$$

由 $S_{j,k}=S_{k,j}$（求和顺序可交换）可知系数矩阵是对称矩阵。

特别地，如果 $\{{\phi }_i(x){\}}_{i=0}^N$ 关于离散点集 $\{x_i\}$ 构成正交系（即 ${\sum }_i{\phi }_j(x_i){\phi }_k(x_i)=0,j\ne k$），则 $S_{j,k}=0(j\ne k)$，系数矩阵退化为对角阵，此时

$$a_k=\frac{T_k}{S_{k,k}}=\frac{{\sum }_{i=1}^M{y}_i{\phi }_k(x_i)}{{\sum }_{i=1}^M{\phi }_k(x_i{)}^2}$$

各系数可以独立、直接算出，无需求解线性方程组——这正是选用正交多项式（如 Chebyshev 多项式、Legendre 多项式）作基函数进行最小二乘拟合的根本原因。

加权最小二乘逼近

若不同数据点的可信度不同（例如测量精度不同），可以引入权因子 ${\omega }_i\ge 0,i=1,\cdots ,M$，定义加权残差平方和：

$$F(a_0,a_1,\cdots ,a_N)=\sum _{i=1}^M{\omega }_i{\left(\sum _{j=0}^N{a}_j{\phi }_j(x_i)-y_i\right)}^2$$

权重越大的点，在拟合中被赋予越高的重要性（迫使曲线更贴近该点）。

同样对 $a_k$ 求偏导数并令其为零：

$$\sum _{i=1}^M\sum _{j=0}^N{\omega }_i{a}_j{\phi }_j(x_i){\phi }_k(x_i)=\sum _{i=1}^M{\omega }_i{y}_i{\phi }_k(x_i),k=0,1,\cdots ,N$$

令

$$S_{j,k}=\sum _{i=1}^M{\omega }_i{\phi }_j(x_i){\phi }_k(x_i),T_k=\sum _{i=1}^M{\omega }_i{y}_i{\phi }_k(x_i),k=0,1,\cdots ,N$$

则法方程与矩阵形式和无权重情形完全相同：

$$\sum _{j=0}^N{a}_j{S}_{j,k}=T_k,k=0,1,\cdots ,N$$

取 ${\omega }_i\equiv 1$ 即退化为一般基函数系的普通最小二乘情形，两者是统一的框架。

---

曲面插值问题

将曲线插值推广到二元函数（曲面）上，最直接的方式是利用直积（张量积）结构。

双三次样条函数 由一元三次样条函数的直积构成：本质上仍是（两次）一元曲线问题的叠加，只是需要特别注意型值点必须分布在矩形网格点上（而不能是散乱点）。

设平面上的网格分割：

$$x_0<x_1<\cdots <x_n,y_0<y_1<\cdots <y_m$$

已知网格节点处的型值点（函数值）：

$$(x_i,y_j,z_{ij}),i=0,\cdots ,n,j=0,\cdots ,m$$

要求在矩形区域 $[x_0,x_n]\times [y_0,y_m]$ 上构造曲面 $z=z(x,y)$，满足插值条件

$$z(x_i,y_j)=z_{ij}$$

构造思路（直积法）： 先沿 $x$ 方向对每一行 $y=y_j$ 固定时的数据 $\{(x_i,z_{ij}){\}}_i$ 做一元三次样条插值，得到一族依赖参数 $y_j$ 的样条曲线；再沿 $y$ 方向对样条系数（如 $M_i$ 或基函数系数）做插值，即可将一元三次样条的基表示

$$s(x)=\sum _i{y}_i{\phi }_i(x)$$

自然地推广为二元的张量积形式：

$$z(x,y)=\sum _{i=0}^n\sum _{j=0}^m{z}_{ij}{\phi }_i(x){\psi }_j(y)$$

其中 ${\phi }_i(x)$、${\psi }_j(y)$ 分别是对应 $x$、$y$ 方向节点分割的一元三次样条基函数。这样构造出的曲面在每个网格小矩形 $[x_{i-1},x_i]\times [y_{j-1},y_j]$ 上是关于 $x,y$ 的双三次多项式（即形如 ${\sum }_{p,q=0}^3{c}_{pq}{x}^p{y}^q$），并在网格线上保持相应阶数的连续性。

这种方法结构简单、计算量小（可以逐行逐列地用一维追赶法求解），但要求数据必须分布在矩形网格上，这是它的局限所在。

---

散乱数据插值方法

实际应用中（如地形测绘、气象观测、有限元后处理等），数据点往往不在规则网格上，而是任意散布在平面区域内（散乱数据，scattered data）。这时直积法不再适用，需要专门的散乱数据插值方法。

问题： 给定平面上的区域 $D$，给定 $D$ 中的 $n$ 个点

$$(x_i,y_i)\in D,i=1,2,\cdots ,n$$

已知各点处的函数值 $f_i=f_i(x_i,y_i),i=1,2,\cdots ,n$。要求构造定义在 $D$ 上的函数 $F(x,y)$，$(x,y)\in D$，满足插值条件

$$F(x_i,y_i)=f_i(x_i,y_i),i=1,\cdots ,n$$

Shepard 方法（反距离加权法）

基本思想：某一点 $(x,y)$ 处的插值应主要由离它最近的数据点决定，距离越远的数据点权重应越小。

对任意 $(x,y)\in D$，记该点到各数据点的欧氏距离：

$${\rho }_i(x,y)=\sqrt{(x-x_i{)}^2+(y-y_i{)}^2},i=1,2,\cdots ,n$$

定义 $D$ 上的插值函数：

$$F(x,y)=\{\begin{array}{ll}{f}_i,& {\rho }_i(x,y)=0\\ \frac{1}{\sum _{j=1}^n\frac{1}{{\rho }_j^2(x,y)}}\sum _{i=1}^n\frac{f_i}{{\rho }_i^2(x,y)}=\sum _{i=1}^n{\omega }_i(x,y)f_i,& \text{其它}\end{array}$$

其中权函数为

$${\omega }_i(x,y)=\frac{\frac{1}{{\rho }_i^2(x,y)}}{\sum _{j=1}^n\frac{1}{{\rho }_j^2(x,y)}}$$

基本性质：

• 单位分解性： 由权函数定义直接可得

  $$\sum _{i=1}^n{\omega }_i(x,y)=1$$

  因而 $F(x,y)$ 是所有 $f_i$ 的凸组合，取值落在 ${\min }_i{f}_i$ 与 ${\max }_i{f}_i$ 之间。
• 插值性： 当 $(x,y)\to (x_k,y_k)$ 时，${\rho }_k\to 0$，权重 ${\omega }_k\to 1$、其余 ${\omega }_i\to 0$（因为 $1/{\rho }_k^2\to \infty$ 而其余项有界），配合 ${\rho }_i=0$ 处的特殊定义 $F=f_i$，可验证 $F$ 在每个数据点处连续且 $F(x_i,y_i)=f_i$。

Shepard 方法实现简单、计算量小（对每个新点只需 $O(n)$ 次运算），但其插值曲面在数据点附近容易出现"平顶"或"凹坑"现象（因为在数据点处 $F$ 的梯度通常为零），光滑性不如样条方法。

推广到向量场情形

Shepard 方法可以直接推广到向量值函数（例如有限元中的位移场插值）。

给定区域 $D$，给定点集 ${\vec{x}}_i\in D,i=1,2,\cdots ,n$，已知各点处的向量值 ${\vec{v}}_i={\vec{v}}_i({\vec{x}}_i)$。要求构造定义在 $D$ 上的向量函数 $\vec{V}(\vec{x})$，$\vec{x}\in D$，满足

$$\vec{V}({\vec{x}}_i)={\vec{v}}_i({\vec{x}}_i)$$

对任意 $\vec{x}\in D$，记

$${\rho }_i(\vec{x})=\Vert \vec{x}-{\vec{x}}_i\Vert ,i=1,2,\cdots ,n$$

定义 $D$ 上的向量值插值函数：

$$\vec{V}(\vec{x})=\{\begin{array}{ll}{\vec{v}}_i,& {\rho }_i(\vec{x})=0\\ \frac{1}{\sum _{j=1}^n\frac{1}{{\rho }_j^2(\vec{x})}}\sum _{i=1}^n\frac{{\vec{v}}_i}{{\rho }_i^2(\vec{x})}=\sum _{i=1}^n{\omega }_i(\vec{x}){\vec{v}}_i,& \text{其它}\end{array}$$

其中权函数 ${\omega }_i(\vec{x})$ 的定义与标量情形完全相同。由于各分量权重相同，这相当于对向量的每个分量分别做一次标量 Shepard 插值。

Hardy 插值方法（多二次曲面法，MQI）

Shepard 方法的权函数是分片有理函数，整体光滑性较差；Hardy 提出的多二次曲面法（Multiquadric Interpolation, MQI）通过引入光滑的径向基函数来改善这一点，是径向基函数（RBF）插值方法的早期代表。

设插值函数取如下径向基函数的线性组合形式：

$$\vec{u}(\vec{x})=\sum _{j=1}^n{\vec{\alpha }}_j{\varphi }_j(\vec{x})$$

其中基函数为

$${\varphi }_j(\vec{x})=\sqrt{{\Vert \vec{x}-{\vec{x}}_j\Vert }^2+R^2}$$

$R$ 是可调节的正常数（形状参数），$R$ 越大基函数越"平缓"。由于 ${\varphi }_j$ 是 $\Vert \vec{x}-{\vec{x}}_j\Vert$ 的光滑函数（不同于 Shepard 权重在数据点处不可导），Hardy 插值得到的曲面整体光滑性更好。

求解待定系数 ${\vec{\alpha }}_j$： 由插值条件 $\vec{u}({\vec{x}}_i)={\vec{v}}_i,i=1,\cdots ,n$，代入得到关于 ${\vec{\alpha }}_1,\cdots ,{\vec{\alpha }}_n$ 的线性方程组：

$$\{\begin{array}{l}{\varphi }_1({\vec{x}}_1){\vec{\alpha }}_1+{\varphi }_2({\vec{x}}_1){\vec{\alpha }}_2+\cdots +{\varphi }_n({\vec{x}}_1){\vec{\alpha }}_n={\vec{v}}_1\\ {\varphi }_1({\vec{x}}_2){\vec{\alpha }}_1+{\varphi }_2({\vec{x}}_2){\vec{\alpha }}_2+\cdots +{\varphi }_n({\vec{x}}_2){\vec{\alpha }}_n={\vec{v}}_2\\ ⋮\\ {\varphi }_1({\vec{x}}_n){\vec{\alpha }}_1+{\varphi }_2({\vec{x}}_n){\vec{\alpha }}_2+\cdots +{\varphi }_n({\vec{x}}_n){\vec{\alpha }}_n={\vec{v}}_n\end{array}$$

系数矩阵的对称性： 容易验证

$${\varphi }_j({\vec{x}}_i)=\sqrt{{\Vert {\vec{x}}_i-{\vec{x}}_j\Vert }^2+R^2}=\sqrt{{\Vert {\vec{x}}_j-{\vec{x}}_i\Vert }^2+R^2}={\varphi }_i({\vec{x}}_j)$$

（因为欧氏距离本身对称：$\Vert {\vec{x}}_i-{\vec{x}}_j\Vert =\Vert {\vec{x}}_j-{\vec{x}}_i\Vert$），所以插值方程组的系数矩阵是对称矩阵。

事实上可以进一步证明该系数矩阵（Hardy 矩阵）是正定的（当 ${\vec{x}}_i$ 两两不同时），因此方程组存在唯一解 ${\vec{\alpha }}_1,{\vec{\alpha }}_2,\cdots ,{\vec{\alpha }}_n$，代回即得到插值位移向量函数 $\vec{u}(\vec{x})$。

与 Shepard 方法相比，Hardy 方法（及一般的径向基函数方法）需要求解一个 $n\times n$ 的稠密线性方程组（计算量更大，当 $n$ 很大时通常需要局部化或快速多极子等加速技术），但换来的是更好的整体光滑性和逼近精度，因此在散乱数据曲面重建、变形场插值等工程问题中被广泛使用。

---

内容小结

方法	适用场景	光滑性	是否插值	主要求解手段
Lagrange 多项式插值	少量点、理论分析	$C^{\infty }$（全局）	是	直接公式，但数值不稳定（Runge 现象）
三次样条插值	有序型值点	$C^2$（分段三次）	是	三对角方程组（追赶法）
B 样条曲线	曲线/曲面设计，控制点驱动	可调（$C^{p-1}$）	一般不插值，逼近控制多边形	递推基函数
Hermite 插值	已知导数信息、需要局部性	$C^1$（分段三次）	是	逐段直接构造，无需解方程组
参数样条	打圈/自交曲线、空间曲线	同三次样条	是	弦长参数化 + 三次样条
多项式/一般基最小二乘	大量带噪声数据	由基函数决定	否（整体拟合）	法方程
双三次样条	网格型曲面数据	$C^2\times C^2$	是	一元样条的直积
Shepard 法	散乱数据，快速插值	数据点处不光滑	是	反距离加权，无需解方程组
Hardy（MQI）法	散乱数据，要求光滑曲面	全局光滑	是	稠密对称正定线性方程组

说明： 原始笔记中定理 3 常数 $C_0$ 处笔误为 $\frac{5}{385}$，按照三次样条误差估计的经典结果（Hall–Meyer 定理）已订正为 $\frac{5}{384}$。

完整讲义