码迷,mamicode.com
首页 > 编程语言 > 详细

建模算法(八)——插值与拟合

时间:2015-02-04 12:49:18      阅读:457      评论:0      收藏:0      [点我收藏+]

标签:

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

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

(一)插值方法

一、拉格朗日多项式插值

1、插值多项式

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

技术分享

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

技术分享

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

        还有些指标

技术分享

技术分享

2、拉格朗日插值多项式

技术分享

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、差商

技术分享

2、牛顿插值公式

技术分享

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

PS:

技术分享

3、差分

技术分享

技术分享

4、等距节点插值公式

技术分享

三、分段线性插值

1、插值多项式的振荡

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

2、分段线性插值

技术分享

技术分享

3、MATLAB实现

技术分享

四、Hermite插值

1、Hermite插值多项式

技术分享

技术分享

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、概念

技术分享

技术分享

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

二、二次样条函数插值

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

技术分享

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

三、三次样本函数插值

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

技术分享

技术分享

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

技术分享

技术分享

五、demo

技术分享

%转载= =
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(一阶,二阶);

技术分享

七、二维插值

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

1、插值节点为网络节点

技术分享

       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、技术分享

     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‘);

技术分享

二、插值节点为散乱节点

1、定义

技术分享

    MATLAB提供了一个函数

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

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

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

2、demo

技术分享

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‘);
技术分享

建模算法(八)——插值与拟合

标签:

原文地址:http://www.cnblogs.com/BlueMountain-HaggenDazs/p/4271977.html

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!