插值相关总结
多项式插值
$$\det (A) = \left| {\begin{array}{*{20}{c}}
1&{{x_0}}&{x_0^2}& \cdots &{x_0^n}&{}\\
1&{{x_1}}&{x_1^2}& \cdots &{x_1^n}&{}\\
{}& \cdots & \cdots & \cdots & \cdots & \cdot \\
1&{{x_n}}&{x_n^2}& \cdots &{x_n^n}&{}
\end{array}} \right|{\rm{ }}! = {\rm{ 0}}$$
拉格朗日插值
1779年英国数学家爱德华·华林于最早发现拉格朗日插值法。不久后1783年，由莱昂哈德·欧拉再次发现。直到1795年，拉格朗日在其著作《师范学校数学基础教程》中发表了这个插值方法。
$${l_i}(x) = \frac{{\left( {x - {x_0}} \right) \cdots \left( {x - {x_{i - 1}}} \right)\left( {x - {x_{i + 1}}} \right) \cdots \left( {x - {x_n}} \right)}}{{\left( {{x_i} - {x_0}} \right) \cdots \left( {{x_i} - {x_{i - 1}}} \right)\left( {{x_i} - {x_{i + 1}}} \right) \cdots \left( {{x_i} - {x_n}} \right)}} = \prod\limits_{\begin{array}{*{20}{c}}
{j = 0}\\
{j \ne i}
\end{array}}^n {\frac{{x - {x_j}}}{{{x_i} - {x_j}}}} ,\quad (i = 0,1, \cdots ,n)$$
注：$\frac{{x - {x_j}}}{{{x_i} - {x_j}}}$就是求解$p(x)$的$n+1$个基函数。如果对推导过程感兴趣，参考链接里有，很简单的。
$${l_i}\left( {{x_j}} \right) = \left\{ {\begin{array}{*{20}{l}}
0&{j \ne i}\\
1&{j = i}
\end{array}} \right.$$
令
$${L_n}(x) = \sum\limits_{i = 0}^n {{y_i}} {l_i}(x) = \sum\limits_{i = 0}^n {{y_i}} \left( {\prod\limits_{\begin{array}{*{20}{c}}
{j = 0}\\
{j = i}
\end{array}}^n {\frac{{x - {x_j}}}{{{x_i} - {x_j}}}} } \right)$$
这就是$n$次的拉格朗日插值多项式，且$n+1$的$n$次拉格朗日插值多项式是惟一的。
function [ yi ] = lagrange_interp (X,Y,xi) n=length(X); %得到已知数据长度 m=length(xi); %得到待插值数据长度 yi=zeros(size(xi)); for j=1:m %待插值数据有m个，计算每个插值结果 for i=1:n %已知的n个数据构造中间值 temp=1; %temp用于存储中间值 for k=1:n if(i~=k) %和自身标号相同的不相乘 temp=temp*(xi(j)-X(k))/(X(i)-X(k)); end end yi(j)=Y(i)*temp+yi(j); end end end
牛顿插值
1. 差商
2. 牛顿插值公式
function yi=newton_interp(X,Y,xi) syms t; %定义自变量t，用于字符公式 if(length(X)==length(Y)) n=length(X); c(1:n)=0.0; else disp('X和Y的维数不相等！'); return; end f=Y(1); %f用来记录得到的牛顿插值公式的字符串表达式 l=1; for i=1:n-1 y1=zeros(1,n-i); for j=i+1:n y1(j)=(Y(j)-Y(i))/(X(j)-X(i)); end c(i)=y1(i+1); %c记录差分 l=l*(t-X(i)); %l记录(x-x0)(x-x1)……的值 f=f+c(i)*l; %累加得到差分公式 Y=y1; end f=simplify(f); %简化得到的牛顿插值公式 m=length(xi); %开始输出 for i=1:m yi(i)=subs(f,'t',xi(i)); % 根据公式计算需要的值 end yi=double(yi); % 转换为数值型，为返回值 end
对比拉格朗日插值法：与拉格朗日插值法相比，牛顿插值法增加节点时，不需要重新计算，在某些情景下比拉格朗日插值法更好用，而且在计算余项时，牛顿插值的余项由于不需要导数，故$f(x)$是由离散点或者导数不存在时仍然适用，这是拉格朗日余项计算所不能比拟的。差商与导数的关系为：$$f\left[x_{0}, x_{1}, \cdots, x_{n}\right]=\frac{f^{(n)}(\xi)}{n !}$$其中，$\xi \in(\alpha, \beta), \alpha=\min _{0 \leq i \leq n}\left\{x_{i}\right\}, \beta=\max _{0 \leq i \leq n}\left\{x_{i}\right\}$
泰勒插值
埃尔米特插值
不少实际的插值问题不但要求在节点上的函数值相等，而且还要求对应的导数值也相等，甚至要求高阶导数也相等，满足这种要求的插值多项式就是埃尔米特插值多项式。
这意味着，拟合曲线与实际曲线不仅都过点$\left(x_{i}, y_{i}\right),\left(x_{i+1}, y_{i+1}\right)$，而且两点处还有相同的切线。
设在节点$a \leq x_{0}<x_{1}<\cdots<x_{n} \leq b$上，$y_{j}=f\left(x_{j}\right)$，$m_{j}=f^{\prime}\left(x_{j}\right) \quad(j=0,1, \cdots, n)$
即满足条件：$H\left(x_{j}\right)=y_{j}, \quad H^{\prime}\left(x_{j}\right)=m_{j}, \quad(j=0,1, \cdots, n)$
这里共有$2n+2$个插值条件，可唯一确定一个次数不超过$2n+1$的多项式$H_{2 n+1}(x)=H(x)$，其形式为：
$$H_{2 n+1}(x)=a_{0}+a_{1} x+\cdots+a_{2 n+1} x^{2 n+1}$$
采用拉格朗日插值多项式的基函数方法的思想：
先求出$2n+2$个插值基函数，每个基函数都是$2n+1$次多项式，且满足条件：
$$\begin{array}{l}{\alpha_{j}\left(x_{k}\right)=\delta_{j k}=\left\{\begin{array}{ll}{0,} & {j \neq k} \\ {1,} & {j=k}\end{array}, \quad \alpha_{j}^{\prime}\left(x_{k}\right)=0\right.} \\ {\beta_{j}\left(x_{k}\right)=0, \quad \beta_{j}^{\prime}\left(x_{k}\right)=\delta_{j k} \quad(j, k=0,1, \cdots, n)}\end{array}$$
用基函数表示插值多项式：
$$H_{2 n+1}(x)=\sum_{j=0}^{n}\left[y_{j} \alpha_{j}(x)+m_{j} \beta_{j}(x)\right]$$
由插值基函数所满足的条件，有$H_{2 n+1}\left(x_{k}\right)=y_{k}, \quad H_{2 n+1}^{\prime}\left(x_{k}\right)=m_{k}, \quad(k=0,1, \cdots, n)$
然后求出基函数$\alpha_{j}(x)$与$\beta_{j}(x)$
利用拉格朗日插值基函数$l_{j}(x)$，
令$\alpha_{j}(x)=(a x+b) l_{j}^{2}(x)$
得到
$$\begin{array}{l}{\alpha_{j}\left(x_{j}\right)=\left(a x_{j}+b\right) l_{j}^{2}\left(x_{j}\right)=1} \\ {\alpha_{j}^{\prime}\left(x_{j}\right)=l_{j}\left(x_{j}\right)\left[a l_{j}\left(x_{j}\right)+2\left(a x_{j}+b\right) l_{j}^{\prime}\left(x_{j}\right)\right]=0}\end{array}$$
整理得
$$\left\{\begin{array}{l}{a x_{j}+b=1} \\ {a+2 l_{j}^{\prime}\left(x_{j}\right)=0}\end{array}\right.$$
解出
$$a=-2 l_{j}^{\prime}\left(x_{j}\right), \quad b=1+2 x_{j} l_{j}^{\prime}\left(x_{j}\right)$$
由于
$$l_{j}(x)=\frac{\left(x-x_{0}\right) \cdots\left(x-x_{j-1}\right)\left(x-x_{j+1}\right) \cdots\left(x-x_{n}\right)}{\left(x_{j}-x_{0}\right) \cdots\left(x_{j}-x_{j-1}\right)\left(x_{j}-x_{j+1}\right) \cdots\left(x_{j}-x_{n}\right)}$$
两边取对数再求导，得
$$l_{j}^{\prime}\left(x_{j}\right)=\sum_{k=0 \atop k \neq j}^{n} \frac{1}{x_{j}-x_{k}}$$
于是，
$$\alpha_{j}(x)=\left(1-2\left(x-x_{j}\right) \sum_{k=0 \atop k \neq j}^{n} \frac{1}{x_{j}-x_{k}}\right) l_{j}^{2}(x)$$
同理，得到
$$\beta_{j}(x)=\left(x-x_{j}\right) l_{j}^{2}(x)$$
综上所述，埃尔米特插值多项式为：
$$\left\{ {\begin{array}{*{20}{l}}
{{H_{2n + 1}}(x) = \sum\limits_{j = 0}^n {\left[ {{y_j}{\alpha _j}(x) + {m_j}{\beta _j}(x)} \right]} }\\
{{y_j} = f\left( {{x_j}} \right),{m_j} = {f^\prime }\left( {{x_j}} \right)}\\
{{\alpha _j}(x) = \left( {1 - 2\left( {x - {x_j}} \right)\sum\limits_{k = 0}^n {\frac{1}{{{x_j} - {x_k}}}} } \right)l_j^2(x),k \ne j}\\
{{\beta _j}(x) = \left( {x - {x_j}} \right)l_j^2(x)}
\end{array}} \right..$$
仿照拉格朗日插值余项的证明方法，得到埃尔米特插值余项：
$$R(x)=f(x)-H_{2 n+1}(x)=\frac{f^{(2 n+2)}(\xi)}{(2 n+2) !} \omega_{n+1}^{2}(x)$$
$n$取1时，是一个非常重要的特例，即两点三次埃尔米特插值，比较典型和常用(4个基函数均为三次多项式，$n=1$)
插值基函数满足条件：
$$\begin{array}{l}{\alpha_{k}\left(x_{k}\right)=1, \quad \alpha_{k}\left(x_{k+1}\right)=0, \quad \alpha_{k}^{\prime}\left(x_{k}\right)=\alpha_{k}^{\prime}\left(x_{k+1}\right)=0} \\ {\alpha_{k+1}\left(x_{k}\right)=0, \quad \alpha_{k+1}\left(x_{k+1}\right)=1, \quad \alpha_{k+1}^{\prime}\left(x_{k}\right)=\alpha_{k+1}^{\prime}\left(x_{k+1}\right)=0} \\ {\beta_{k}\left(x_{k}\right)=\beta_{k}\left(x_{k+1}\right)=0, \beta_{k}^{\prime}\left(x_{k}\right)=1, \quad \beta_{k}^{\prime}\left(x_{k+1}\right)=0} \\ {\beta_{k+1}\left(x_{k}\right)=\beta_{k+1}\left(x_{k+1}\right)=0, \beta_{k+1}\left(x_{k}\right)=0, \quad \beta_{k+1}^{\prime}\left(x_{k+1}\right)=1}\end{array}$$
两点三次埃尔米特插值多项式为：
$$H_{3}(x)=y_{k} \alpha_{k}(x)+y_{k+1} \alpha_{k+1}(x)+m_{k} \beta_{k}(x)+m_{k+1} \beta_{k+1}(x)$$
余项为：
$$R_{3}(x)=\frac{1}{4 !} f^{(4)}(\xi)\left(x-x_{k}\right)^{2}\left(x-x_{k+1}\right)^{2}, \quad \xi \in\left(x_{k}, x_{k+1}\right)$$
function y-hermite(x0,y0,y1,x); %x0,y0为样本点数据，y1为导数指，m个插值点以数组x输入，输出数组y为m个插值 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=i: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)-x0(k))*(2*a*y0(i)-y(i))+y0(i)); end y(k)=yy; end
分段插值
分段线性插值
$${l_i}(x) = \left\{ {\begin{array}{*{20}{l}}
{\frac{{x - {x_{i - 1}}}}{{{x_i} - {x_{i - 1}}}},}&{x \in \left[ {{x_{i - 1}},{x_i}} \right]{\rm{ }}i = 0}\\
{\frac{{x - {x_{i + 1}}}}{{{x_i} - {x_{i + 1}}}},}&{x \in \left[ {{x_i},{x_{i + 1}}} \right]{\rm{ }}i = n}\\
{0,}&{else}
\end{array}{\rm{ }}} \right.$$
MATLAB：
method：
| nearest | 最近项插值 |
| linear | 线性插值 |
| spline | 逐段3次样条插值 |
| cubic | 保凹凸性3次插值 |
n=11; m=61; x=-5:10/(m-1):5; y=1./(1+x.^2); z=0*x; x0=-5:10/(n-1):5; y0=1./(1+x0.^2); y1=interp1(x0, y0, x);%interp1(x0,y0,x)为Matlab中现成的分段线性插值程序 plot(x, z, 'r', x, y, 'k:', x, y1, 'r') gtext('Piece. –linear.'), gtext('y=1/(1+x^2)') title('Piecewise Linear')
分段三次埃尔米特插值
根据两点三次埃尔米特插值多项式，得到${I_h}(x)$在$\left[x_{k}, x_{k+1}\right]$的表达式为：
若在整个区间$[a, b]$上定义一组分段三次插值基函数$\alpha_{j}(x)$及$\beta_{j}(x)$,$(j=0,1, \cdots, n)$，则$I_{h}(x)$可表示为：
$${\alpha _j}(x) = \left\{ {\begin{array}{*{20}{c}}
{{{\left( {\frac{{x - {x_{j - 1}}}}{{{x_j} - {x_{j - 1}}}}} \right)}^2}\left( {1 + 2\frac{{x - {x_j}}}{{{x_{j - 1}} - {x_j}}}} \right){x_{j - 1}} \le x \le {x_j}}\\
\begin{array}{l}
{\left( {\frac{{x - {x_{j + 1}}}}{{{x_j} - {x_{j + 1}}}}} \right)^2}\left( {1 + 2\frac{{x - {x_j}}}{{{x_{j + 1}} - {x_j}}}} \right){x_j} \le x \le {x_{j + 1}}\\
{\rm{ }}0
\end{array}
\end{array}} \right.$$
$${\beta _j}(x) = \left\{ {\begin{array}{*{20}{c}}
{{{\left( {\frac{{x - {x_{j - 1}}}}{{{x_j} - {x_{j - 1}}}}} \right)}^2}\left( {x - {x_j}} \right){x_{j - 1}} \le x \le {x_j}}\\
\begin{array}{l}
{\left( {\frac{{x - {x_{j + 1}}}}{{{x_j} - {x_{j + 1}}}}} \right)^2}\left( {x - {x_j}} \right){x_j} \le x \le {x_{j + 1}}\\
{\rm{ }}0
\end{array}
\end{array}} \right.$$
MATLAB：来自官方文档 https://ww2.mathworks.cn/help/matlab/ref/pchip.html
spline 和 pchip 为两个不同函数生成的插值结果进行比较。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); %逐段3次样条插值
plot(x,y,'o',xq1,p,'-',xq1,s,'-.')
legend('Sample Points','pchip','spline','Location','SouthEast')
三次样条插值
分段插值的优点是误差小，稳定性高，但是曲线的光滑性不好，而许多实际问题需要插值曲线有较高阶的整体光滑性，理论上是可以使用高次埃尔米特插值或分段高次埃尔米特插值，但是必须知道每一个点的导数是不现实的。所以，样条插值诞生了。理论推导参考：https://www.cnblogs.com/duye/p/8671820.html1、y=interp1(${x_0}$,${y_0}$,$x$,'spline');
% spline改成linear，则变成线性插值
2、y=spline(${x_0}$,${y_0}$,$x$);
% 这个是根据己知的x，y数据，用样条函数插值出${x_i}$处的值。
% 由己知的x，y数据，求出它的样条函数表达式，不过该表达式不是用矩阵直接表示，要求点x`的值，要用函数
3、pp=csape(x,y,conds); y=ppval(pp,$x$);
% 各种边界条件的三次样条插值
conds是指边界条件，边界类型可为:
①'complete'： 给定边界一阶导数，默认的边界条件
②'not-a-knot'：非扭结条件，不用给边界值.
③'periodic'： 周期性边界条件，不用给边界值.
④'second'： 给定边界二阶导数.
⑤'variational'： 自然样条(边界二阶导数为[0,0]。
clear clc x0=[0,3,5,7,9,11,12,13,14,15]; y0=[0,1.2,1.7,2.0,2.1,2.0,1.8,1.2,1.0,1.6]; t=0:0.05:15; y1=spline(x0,y0,t); dy1=(spline(x0,y0,0.0001)-spline(x0,y0,0))/0.0001%x=0处斜率 min1=min(spline(x0,y0,13:0.001:15))%13到15最小值 figure plot(x0,y0,'ro',t,y1);%画出曲线 title('三次样条插值');
二维插值
网络节点插值
z=interp2(x0,y0,z0,x,y,'method')
x0，y0是节点坐标(分别为行向量和列向量)；
z0是节点的值；
method和前面所述相同。
$x$，$y$为插值点坐标，$z$为函数值。
如果是三次样条插值
pp=csape({x0,y0},z0,conds)
z=fnval(pp,{x,y})
“conds”与一维相同，一般默认。
clear clc x=100:100:400; y=100:100:400; z=[636,697,624,478; 712,630,478,420; 674,598,412,400; 626,552,334,310]; pp=csape({x,y},z'); %z矩阵的行列对应向量x,y xi=100:10:500; yi=100:10:400; cz1=fnval(pp,{xi,yi}); cz2=interp2(x,y,z,xi,yi','spline'); [i,j]=find(cz1==max(max(cz1))) subplot(1,2,1); surf(xi,yi,cz1'); shading interp; %插入颜色插值 axis equal; title('cz1'); subplot(1,2,2); surf(xi,yi,cz2); shading interp; axis equal; title('cz2');
散乱节点插值
zi=griddata(x,y,z,xi,yi)
clear clc %样本点信息 x=[129,140,103.5,88,185.5,95,105,157.5,107.5,77,81,162,162,117.5]; y=[7.5,141.5,23,147,22.5,137.5,85.5,-6.5,-81,3,56.5,-66.5,84,-33.5]; z=-[4,8,6,8,6,8,8,9,9,8,8,9,4,9]; xi=75:200; yi=-50:150; zi=griddata(x,y,z,xi,yi','cubic'); subplot(1,2,1); plot(x,y,'*'); title('xy'); subplot(1,2,2); mesh(xi,yi,zi); shading interp; axis equal; title('xyz');
