MTLAB------2019/7/20

插值与拟合

1. MATLAB插值工具箱：、

• 一维插值（节点为一维变量，插值函数为一元函数（曲线））：
  （1） MATLAB中现有一维插值函数interpl,语法：
  y=interpl(x0,y0,x,‘method’)
  其中：method指定插值的办法，默认为线性插值。其值可为
  ‘nearest’   最近项插值
  ‘linear’     线性插值
  ‘spline’    立方样条插值
  ‘cubic’     立方插值
  
  注:所有插值方法要求x0是单调的
  （2）三次样条插值有如下函数：
  
  y=interpl(x0,y0,x,‘spline’);
  y=spline(x0,y0,x);
  pp=csape(x0,y0,conds);
  pp=csape(x0,y0,conds,valconds);y=fnva(pp，x）
  其中：x0,y0是已知数据点；x是插值点，y是插值点的函数值
  
  这几种函数的区别：
  1）三次样条曲线有4种边界条件。
  自然边界条件：二阶导数在边界处为0，可视为简支梁，是最常用的边界条件。
  第一边界条件：二阶导在边界处已知的边界条件。自然边界条件可视为特例。
  第二边界条件：一阶导在边界处已知的边界条件。
  循环边界条件：一阶导与二阶导在边界处相等的边界条件，适用于封闭或循环的图形。
  （原文：https://blog.csdn.net/weixin_42943114/article/details/84843632 ）

  2）spline函数只能实现自然边界条件和第二边界条件；
  csape函数可以实现三次样条曲线的各种条件，包括第一边界、第二边界、循环边界、混合边界等。可以实现一维至多维的各种情况。
p88页例5.1

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,'linear');
y2=interp1(x0,y0,x,'spline');
pp1=csape(x0,y0);
y3=fnval(pp1,x);
pp2=csape(x0,y0,'second')  %'second'边界为二阶导数，
                           %因在pp=csape(x0,y0,conds,valconds)中忽略二阶导数的值：
                           %valconde故二阶导数的值默认为[0.0];
                           
y4=fnval(pp2,x);
[x',y1',y2',y3',y4'];
subplot(1,3,1);
plot(x0,y0,'+',x,y1);
title('piecewisw linear')
subplot(1,3,2);
plot(x0,y0,'+',x,y2);
title('spline1');
subplot(1,3,3);
plot(x0,y0,'+',x,y3);
title('spline2');
dx=diff(x);                        %x的导数
dy=diff(y3);                       %y3的导数
dy_dx=dy./dx;                      %曲线斜率
dy_dx0=dy_dx(1)              %当x=0（因为下标从1开始，所以写dy_dx（1））时的曲线斜率
ytemp=y3(13:15);              %x13<=x<=15的范围内y的值
ymin=min(ytemp);                  %x13<=x<=15的范围内y的最小值
index=find(y3==ymin);           %找出满足y3=ymin的x的位置
xmin=x(index) ;                        %找出满足y3=ymin的x的位置所对应的x的值
[xmin,ymin]

例5.2（利用三次样条插值法求积分函数）

x0=[0.15,0.16,0.17,0.18];
y0=[3.5,1.5,2.5,2.8];
pp=csape(x0,y0)
xishu=pp.coefs                        %显示每个区间上三次多项式的系数
s=quadl(@(t)ppval(pp,t),0.15,0.18)  % ppval函数: 给出三次样条插值pp在t处对应的函数值;
                                    % quadl函数用于求积分

• 二维插值（插值函数为二元函数即曲面，节点为二维）
  (1)插值节点为网状节点
  已知$m\times n$个节点：$（x_i,y_j,z_{ij}）$,求点$(x,y)$处的插值$z$
  MATLAB命令：
  z=interp2(x0,y0,z0,x,y,‘method’)
  其中：x,y为一维数组，x为行向量，y为列向量，z0为$n\times m$矩阵
  pp=csape({x0,y0},z0,conds,valconds), z=fnval(pp,{x,y});
  
  例5.3

x0=100:100:500;
y0=100:100:400;
z0=[636 697 624 478 450
    698 712 630 478 420
    680 674 598 412 400
    662 626 552 334 310];
z0=z0'                    %将z0转置成n×m矩阵
pp=csape({x0,y0},z0);
xi=100:100:500;
yi=100:100:400;
z=fnval(pp,{xi,yi});
[i,j]=find(z==max(max(z))); %(max(z))是找出z矩阵每一列中的最大值；max（max（z））为找出z矩阵中的最大值，find函数为找出z矩阵中最大值的位置
x=xi(i);
y=yi(j);
zmax=z(i,j);

例5.4

x0=[129 140 103.5 88 185.5 195 105 157.5 107.5 77 81 162 162 117.5];
y0=[7.5 141.5 23 147 22.5 137.5 85.5 -6.5 -81 3 56.5 -66.5 84 -33.5];
z0=-[4 8 6 8 6 8 8 9 9 8 8 9 4 9];
xmm=minmax(x0);    %求x0中的最大，最小值
ymm=minmax(y0);    %求y0中的最大，最小值
xi=xmm(1):xmm(2);  
yi=ymm(1):ymm(2);
yi=yi';
zi1=griddata(x0,y0,z0,xi,yi,'cubic');  %立方插值
zi2=griddata(x0,y0,z0,xi,yi,'nearest');%最近点插值
zi=zi1;   %立方插值与最近点插值的混合插值的初始值
zi(isnan(zi1))=zi2(isnan(zi1))；%把立方值中的不确定值换位最近点插值的结果；isnan函数判断组中的元素是否是无穷大，若是返回1.
subplot(1,2,1),plot(x0,y0,'*');
subplot(1,2,2),mesh(xi,yi,zi);%mesh函数画三维图


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