插值算法

请注意,文章中部分理解可能已被笔者废弃,本文最后更新于:8 个月前

请注意,本文最后更新于2022.2.6,其中一些理解可能已被笔者证伪或废弃。

定义:给定一组离散点列,要求一条曲线将点依次连接,称之为插值[1]

利用已知的点建立合适的插值函数 f(x)f(x) ,未知点 xix_i 由插值函数 f(x)f(x) 可以求得 (xi,f(xi))(x_i,f(x_i)) 近似代替未知点[3]

作用:利用插值曲线可以对数据进行填充——用少量模拟产生一些靠谱的新数据。

分段插值

分段线性插值

这是最简单也最基础的插值方法,折现统计图的连接方式就是用的分段线性插值。使用函数表时一般直接用该方法。

某折线统计图

分段二次插值

取结点 xix_i 以及其左右的三个节点 xi1,xi,xi+1x_{i-1},x_{i},x_{i+1} 进行在区间 [xi1,xi+1][x_{i-1},x_{i+1}] 的二次函数插值,即

f(x)L2(x)=k=i1i+1[ykj=i1jki+1(xxj)(xkxi)]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)y=f(x) [2]。这个方法能保证函数的连续型与一定看起来的平滑性,但是不保证函数光滑,即不保证函数可导。

分段二维插值可以很好的避免多项式插值带来的龙格效应,同时简单易于理解又能看起来平滑,是很常用的方法。

多项式插值

多项式插值的本质是对于 n+1n+1 个互不相同的节点 (xi,yi)(i=0,1,2,,n)(x_i,y_i) (i=0,1,2,\dots,n) 求得一个唯一的多项式:

Ln(x)=a0+a1x+a2x2++anxnL_{n}(x)=a_{0}+a_{1} x+a_{2} x^{2}+\ldots+a_{n} x^{n}

使得 Ln(xi)=yiL_n(x_i)=y_i ,即

[1x0x0n1x1x1n1xnxnn][a0a1an]=[y0y1yn]\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]

对于上述范德蒙行列式,易得,[a1,a2,][a_1,a_2,\dotsb] 有且仅有一组解。

拉格朗日插值法

上多项式组易解得

L(x)=y1(xx2)(xx3)...(xxn)(x1x2)(x1x3)...(x1xn)+y2(xx1)(xx3)...(xxn)(x2x1)(x2x3)...(x2xn)...+yn(xx1)(xx2)...(xxn1)(xnx1)(xnx2)...(xnxn1)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)=i=1nyij=1,jinxxjxixjL(x)=\sum_{i=1}^{n}y_i \prod_{j=1,j \neq i}^n \frac{x-x_j}{x_i-x_j}

为拉格朗日插值法。

:拉格朗日插值公式在理论分析理解上很容易理解,但是若插值节点发生改变时,插值公式随之就要重新计算生成,在实际计算中会占用大量的计算量。牛顿法的出现正是克服这个问题[3]

牛顿插值法

f(x)=f(x0)+f[x0,x1](xx0)+f[x0,x1,x2](xx0)(xx1)++f[x0,x1,,xn2,xn1](xx0)(xx1)(xxn3)(xxn2)+f[x0,x1,,xn1,xn](xx0)(xx1)(xxn2)(xxn1)\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[xi,xj]=f(xi)f(xj)xixj(ij,xixj)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[xi,xj,xk]=f[xi,xj]f[xj,xk]xixk(ik)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[x0,x1,Λxn]=f[x0,x1,Λ,xn1]f[x1,x2,Λxn]x0xnf\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}}

由以上定义我们的到差商表如下:

xiyi 一阶差商  二阶差商 n 阶差商 x0f(x0)x1f(x1)f[x0,x1]x2f(x2)f[x1,x2]f[x0,x1,x2]xn1f(xn1)xnf(xn)f[xn1,xn]f[xn2,xn1,xn]f[xn,xn]\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}

我们可以得到下面公式:

f(x)=f(xθ)+f[x,xθ](xxθ),f[x,x0]=f[x0,x1]+f[x,x0,x1](xx1),f[x,x0,x1]=f[x0,x1,x2]+f[x,x0,x1,x2](xx2),f[x,xθ,,xm1]=f[xθ,x1,,xn]+J[x,xθ,,xn](xxn)\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)f(x) 在揷值区间 [p,q][p, q] 上的 n+1n+1 个互不相同的节点 xi(i=0,1,,n)x_{i}(i=0,1, \ldots, n) 处满足 f(xi)=fif\left(x_{i}\right)=f_{i}f(xi)=fi(i=0,1,2,,n)f^{\prime}\left(x_{i}\right)=f_{i}^{\prime}(i=0,1,2, \ldots, n) , 如果函数 G(x)G(x) 的存在满足下列条件:

  1. G(x)G(x) 在每个小区间上的多项式次数为 3 ;
  2. G(x)C1[a,b]G(x) \in C^{1}[a, b] ;
  3. G(xi)=f(xi),G(xi)=f(xi),i=(0,1,,n)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)G(x)f(x)f(x)n+1n+1 个节点 xix_{i} 上的分段三次埃尔米特插值多项式。

G(x)=hkyk(x)+hk+1yk+1(x)+Hk(x)yk+Hk+1(x)yk+1=(1+2xxkxk+1xk)(xxk+1xkxk+1)2yk+(1+2xxk+1xkxk+1)(xxkxk+1xk)2yk+1+(xxk)(xxk+1xkxk+1)2yk+(xxk+1)(xxkxk+1xk)2yk+1\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)=i=0nhi[(xix)(2aiyiyi)+yi]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]

其中:

hi=j=0jin(xxjxixj)2,ai=j=0jin1xixjh_{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 的数组下标从 11 开始),mm 个插值点以数组 x 输入,输出数组 ymm 个插值。编

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;
end

调用Matlab封装函数

x = -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]

数学上将具有一定光滑性的分段多项式称为样条函数。

对于 Δ:a=x0<x1<<xn1<xn=b\Delta : \quad a=x_{0}<x_{1}<\cdots<x_{n-1}<x_{n}=b 满足:

  1. 在每个小区间 [xi,xi1](i=0,1,,n1)\left[x_{i}, x_{i-1}\right](i=0,1, \cdots, n-1)s(x)s(x)kk 次多项式。
  2. s(x)s(x)[a,b][a, b] 上具有 k1k-1 阶连续导数。即 Si(j)(xi+1)=Si+1(j)(xi+1)S_{i}^{(j)}(x_{i+1})=S_{i+1}^{(j)}(x_{i+1}) (j=0,1,,k1)(j=0,1,\dots,k-1)

一般形式为:

sk(x)=i=0kαixii!+j=1n1βjk!(xxj)+ks_{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') %标注显示在东南方向

参考