数学建模（8）——插值

插值与拟合均是根据一组数据构造一个函数作为近似，不同的是，插值要求求得的曲线通过所有的数据点，而拟合则是要求在某种意义下它在这些点上的总偏差最小

插值

1、一维插值函数

y = interp1（x0,y0,x,'method'）

mothod指定插值的方法，默认为线性插值

‘nearest’最近项插值

‘linear’线性插值

‘spline’立方样条插值

‘cubic’立方插值

2、三次样条插值

pp = csape（x0,y0,conds,valconds）;
y = fnval（pp,x）

‘conds’指定边界条件

‘complete’边界为一阶导数，一阶导数的值在valconds中给出，若忽略valconds，则使用Lagrange边界条件

‘not-a-knot’非扭结条件

‘periodic’周期条件

‘second’边界为二阶导数，二阶导数值在valconds中给出，若忽略valconds，默认值为【0，0】

‘variational’设置边界的二阶导数值为【0，0】

例题1

%此处使用四种插值方式作为比较
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];
x=0:0.1:15;
y1=lagrange(x0,y0,x);  %调用前面编写的Lagrange插值函数
y2=interp1(x0,y0,x);
y3=interp1(x0,y0,x,'spline');
pp1=csape(x0,y0);
y4=fnval(pp1,x);

subplot(2,2,1),
plot(x0,y0,'+',x,y1),
title('拉格朗日插值')
subplot(2,2,2),
plot(x0,y0,'+',x,y2),
title('分段线性插值')
subplot(2,2,3),
plot(x0,y0,'+',x,y3),
title('立方样条插值')
subplot(2,2,4),
plot(x0,y0,'+',x,y4),
title('三次样条插值')

dx=diff(x);
dy=diff(y4);
dy_dx=dy./dx;
dy_dx0=dy_dx(1) %x=0出的曲线斜率

ytemp=y4(131:151);
ymin=min(ytemp);
index=find(y4==ymin);
xmin=x(index);
[xmin,ymin] %13~15范围内y的最小值

function y=lagrange(x0,y0,x);
%拉格朗日插值
n=length(x0);m=length(x);
for i=1:m
z=x(i);
s=0.0;
for k=1:n
p=1.0;
for j=1:n
if j~=k
p=p*(z-x0(j))/(x0(k)-x0(j));
end
end
s=p*y0(k)+s;
end
y(i)=s;
end

clear,clc
x0=[0.15 0.16 0.17 0.18];
y0=[3.5 1.5 2.5 2.8];
pp=csape(x0,y0);
format long g
xishu = pp.coefs; %显示每个区间上三次多项式的系数
s = quadl(@(t)ppval(pp,t),0.15,0.18); %求积分,t为自变量
format %恢复端小数的显示格式

3、网格节点二维插值

pp = csape（{x0,y0},z0,conds,valconds）;
z = fnval（pp,{x,y}）

clear,clc
x=100:100:500;
y=100:100:400;
z=[636 697 624 478 450
698 712 630 478 420
680 674 598 412 400
662 626 552 334 310];
pp = csape({x,y},z');
xi = 100:10:500;yi=100:10:400;
cz=fnval(pp,{xi,yi});
[i,j]=find(cz==max(max(cz))) %找最高点的地址
x=xi(i),y=yi(j),zmax=cz(i,j) %求最高点坐标

4、散乱节点插值

zi = griddata( x,y,z,xi,yi )

clear,clc
x=[129  140  103.5  88  185.5  195  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:1:200;
yi=-50:1:150;
zi1=griddata(x,y,z,xi,yi','cubic') ;%立方插值
zi2=griddata(x,y,z,xi,yi','nearest');%最近点插值
zi=zi1;%立方插值与最近点插值的混合插值的初始值
zi(isnan(zi1))=zi2(isnan(zi1));%把立方插值中的不确定值换成最近点插值的结果
subplot(1,2,1), plot(x,y,'*')
subplot(1,2,2), mesh(xi,yi,zi)

Matlab插值时外插值是不确定的，这里使用混合插值，把不确定值换成最近点插值的结果