插值算法
请注意，本文最后更新于2022.2.6，其中一些理解可能已被笔者证伪或废弃。
定义：给定一组离散点列，要求一条曲线将点依次连接，称之为插值[1]。
利用已知的点建立合适的插值函数 \(f(x)\) ,未知点 \(x_i\) 由插值函数 \(f(x)\) 可以求得 \((x_i,f(x_i))\) 近似代替未知点[3]。
作用：利用插值曲线可以对数据进行填充——用少量模拟产生一些靠谱的新数据。
分段插值
分段线性插值
这是最简单也最基础的插值方法，折现统计图的连接方式就是用的分段线性插值。使用函数表时一般直接用该方法。
分段二次插值
取结点 \(x_i\) 以及其左右的三个节点 \(x_{i-1},x_{i},x_{i+1}\) 进行在区间 \([x_{i-1},x_{i+1}]\) 的二次函数插值，即
\[f(x) \approx L_{2}(x)=\sum_{k=i-1}^{i+1}\left[y_{k} \prod_{j=i-1 \atop j \neq k}^{i+1} \frac{\left(x-x_{j}\right)}{\left(x_{k}-x_{i}\right)}\right]\]
其意义是用分段抛物线代替 \(y=f(x)\) [2]。这个方法能保证函数的连续型与一定看起来的平滑性，但是不保证函数光滑，即不保证函数可导。
分段二维插值可以很好的避免多项式插值带来的龙格效应，同时简单易于理解又能看起来平滑，是很常用的方法。
多项式插值
多项式插值的本质是对于 \(n+1\) 个互不相同的节点 \((x_i,y_i) (i=0,1,2,\dots,n)\) 求得一个唯一的多项式： \[L_{n}(x)=a_{0}+a_{1} x+a_{2} x^{2}+\ldots+a_{n} x^{n}\] 使得 \(L_n(x_i)=y_i\) ,即 \[\left[\begin{array}{cccc} 1 & x_{0} & \cdots & x_{0}^{n} \\ 1 & x_{1} & \cdots & x_{1}^{n} \\ \vdots & \vdots & \cdots & \vdots \\ 1 & x_{n} & \cdots & x_{n}^{n} \end{array}\right] \left[\begin{array}{c} a_{0} \\ a_{1} \\ \vdots \\ a_{n} \end{array}\right] =\left[\begin{array}{c} y_{0} \\ y_{1} \\ \vdots \\ y_{n} \end{array}\right]\]
对于上述范德蒙行列式，易得，\([a_1,a_2,\dotsb]\) 有且仅有一组解。
拉格朗日插值法
上多项式组易解得
\[L(x)=y_1\frac{(x-x_2)(x-x_3)...(x-x_n)}{(x_1-x_2)(x_1-x_3)...(x_1-x_n)}+y_2\frac{(x-x_1)(x-x_3)...(x-x_n)}{(x_2-x_1)(x_2-x_3)...(x_2-x_n)}...+y_{n}\frac{(x-x_1)(x-x_2)...(x-x_{n-1})}{(x_{n}-x_1)(x_{n}-x_2)...(x_{n}-x_{n-1})}\] 即 \[L(x)=\sum_{i=1}^{n}y_i \prod_{j=1,j \neq i}^n \frac{x-x_j}{x_i-x_j} \] 为拉格朗日插值法。
注：拉格朗日插值公式在理论分析理解上很容易理解，但是若插值节点发生改变时，插值公式随之就要重新计算生成，在实际计算中会占用大量的计算量。牛顿法的出现正是克服这个问题[3]。
牛顿插值法
\[\begin{aligned} f(x)=& f\left(x_{0}\right)+f\left[x_{0}, x_{1}\right]\left(x-x_{0}\right) \\ &+f\left[x_{0}, x_{1}, x_{2}\right]\left(x-x_{0}\right)\left(x-x_{1}\right)+\cdots \\ &+f\left[x_{0}, x_{1}, \cdots, x_{n-2}, x_{n-1}\right]\left(x-x_{0}\right)\left(x-x_{1}\right) \cdots\left(x-x_{n-3}\right)\left(x-x_{n-2}\right) \\ &+f\left[x_{0}, x_{1}, \cdots, x_{n-1}, x_{n}\right]\left(x-x_{0}\right)\left(x-x_{1}\right) \cdots\left(x-x_{n-2}\right)\left(x-x_{n-1}\right) \end{aligned}\]
商差
- 一阶差商
\[f\left[x_{i}, x_{j}\right]=\frac{f\left(x_{i}\right)-f\left(x_{j}\right)}{x_{i}-x_{j}} \quad\left(i \neq j, x_{i} \neq x_{j}\right)\]
- 二阶差商
\[f\left[x_{i}, x_{j}, x_{k}\right]=\frac{f\left[x_{i}, x_{j}\right]-f\left[x_{j}, x_{k}\right]}{x_{i}-x_{k}} \quad(i \neq k)\]
- n阶差商
\[f\left[x_{0}, x_{1}, \Lambda x_{n}\right]=\frac{f\left[x_{0}, x_{1}, \Lambda, x_{n-1}\right]-f\left[x_{1}, x_{2}, \Lambda x_{n}\right]}{x_{0}-x_{n}}\]
由以上定义我们的到差商表如下：
\[\begin{array}{cccrcc} x_{i} & y_{i} & \text { 一阶差商 } & \text { 二阶差商 } & \cdots \cdots & n \text { 阶差商 } \\ & & & & \\ x_{0} & f\left(x_{0}\right) & & & \\ x_{1} & f\left(x_{1}\right) & f\left[x_{0}, x_{1}\right] & & \\ x_{2} & f\left(x_{2}\right) & f\left[x_{1}, x_{2}\right] & f\left[x_{0}, x_{1}, x_{2}\right] & \\ & \cdots \cdots & \cdots \cdots & \cdots \cdots & \\ x_{n-1} & f\left(x_{n-1}\right) & \cdots \cdots & \cdots \cdots \\ x_{n} & f\left(x_{n}\right) & f\left[x_{n-1}, x_{n}\right] & f\left[x_{n-2}, x_{n-1}, x_{n}\right] &\dots &f\left[x_{n} \ldots, x_{n}\right] \\ \end{array}\] 我们可以得到下面公式： \[\begin{array}{l} f(x)=f\left(x_{\theta}\right)+f\left[x, x_{\theta}\right]\left(x-x_{\theta}\right),\\ f\left[x, x_{0}\right]=f\left[x_{0}, x_{1}\right]+f\left[x, x_{0}, x_{1}\right]\left(x-x_{1}\right),\\ f\left[x, x_{0}, x_{1}\right]=f\left[x_{0}, x_{1}, x_{2}\right]+f\left[x, x_{0}, x_{1}, x_{2}\right]\left(x-x_{2}\right),\\ \cdots \\ f\left[x, x_{\theta}, \cdots, x_{m-1}\right]=f\left[x_{\theta}, x_{1}, \cdots, x_{n}\right]+J\left[x, x_{\theta}, \cdots, x_{n}\right]\left(x-x_{n}\right) \end{array}\] 可得到牛顿插值法。
拉格朗日插值与牛顿插值
牛顿插值法和拉格朗日插值法两者都是多项式插值法，本质上，两者的结果相同，只不过表示的形式不同。牛顿插 值法的计算过程具有继承性有易于变动的特点。
龙格现象
当次数渐高，很容易在多项式插值的边界处产生巨大波动，称为龙格现象。
这里是龙格现象产生原因
埃尔米特(Hermite)插值
不但要求在节点处上的函数值相等，而且还要求对应的导数值也相等，甚至要求高阶导数也相等，满足这种要求的插值多项式就是“Hermite插值多项式”[4]。
直接使用Hermite插值得到的多项式次数较高，也存在着“龙格现象（Runge phenomenon）”。因此，在实际应用中，往往使用分段三次Hermite插值多项式（PCHIP），来提高“模拟数据的准确性”[4]。
原理
假定已知函数 \(f(x)\) 在揷值区间 \([p, q]\) 上的 \(n+1\) 个互不相同的节点 \(x_{i}(i=0,1, \ldots, n)\) 处满足 \(f\left(x_{i}\right)=f_{i}\) 及 \(f^{\prime}\left(x_{i}\right)=f_{i}^{\prime}(i=0,1,2, \ldots, n)\) , 如果函数 \(G(x)\) 的存在满足下列条件: 1. \(G(x)\) 在每个小区间上的多项式次数为 3 ; 2. \(G(x) \in C^{1}[a, b]\) ; 3. \(G\left(x_{i}\right)=f\left(x_{i}\right), \quad G^{\prime}\left(x_{i}\right)=f^{\prime}\left(x_{i}\right), \quad i=(0,1, \ldots, n)\) 就称 \(G(x)\) 是 \(f(x)\) 在 \(n+1\) 个节点 \(x_{i}\) 上的分段三次埃尔米特插值多项式。 \[ \begin{aligned} G(x)=& h_{k} y_{k}(x)+h_{k+1} y_{k+1}(x)+H_{k}(x) y_{k}^{\prime}+H_{k+1}(x) y_{k+1}^{\prime} \\ =&\left(1+2 \frac{x-x_{k}}{x_{k+1}-x_{k}}\right)\left(\frac{x-x_{k+1}}{x_{k}-x_{k+1}}\right)^{2} y_{k}+\left(1+2 \frac{x-x_{k+1}}{x_{k}-x_{k+1}}\right)\left(\frac{x-x_{k}}{x_{k+1}-x_{k}}\right)^{2} y_{k+1} \\ &+\left(x-x_{k}\right)\left(\frac{x-x_{k+1}}{x_{k}-x_{k+1}}\right)^{2} y_{k}^{\prime}+\left(x-x_{k+1}\right)\left(\frac{x-x_{k}}{x_{k+1}-x_{k}}\right)^{2} y_{k+1}^{\prime} \end{aligned} \]
推广
Hermite 揷值多项式为: \[ H(x)=\sum_{i=0}^{n} h_{i}\left[\left(x_{i}-x\right)\left(2 a_{i} y_{i}-y_{i}^{\prime}\right)+y_{i}\right]\] 其中: \[h_{i}=\prod_{j=0 \atop j \neq i}^{n}\left(\frac{x-x_{j}}{x_{i}-x_{j}}\right)^{2}, \quad a_{i}=\sum_{j=0 \atop j \neq i}^{n} \frac{1}{x_{i}-x_{j}}\]
代码
手写函数
设n个节点的数据以数组
x0（已知点的横坐标），y0（函数值），y1（导数值）输入（注意
Matlat 的数组下标从 \(1\) 开始），\(m\) 个插值点以数组 x
输入，输出数组 y 为 \(m\)
个插值。编
function y=hermite(x0,y0,y1,x);
n=length(x0);m=length(x);
for k=1:m
yy=0.0;
for i=1:n
h=1.0;
a=0.0;
for j=1:n
if j~=i
h=h*((x(k)-x0(j))/(x0(i)-x0(j)))^2;
a=1/(x0(i)-x0(j))+a;
end
end
yy=yy+h*((x0(i)-x(k))*(2*a*y0(i)-y1(i))+y0(i));
end
y(k)=yy;
endx = -3:3;
y = [-1 -1 -1 0 1 1 1];
xq1 = -3:.01:3;
p = pchip(x,y,xq1);
s = spline(x,y,xq1);
m = makima(x,y,xq1);
plot(x,y,'o',xq1,p,'-',xq1,s,'-.',xq1,m,'--')
legend('Sample Points','pchip','spline','makima','Location','SouthEast')
| 
---|---
样条插值
所谓样条（Spline）本来是工程设计中使用的一种绘图工具，它是富有弹性的细木条或细金属条。绘图员利用它把一些已知点连接成一条光滑曲线（称为样条曲线），并使连接点处有连续的曲率[4]。
数学上将具有一定光滑性的分段多项式称为样条函数。
对于 \(\Delta : \quad a=x_{0}<x_{1}<\cdots<x_{n-1}<x_{n}=b\) 满足：
- 在每个小区间 \(\left[x_{i}, x_{i-1}\right](i=0,1, \cdots, n-1)\) 上 \(s(x)\) 是 \(k\) 次多项式。
- \(s(x)\) 在 \([a, b]\) 上具有 \(k-1\) 阶连续导数。即 \(S_{i}^{(j)}(x_{i+1})=S_{i+1}^{(j)}(x_{i+1})\) \((j=0,1,\dots,k-1)\)
一般形式为：
\[s_{k}(x)=\sum_{i=0}^{k} \frac{\alpha_{i} x^{i}}{i !}+\sum_{j=1}^{n-1} \frac{\beta_{j}}{k !}\left(x-x_{j}\right)_{+}^{k}\]
在实际中最常用的是二次样条函数和三次样条函数。
三次样条生成的曲线更加光滑，相较Hermite更适合曲线(三角函数)。
三次样条插值代码
x = ‐pi:pi;
y = sin(x);
new_x = ‐pi:0.1:pi;
p1 = pchip(x,y,new_x); %分段三次埃尔米特插值
p2 = spline(x,y,new_x); %三次样条插值
pp = csape(x0,y0_ext,conds); %另一个三次样条插值？
plot(x,y,'o',new_x,p1,'r‐',new_x,p2,'b‐')
legend('样本点','三次埃尔米特插值','三次样条插值','Location','SouthEast') %标注显示在东南方向