建模算法（八）——插值

忘是亡心i 2021-12-14 07:05 380阅读 0赞

插值：求过已知有限个数据点的近似函数

拟合：已知有限个数据点，求近似函数，不要求过已知数据点，只要求在某种意义下在这些点的误差最小

# （一）插值方法 #

## 一、拉格朗日多项式插值 ##

### 1、插值多项式 ###

就是做出一个多项式函数，经过给出的n个节点，并尽可能的接近原函数，将点带入多项式函数得到一个线性方程组

[![image][]][image 1]

当系数矩阵满秩时，有唯一解。而，系数矩阵的行列式为

[![image][image 2]][image_image 2]

这是一个范德蒙德行列式，只要各个节点不同时，行列式就不为0，因此可得，一定能够解出系数方程

还有些指标

[![image][image 3]][image_image 3]

[![image][image 4]][image_image 4]

### 2、拉格朗日插值多项式 ###

[![image][image 5]][image_image 5]

### 3、MATLAB实现 ###

function y=lagrange(x0,y0,x)

%n个数据以数组x0,y0输入，m个插值点以数组x输入，输出数组y为m个插值。
    n=length(x0);m=length(x);
    for i=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

## 二、牛顿插值 ##

### 1、差商 ###

[![image][image 6]][image_image 6]

### 2、牛顿插值公式 ###

[![image][image 7]][image_image 7]

有点就是，多一个数据点，只多一项

PS：

[![image][image 8]][image_image 8]

### 3、差分 ###

[![image][image 9]][image_image 9]

[![image][image 10]][image_image 10]

### 4、等距节点插值公式 ###

[![image][image 11]][image_image 11]

## 三、分段线性插值 ##

### 1、插值多项式的振荡 ###

即如果插值多项式的次数越高，越容易发生振荡，不能很好的拟合。

### 2、分段线性插值 ###

[![image][image 12]][image_image 12]

[![image][image 13]][image_image 13]

### 3、MATLAB实现 ###

[![image][image 14]][image_image 14]

## 四、Hermite插值 ##

### 1、Hermite插值多项式 ###

[![image][image 15]][image_image 15]

[![image][image 16]][image_image 16]

### 2、MATLAB实现 ###

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

## 五、样条插值 ##

### 1、概念 ###

[![image][image 17]][image_image 17]

[![image][image 18]][image_image 18]

实际中最常用的是k=2和k=3的情况

## 二、二次样条函数插值 ##

### 二次样条函数有n+2个待定常数，所以要有n+2个条件，才能有唯一解 ###

[![image][image 19]][image_image 19]

一定要有一个条件为一阶导数

## 三、三次样本函数插值 ##

### 二次样条函数有n+3个待定常数，所以要有n+3个条件，才能有唯一解 ###

[![image][image 20]][image_image 20]

[![image][image 21]][image_image 21]

## 四、三次插值在MATLAB中的实现 ##

部分转载

