zoukankan      html  css  js  c++  java
  • 建模算法(八)——插值

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

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

    (一)插值方法

    一、拉格朗日多项式插值

    1、插值多项式

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

    image

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

    image

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

            还有些指标

    image

    image

    2、拉格朗日插值多项式

    image

    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

    2、牛顿插值公式

    image

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

    PS:

    image

    3、差分

    image

    image

    4、等距节点插值公式

    image

    三、分段线性插值

    1、插值多项式的振荡

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

    2、分段线性插值

    image

    image

    3、MATLAB实现

    image

    四、Hermite插值

    1、Hermite插值多项式

    image

    image

    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

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

    二、二次样条函数插值

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

    image

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

    三、三次样本函数插值

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

    image

    image

    四、三次插值在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

    五、demo

    image

    %转载= =
    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

    七、二维插值

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

    1、插值节点为网络节点

    image

           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

         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

    二、插值节点为散乱节点

    1、定义

    image

        MATLAB提供了一个函数

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

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

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

    2、demo

    image

    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
  • 相关阅读:
    思维发散的双刃剑
    一个请求过来都经过了什么?(Thrift版)
    业务开发转基础开发,这三种「高可用」架构你会么?
    程序常用的设计技巧
    美团分布式服务通信框架及服务治理系统OCTO
    《程序员修炼之道》解读
    面试官说:你真的不是不优秀只是不合适
    架构视角-文件的通用存储原理
    那些影响深远的弯路
    iOS sqlite ORM框架-LKDBHelper
  • 原文地址:https://www.cnblogs.com/BlueMountain-HaggenDazs/p/4271977.html
Copyright © 2011-2022 走看看