Matalab插值详解和源码

 转载：Matalab插值详解和源码 - 知乎 (zhihu.com)

插值法

插值法又称“内插法”，是利用函数f (x)在某区间中已知的若干点的函数值，作出适当的特定函数，在区间的其他点上用这特定函数的值作为函数f (x)的近似值，这种方法称为插值法。如果这特定函数是多项式，就称它为插值多项式。常见分段线性插值法和样条差值，样条插值误差更小。

1 线性插值法

线性插值法是指使用连接两个已知量的直线来确定在这两个已知量之间的一个未知量的值的方法。 假设我们已知坐标(x0,y0)与(x1,y1),要得到[x0,x1]区间内某一位置x在直线上的y值。根据图中所示，我们得到两点式直线方程：

假设方程两边的值为α,则：

则：y = y0 + α(y1 − y0)(即已知x就可得a值 然后再得y值。

实际上，即使x不在x0到x1之间并且α也不是介于0到1之间，这个公式也是成立的。在这种情况下，这种方法叫作线性外插—参见外插值。

已知y求x的过程与以上过程相同。

2 样条插值法

常用三次样条插值法，具体概念比较复杂，可以自行百度。

3 拉格朗日差值法

略。

• MATLAB中的使用-interp1函数：

MATLAB中的插值函数为interp1，其调用格式为： yi=interp1(x,y,xi,'method') 其中x，y为插值点，yi为在被插值点xi处的插值结果；x,y为向量，
'method'表示采用的插值方法，MATLAB提供的插值方法有几种：
'nearest'是最邻近插值， 'linear'线性插值； 'spline'三次样条插值； 'pchip'立方插值．缺省时表示线性插值
注意：所有的插值方法都要求x是单调的，并且xi不能够超过x的范围。

x = 0:2*pi;    
y = sin(x);    
xx = 0:0.5:2*pi;    
  
% interp1对sin函数进行分段线性插值，调用interp1的时候，默认的是分段线性插值    
y1 = interp1(x,y,xx,'linear');    
subplot(2,2,1);  
plot(x,y,'o',xx,y1,'r')    
title('分段线性插值')    
    
% 临近插值    
y2 = interp1(x,y,xx,'nearest');    
subplot(2,2,2);  
plot(x,y,'o',xx,y2,'r');    
title('临近插值')    
    
%球面线性插值    
y3 = interp1(x,y,xx,'spline');    
subplot(2,2,3);  
plot(x,y,'o',xx,y3,'r')    
title('球面插值')    
    
%三次多项式插值法    
y4 = interp1(x,y,xx,'pchip');    
subplot(2,2,4);  
plot(x,y,'o',xx,y4,'r');    
title('三次多项式插值')    

例：环境温度

例如：在一 天24小时内，从零点开始每间隔2小时测得的环境温度数据分别为
12，9，9，1，0，18 ，24，28，27，25，20，18，15，13，
推测中午1点（即13点）时的温度．

x = 0:2:24;  
y = [12   9   9   10   18  24   28   27   25   20  18  15  13];  
a = 13;  
y1 = interp1(x,y,a,'spline')  
% 结果为：  27.8725  
 
% 若要得到一天24小时的温度曲线，则：  
xi = 0:1/3600:24;  
% 插值点可以是向量，则返回的也就是对应的向量  
yi = interp1(x,y,xi, 'spline');  
plot(x,y,'o' ,xi,yi);  

一天24小时温度曲线

• Csape函数

csape函数的用法如下：
pp = csape(x,y,conds,valconds)
其中(x,y)为数据向量，conds表示变界类型， valconds表示边界值。
边界类型(conds)可为：
'complete'，给定边界 一阶导数.
'not-a-knot'，非扭结条件,不用给边界值.
'periodic'， 周期性边界条件,不用给边界值.
'second'，给定边界二阶导数.
'variational'，自然样条(边界二阶导数为0)
边界类型(valconds)可为：
'complete'，给定边界 一阶导数.
'not-a-knot'，非扭结条件,不用给边界值.
'periodic'， 周期性边界条件,不用给边界值.
'second'，给定边界二阶导数.
'variational'，自然样条(边界二阶导数为0)

例：机床加工

待加工零件的外形根据工艺要求由一组数据(x,y)给出，用程控铣床加工时每一刀只能沿x方向和y方向走非常小的一步，这就需要从已知数据得到加工所要求的步长很小的(x,y)坐标。

依据下表给出的(x,y)数据求x坐标每改变0.1时的y坐标。试完成加工所需数据，画出曲线，并求出x=0处的曲线斜率和13<=x<=15范围内y的最小值。

x 0 3 5 7 9 11 12 13 14 15
y 0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6

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];
x=0:0.1:15;
y1=interp1(x0,y0,x);%分段线性插值 得到x=0:0.1:15相应的函数值向量y1
y2=interp1(x0,y0,x,'spline');%三次样条插值 得到x=0:0.1:15相应的函数值向量y2
pp1=csape(x0,y0);%csape函数三次样条插值，默认边界条件为拉格朗日边界条件 
y3=fnval(pp1,x);%得到x=0:0.1:15相应的函数值向量y3
pp2=csape(x0,y0,'second');%设置边界条件为二阶导数，默认为[0,0]
y4=fnval(pp2,x);%得到x=0:0.1:15相应的函数值向量y4
[x',y1',y2',y3',y4'];%得到一个151*4的矩阵可以观察y1,y2,y3,y4随x的变化
subplot(1,3,1)%将当前图窗分为1*3网络并在1处创建坐标区（即当前选定1区绘图）
plot(x0,y0,'+',x,y1);%(x0,y0)所代表的点用+表示
title('分段线性插值');
subplot(1,3,2)%选定2区绘图
plot(x0,y0,'+',x,y2);
title('三次样条插值-interp1函数')
subplot(1,3,3)%选定3区绘图
plot(x0,y0,'+',x,y3);
title('三次样条插值-scape函数')
%求x=0处斜率
dx=diff(x);%dx即∆x
dy=diff(y3);%dy即∆y
dy_dx=dy./dx;%(./ 对应元素相除)
dy_dx0=dy_dx(1);%x=0处的导数就是dy./dx的第一个元素:dy0/dx0
%求13<=x<=15范围内y的最小值
ytemp=y3(131:151);%取出13<=x<=15区间内对应的y3值（因为y3误差更小，数据分析更准确）
ymin=min(ytemp);
index=find(y3==ymin);%y3中最小值对应的序列
xmin=x(index);%并由序列找到对应的x值
[xmin,ymin]%输出

运行程序：

可以看出调用interp1和csape三次样条插值结果相同，同时三次样条插值比分段线性插值更精确。

最小值：
ans =  13.8000    0.9851

csape 和interp1有什么相同和区别？

csape和interp1都是插值函数。
csape可以选择样条的边界条件，interp1无法使用边界条件；
csape只是Cubic spline插值，interp1可以选择几种不同的插值方法。

注：本文interp1函数和代码使用部分引用于matlab自带的插值函数interp1的几种插值方法 。个人在原文上修改了一小部分且补充了一点csape函数使用方面，自以为豪得是本文注释讲解做的格外用心.

还有就是上面的案例都是一维插值(插值函数为一元函数)，二维插值(插值函数为二元函数-曲面)使用interp2函数或仍使用csape函数，涉及到具体问题得话可以自寻查找相关内容。

参考文献：《数学建模算法与应用》.