多项式插值
$$det (A) = left| {egin{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}}
ight|{
m{ }}! = {
m{ 0}}$$
拉格朗日插值
1779年英国数学家爱德华·华林于最早发现拉格朗日插值法。不久后1783年，由莱昂哈德·欧拉再次发现。直到1795年，拉格朗日在其著作《师范学校数学基础教程》中发表了这个插值方法。
$${l_i}(x) = frac{{left( {x - {x_0}}
ight) cdots left( {x - {x_{i - 1}}}
ight)left( {x - {x_{i + 1}}}
ight) cdots left( {x - {x_n}}
ight)}}{{left( {{x_i} - {x_0}}
ight) cdots left( {{x_i} - {x_{i - 1}}}
ight)left( {{x_i} - {x_{i + 1}}}
ight) cdots left( {{x_i} - {x_n}}
ight)}} = prodlimits_{egin{array}{*{20}{c}}
{j = 0}\
{j
e 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}}
ight) = left{ {egin{array}{*{20}{l}}
0&{j
e i}\
1&{j = i}
end{array}}
ight.$$
令
$${L_n}(x) = sumlimits_{i = 0}^n {{y_i}} {l_i}(x) = sumlimits_{i = 0}^n {{y_i}} left( {prodlimits_{egin{array}{*{20}{c}}
{j = 0}\
{j = i}
end{array}}^n {frac{{x - {x_j}}}{{{x_i} - {x_j}}}} }
ight)$$
这就是$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)$是由离散点或者导数不存在时仍然适用，这是拉格朗日余项计算所不能比拟的。差商与导数的关系为：$$fleft[x_{0}, x_{1}, cdots, x_{n} ight]=frac{f^{(n)}(xi)}{n !}$$其中，$xi in(alpha, eta), alpha=min _{0 leq i leq n}left{x_{i} ight}, eta=max _{0 leq i leq n}left{x_{i} ight}$
泰勒插值
埃尔米特插值
不少实际的插值问题不但要求在节点上的函数值相等，而且还要求对应的导数值也相等，甚至要求高阶导数也相等，满足这种要求的插值多项式就是埃尔米特插值多项式。
这意味着，拟合曲线与实际曲线不仅都过点$left(x_{i}, y_{i} ight),left(x_{i+1}, y_{i+1} ight)$，而且两点处还有相同的切线。
设在节点$a leq x_{0}<x_{1}<cdots<x_{n} leq b$上，$y_{j}=fleft(x_{j} ight)$，$m_{j}=f^{prime}left(x_{j} ight) quad(j=0,1, cdots, n)$
即满足条件：$Hleft(x_{j} ight)=y_{j}, quad H^{prime}left(x_{j} ight)=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$次多项式，且满足条件：
$$egin{array}{l}{alpha_{j}left(x_{k} ight)=delta_{j k}=left{egin{array}{ll}{0,} & {j eq k} \ {1,} & {j=k}end{array}, quad alpha_{j}^{prime}left(x_{k} ight)=0 ight.} \ {eta_{j}left(x_{k} ight)=0, quad eta_{j}^{prime}left(x_{k} ight)=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} eta_{j}(x) ight]$$
由插值基函数所满足的条件，有$H_{2 n+1}left(x_{k} ight)=y_{k}, quad H_{2 n+1}^{prime}left(x_{k} ight)=m_{k}, quad(k=0,1, cdots, n)$
然后求出基函数$alpha_{j}(x)$与$eta_{j}(x)$
利用拉格朗日插值基函数$l_{j}(x)$，
令$alpha_{j}(x)=(a x+b) l_{j}^{2}(x)$
得到
$$egin{array}{l}{alpha_{j}left(x_{j} ight)=left(a x_{j}+b ight) l_{j}^{2}left(x_{j} ight)=1} \ {alpha_{j}^{prime}left(x_{j} ight)=l_{j}left(x_{j} ight)left[a l_{j}left(x_{j} ight)+2left(a x_{j}+b ight) l_{j}^{prime}left(x_{j} ight) ight]=0}end{array}$$
整理得
$$left{egin{array}{l}{a x_{j}+b=1} \ {a+2 l_{j}^{prime}left(x_{j} ight)=0}end{array} ight.$$
解出
$$a=-2 l_{j}^{prime}left(x_{j} ight), quad b=1+2 x_{j} l_{j}^{prime}left(x_{j} ight)$$
由于
$$l_{j}(x)=frac{left(x-x_{0} ight) cdotsleft(x-x_{j-1} ight)left(x-x_{j+1} ight) cdotsleft(x-x_{n} ight)}{left(x_{j}-x_{0} ight) cdotsleft(x_{j}-x_{j-1} ight)left(x_{j}-x_{j+1} ight) cdotsleft(x_{j}-x_{n} ight)}$$
两边取对数再求导，得
$$l_{j}^{prime}left(x_{j} ight)=sum_{k=0 atop k eq j}^{n} frac{1}{x_{j}-x_{k}}$$
于是，
$$alpha_{j}(x)=left(1-2left(x-x_{j} ight) sum_{k=0 atop k eq j}^{n} frac{1}{x_{j}-x_{k}} ight) l_{j}^{2}(x)$$
同理，得到
$$eta_{j}(x)=left(x-x_{j} ight) l_{j}^{2}(x)$$
综上所述，埃尔米特插值多项式为：
$$left{ {egin{array}{*{20}{l}}
{{H_{2n + 1}}(x) = sumlimits_{j = 0}^n {left[ {{y_j}{alpha _j}(x) + {m_j}{eta _j}(x)}
ight]} }\
{{y_j} = fleft( {{x_j}}
ight),{m_j} = {f^prime }left( {{x_j}}
ight)}\
{{alpha _j}(x) = left( {1 - 2left( {x - {x_j}}
ight)sumlimits_{k = 0}^n {frac{1}{{{x_j} - {x_k}}}} }
ight)l_j^2(x),k
e j}\
{{eta _j}(x) = left( {x - {x_j}}
ight)l_j^2(x)}
end{array}}
ight..$$
仿照拉格朗日插值余项的证明方法，得到埃尔米特插值余项：
$$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$)
插值基函数满足条件：
$$egin{array}{l}{alpha_{k}left(x_{k} ight)=1, quad alpha_{k}left(x_{k+1} ight)=0, quad alpha_{k}^{prime}left(x_{k} ight)=alpha_{k}^{prime}left(x_{k+1} ight)=0} \ {alpha_{k+1}left(x_{k} ight)=0, quad alpha_{k+1}left(x_{k+1} ight)=1, quad alpha_{k+1}^{prime}left(x_{k} ight)=alpha_{k+1}^{prime}left(x_{k+1} ight)=0} \ {eta_{k}left(x_{k} ight)=eta_{k}left(x_{k+1} ight)=0, eta_{k}^{prime}left(x_{k} ight)=1, quad eta_{k}^{prime}left(x_{k+1} ight)=0} \ {eta_{k+1}left(x_{k} ight)=eta_{k+1}left(x_{k+1} ight)=0, eta_{k+1}left(x_{k} ight)=0, quad eta_{k+1}^{prime}left(x_{k+1} ight)=1}end{array}$$
两点三次埃尔米特插值多项式为：
$$H_{3}(x)=y_{k} alpha_{k}(x)+y_{k+1} alpha_{k+1}(x)+m_{k} eta_{k}(x)+m_{k+1} eta_{k+1}(x)$$
余项为：
$$R_{3}(x)=frac{1}{4 !} f^{(4)}(xi)left(x-x_{k} ight)^{2}left(x-x_{k+1} ight)^{2}, quad xi inleft(x_{k}, x_{k+1} ight)$$
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{ {egin{array}{*{20}{l}}
{frac{{x - {x_{i - 1}}}}{{{x_i} - {x_{i - 1}}}},}&{x in left[ {{x_{i - 1}},{x_i}}
ight]{
m{ }}i = 0}\
{frac{{x - {x_{i + 1}}}}{{{x_i} - {x_{i + 1}}}},}&{x in left[ {{x_i},{x_{i + 1}}}
ight]{
m{ }}i = n}\
{0,}&{else}
end{array}{
m{ }}}
ight.$$
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} ight]$的表达式为：
若在整个区间$[a, b]$上定义一组分段三次插值基函数$alpha_{j}(x)$及$eta_{j}(x)$,$(j=0,1, cdots, n)$，则$I_{h}(x)$可表示为：
$${alpha _j}(x) = left{ {egin{array}{*{20}{c}}
{{{left( {frac{{x - {x_{j - 1}}}}{{{x_j} - {x_{j - 1}}}}}
ight)}^2}left( {1 + 2frac{{x - {x_j}}}{{{x_{j - 1}} - {x_j}}}}
ight){x_{j - 1}} le x le {x_j}}\
egin{array}{l}
{left( {frac{{x - {x_{j + 1}}}}{{{x_j} - {x_{j + 1}}}}}
ight)^2}left( {1 + 2frac{{x - {x_j}}}{{{x_{j + 1}} - {x_j}}}}
ight){x_j} le x le {x_{j + 1}}\
{
m{ }}0
end{array}
end{array}}
ight.$$
$${eta _j}(x) = left{ {egin{array}{*{20}{c}}
{{{left( {frac{{x - {x_{j - 1}}}}{{{x_j} - {x_{j - 1}}}}}
ight)}^2}left( {x - {x_j}}
ight){x_{j - 1}} le x le {x_j}}\
egin{array}{l}
{left( {frac{{x - {x_{j + 1}}}}{{{x_j} - {x_{j + 1}}}}}
ight)^2}left( {x - {x_j}}
ight){x_j} le x le {x_{j + 1}}\
{
m{ }}0
end{array}
end{array}}
ight.$$
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');