**1、y=interp1(x0,y0,x,\`spline\`);     %** (spline改成linear，则变成线性插值)

**2、y=spline(x0,y0,xi);%**这个是根据己知的x，y数据，用样条函数插值出xi处的值。即由x,y的值计算出xi对应的函数值。

**3、pp=spline(x0,y0);%**是由根据己知的x，y数据，求出它的样条函数表达式，不过该表达式不是用矩阵直接表示，要求点x\`的值，要用函数**y\`=ppval(pp,x\`);**

**4、pp=csape(x,y,'变界类型','边界值conds');**生成各种边界条件的三次样条插值. 其中,(x,y)为数据向量,边界类型可为:

'complete'：给定边界一阶导数,即默认的边界条件,Lagrange边界条件  
             'not-a-knot'：非扭结条件,不用给边界值.  
             'periodic'：周期性边界条件,不用给边界值.  
             'second'：给定边界二阶导数.  
             'variational'：自然样条(边界二阶导数为\[0,0\]。

[![image][image 22]][image_image 22]

[![image][image 23]][image_image 23]

## 五、demo ##

[![image][image 24]][image_image 24]

%转载= =

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=lagrange(x0,y0,t);%调用编写的lagrange函数 
    dy1=(lagrange(x0,y0,0.0001)-lagrange(x0,y0,0))/0.0001%x=0处斜率 
    min1=min(lagrange(x0,y0,13:0.001:15))%13到15最小值 
    subplot(2,2,1); 
    plot(x0,y0,'ro',t,y1);%画出曲线 
    title('拉格朗日插值函数'); 
    %分段线性插值 
    y2=interp1(x0,y0,t,'spline');%注意区分spline与linear 
    Y2=interp1(x0,y0,t);%默认linear 
    dy2=(interp1(x0,y0,0.0001,'spline')-interp1(x0,y0,0,'spline'))/0.0001%x=0处斜率 
    min2=min(interp1(x0,y0,13:0.001:15,'spline'))%13到15最小值 
    subplot(2,2,2); 
    plot(t,y2,'b',t,Y2,'r',x0,y0,'ro');%画出曲线 
    title('分段线性插值'); 
    legend('边条','线性');%显示图形图例 
    %三次线条插值A 
    y3=spline(x0,y0,t); 
    dy3=(spline(x0,y0,0.0001)-spline(x0,y0,0))/0.0001%x=0处斜率 
    min3=min(spline(x0,y0,13:0.001:15))%13到15最小值 
    subplot(2,2,3); 
    plot(x0,y0,'ro',t,y3);%画出曲线 
    title('三次线条插值A'); 
    %三次线条插值B 
    pp1=csape(x0,y0);%默认的边界条件,即给定边界一阶导数 
    pp2=csape(x0,y0,'second');%给定边界二阶导数 
    y4=ppval(pp1,t); 
    Y4=ppval(pp2,t); 
    dy4=(ppval(pp1,0.0001)-ppval(pp1,0))/0.0001%x=0处斜率 
    min4=min(ppval(pp1,13:0.001:15))%13到15最小值 
    subplot(2,2,4); 
    plot(t,y4,'b',t,Y4,'r',x0,y0,'ro');%画出曲线 
    title('三次线条插值B'); 
    legend('一阶','二阶');

[![image][image 25]][image_image 25]

## 七、二维插值 ##

如果节点是二维的，插值函数是二元函数的话（曲面），我们可以画出三维的效果图

1、插值节点为网络节点

[![image][image 26]][image_image 26]

MatLab封装程序

z=interp2(x0,y0,z0,x,y,'method')

a、x0,y0为节点坐标，z0为n\*m维矩阵，表示节点的值

b、x0,y0要求一个为行向量一个为列向量

c、z为矩阵，n=length(y)，m=length(x)       因为MATLAB是列优先

d、[![image][image 27]][image_image 27]

e、x,y为插值点坐标，z为函数值

然后如果是三次样条插值，可以使用命令

pp=csape({x0,y0},z0,conds,valconds),z=fnval(pp,{x,y})

a、x0,y0为节点坐标，z0为n\*m维矩阵，表示节点的值

b、x0,y0要求一个为行向量一个为列向量

c、“conds”与一维相同，一般默认

d、x,y为插值点坐标 ，z为函数值

2、demo

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');  %注意z矩阵的行列所对应的向量
    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');

[![image][image 28]][image_image 28]

## 二、插值节点为散乱节点 ##

1、定义

[![image][image 29]][image_image 29]

MATLAB提供了一个函数

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

a、x,y,z为n维向量，就是数据点

b、xi,yi是给定的网格点横纵坐标（插值点），返回zi的值

2、demo

[![image][image 30]][image_image 30]

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: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');

[![image][image 31]][image_image 31]

转载于:https://www.cnblogs.com/BlueMountain-HaggenDazs/p/4271977.html

[image]: /images/20211214/e16ddc32df7949b2b834a8527d1ef4a6.png
[image 1]: https://images0.cnblogs.com/blog/657755/201502/041146073743760.png
[image 2]: https://images0.cnblogs.com/blog/657755/201502/041146093745104.png
[image_image 2]: https://images0.cnblogs.com/blog/657755/201502/041146087964760.png
[image 3]: https://images0.cnblogs.com/blog/657755/201502/041146105936576.png
[image_image 3]: https://images0.cnblogs.com/blog/657755/201502/041146098591718.png
[image 4]: https://images0.cnblogs.com/blog/657755/201502/041146122817890.png
[image_image 4]: https://images0.cnblogs.com/blog/657755/201502/041146117654034.png
[image 5]: https://images0.cnblogs.com/blog/657755/201502/041146143435721.png
[image_image 5]: https://images0.cnblogs.com/blog/657755/201502/041146132349048.png
[image 6]: https://images0.cnblogs.com/blog/657755/201502/041146155461423.png
[image_image 6]: https://images0.cnblogs.com/blog/657755/201502/041146149064592.png
[image 7]: https://images0.cnblogs.com/blog/657755/201502/041146172342737.png
[image_image 7]: https://images0.cnblogs.com/blog/657755/201502/041146166097379.png
[image 8]: https://images0.cnblogs.com/blog/657755/201502/041146185316466.png
[image_image 8]: https://images0.cnblogs.com/blog/657755/201502/041146179213582.png
[image 9]: https://images0.cnblogs.com/blog/657755/201502/041146201566996.png
[image_image 9]: https://images0.cnblogs.com/blog/657755/201502/041146191096810.png
[image 10]: https://images0.cnblogs.com/blog/657755/201502/041146215465455.png
[image_image 10]: https://images0.cnblogs.com/blog/657755/201502/041146208129597.png
[image 11]: https://images0.cnblogs.com/blog/657755/201502/041146230151471.png
[image_image 11]: https://images0.cnblogs.com/blog/657755/201502/041146225461628.png
[image 12]: https://images0.cnblogs.com/blog/657755/201502/041146240938900.png
[image_image 12]: https://images0.cnblogs.com/blog/657755/201502/041146235785043.png
[image 13]: https://images0.cnblogs.com/blog/657755/201502/041146254994130.png
[image_image 13]: https://images0.cnblogs.com/blog/657755/201502/041146247658271.png
[image 14]: https://images0.cnblogs.com/blog/657755/201502/041146266246574.png
[image_image 14]: https://images0.cnblogs.com/blog/657755/201502/041146260314758.png
[image 15]: https://images0.cnblogs.com/blog/657755/201502/041146283436132.png
[image_image 15]: https://images0.cnblogs.com/blog/657755/201502/041146273596730.png
[image 16]: https://images0.cnblogs.com/blog/657755/201502/041146293906318.png
[image_image 16]: https://images0.cnblogs.com/blog/657755/201502/041146289998732.png
[image 17]: https://images0.cnblogs.com/blog/657755/201502/041146301872962.png
[image_image 17]: https://images0.cnblogs.com/blog/657755/201502/041146297813903.png
[image 18]: https://images0.cnblogs.com/blog/657755/201502/041146311247348.png
[image_image 18]: https://images0.cnblogs.com/blog/657755/201502/041146306871048.png
[image 19]: https://images0.cnblogs.com/blog/657755/201502/041146321092049.png
[image_image 19]: https://images0.cnblogs.com/blog/657755/201502/041146315781421.png
[image 20]: https://images0.cnblogs.com/blog/657755/201502/041146330158192.png
[image_image 20]: https://images0.cnblogs.com/blog/657755/201502/041146325468349.png
[image 21]: https://images0.cnblogs.com/blog/657755/201502/041146339533580.png
[image_image 21]: https://images0.cnblogs.com/blog/657755/201502/041146334848035.png
[image 22]: https://images0.cnblogs.com/blog/657755/201502/041146355157622.png
[image_image 22]: https://images0.cnblogs.com/blog/657755/201502/041146345469693.png
[image 23]: https://images0.cnblogs.com/blog/657755/201502/041146368439595.png
[image_image 23]: https://images0.cnblogs.com/blog/657755/201502/041146362038467.png
[image 24]: https://images0.cnblogs.com/blog/657755/201502/041146387341439.png
[image_image 24]: https://images0.cnblogs.com/blog/657755/201502/041146377496739.png
[image 25]: https://images0.cnblogs.com/blog/657755/201502/041146398593883.png
[image_image 25]: https://images0.cnblogs.com/blog/657755/201502/041146393747269.png
[image 26]: https://images0.cnblogs.com/blog/657755/201502/041227490783978.png
[image_image 26]: https://images0.cnblogs.com/blog/657755/201502/041227487188933.png
[image 27]: https://images0.cnblogs.com/blog/657755/201502/041227500461907.png
[image_image 27]: https://images0.cnblogs.com/blog/657755/201502/041227496402849.png
[image 28]: https://images0.cnblogs.com/blog/657755/201502/041227515156921.png
[image_image 28]: https://images0.cnblogs.com/blog/657755/201502/041227507034507.png
[image 29]: https://images0.cnblogs.com/blog/657755/201502/041227525935351.png
[image_image 29]: https://images0.cnblogs.com/blog/657755/201502/041227521245508.png
[image 30]: https://images0.cnblogs.com/blog/657755/201502/041227537496038.png
[image_image 30]: https://images0.cnblogs.com/blog/657755/201502/041227530933437.png
[image 31]: https://images0.cnblogs.com/blog/657755/201502/041227548284468.png
[image_image 31]: https://images0.cnblogs.com/blog/657755/201502/041227543907166.png