Matlab随笔之插值与拟合(上)

1、拉格朗日插值

新建如下函数:

function y=lagrange(x0,y0,x)
%拉格朗日插值函数
%n 个节点数据以数组 x0, y0 输入(注意 Matlat 的数组下标从1开始),
%m 个插值点以数组 x 输入,输出数组 y 为 m 个插值
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

应用实例:

x0=1:1:20;
y0=x0.^2-20*x0-5;
x=1:0.1:20;
z=lagrange(x0,y0,x);
plot(x,z,‘:‘,x0,y0,‘ko‘);

运行结果:

2、分段线性插值

MATLAB现成的插值函数为interp1,其调用格式为:  yi= interp1(x,y,xi,‘method‘)

其中x,y为插值点,yi为在被插值点xi处的插值结果;x,y为向量, ‘method‘表示采用的插值方法,包括:
‘method‘:是最近项插值;                                                              ‘linear‘:线性插值;(默认)

‘spline‘:逐段3次样条插值; (下面的三次样条插值会用到)             ‘cubic‘:保凹凸性3次插值

‘pchip‘:分段三次Hermite 插值。

例如:在一天24小时内,从零点开始每间隔2小时测得的环境温度数据分别为
            12,9,9,1,0,18 ,24,28,27,25,20,18,15,13,
推测中午12点(即13点)时的温度.

x=0:2:24;
y=[12,9,9,10,18,24,28,27,25,20,18,15,13];
x1=0:0.5:24;
y1=interp1(x,y,x1,‘linear‘);
plot(x,y,‘bo‘,x1,y1,‘r:‘);

运行结果:

3、埃尔米特插值

如果要求插值函数不仅在节点处与函数同值,而且要求它与函数有相同的一阶、二阶甚至高阶导数值,这就是埃尔米特插值问题。

已知f(x)的n+1个节点的函数值f(xi)以及导数值f`(xi),可得一个至多n+1次的多项式H(x),即hermite插值多项式。

新建以下这个函数:

function y = hermite( x0,y0,y1,x )
%埃尔米特插值多项式
%x0为点横坐标
%y0为函数值
%y1为导数值
%m个插值点用数组x输入
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=1: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)-x(k))*(2*a*y0(i)-y1(i))+y0(i));
end
y(k)=yy;
end

4、样条插值

所谓样条( Spline)本来是工程设计中使用的一种绘图工具,它是富有弹性的细木条或细金属条。绘图员利用它把一些已知点连接成一条光滑曲线(称为样条曲线),并使连接点处有连续的曲率。

数学上将具有一定光滑性的分段多项式称为样条函数。

在实际中最常用的是二次样条函数和三次样条函数:

二次样条函数插值


首先,我们注意到s2 (x)中含有 n + 2 个特定常数,故应需要 n + 2 个插值条件,因此,二次样条插值问题可分为两类:

(1)已知插值节点xi 和相应的函数值 yi (i = 0,1,…,n) 以及端点 x0 (或 xn )处的导数值y‘0(或yn)

(2)已知插值节点xi 和相应的导数值 yi (i = 0,1,…,n) 以及端点 x0 (或 xn )处的函数值y0 (或yn

三次样条函数插值

由于 s3 (x)中含有n + 3 个待定系数,故应需要 n + 3 个插值条件,已知插值节点xi 和相应的函数值 f(xi ) = yi (i = 0,1,…,n) ,这里提供了 n + 1 个条件,还需要 2 个边界条件。因此,三次样条插值问题可分为三类:

(1)s‘3 (a) = y‘0 ,s‘3 (b) = yn 。由这种边界条件建立的样条插值函数称为 f(x) 的完备三次样条插值函数。
特别地,y0‘ = yn`= 0时,样条曲线在端点处呈水平状态。
如果 f‘ (x) 不知道,我们可以要求 s‘3 (x) 与 f‘ (x) 在端点处近似相等。这时以x0 , x1 , x2 , x3 为节点作一个三次 Newton 插值多项式 Na (x) ,以 xn, xn−1, xn−2, xn−3 作一个三次 Newton 插值多项式 Nb (x) ,要求s‘ (a) = Na (a), s‘ (b) = Nb (b)
由这种边界条件建立的三次样条称为 f(x) 的 Lagrange 三次样条插值函数。

(2)s"3 (a) = y"0 ,s"3 (b) = y"3 。特别地 y"n = y"n = 0 时,称为自然边界条件。

(3)s‘3 ( a + 0) = s‘3 ( b − 0), s"3 (a + 0) = s"3 (b − 0) , (这里要求 s3 (a + 0) =s3 (b − 0) )此条件称为周期条件。

Matlab实现(三次样条插值)

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]

边界值conds可用1x2矩阵表示,矩阵元素取值为1,2,此时,使用命令
pp=csape(x0,y0_ext,conds)
其中 y0_ext=[left, y0, right],这里 left 表示左边界的取值, right 表示右边界的取值。
conds(i)=j 的含义是给定端点 i j 阶导数, 即 conds 的第一个元素表示左边界的条
件,第二个元素表示右边界的条件, conds=[2,1]表示左边界是二阶导数,右边界是一阶
导数,对应的值由 left 和 right 给出。

例子:

表 1
0  3      5   7     9  11  12  13  14  15
0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6
要求用 Lagrange、分段线性和三次样条三种插值方法计算。

编程实现:

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(‘一阶‘,‘二阶‘);

运行结果:

dy1 =

  -55.2855

min1 =

    0.9391

dy2 =

    0.5023

min2 =

    0.9828

dy3 =

    0.5023

min3 =

    0.9828

dy4 =

    0.5007

min4 =

    0.9851

综上,可以看出,拉格朗日插值函数根本不能应用,分段线性函数的光滑性较差,推荐三次样条插值。

同时,可以看出,interp1(x0,y0,’spline’)等价于spline(x0,y0)。

时间: 2024-08-01 10:25:14

Matlab随笔之插值与拟合(上)的相关文章

Matlab随笔之线性规划

原文:Matlab随笔之线性规划   LP(Linear programming,线性规划)是一种优化方法,在优化问题中目标函数和约束函数均为向量变量的线性函数,LP问题可描述为:min xs.t. A·x b Aeq·x=beq vlb x vub其中 ,b,beq均为向量,A,Aeq为矩阵,x为向量变量.矩阵A和向量b是线性不等式约束条件的系数, Aeq和beq是等式约束条件的系数. 在MATLAB中,用于LP的求解函数为linprog.其调用格式为:[x,fval,lambda]=linp

插值与拟合

1.插值 -->求过已知有限个数据点的近似函数 1)拉格朗日多项式插值 -->n个插值点不同时确定了一个唯一的n次多项式 构造n次拉格朗日插值多项式(不使用解方程n个约束来求解待定系数) 2)牛顿插值 使用差商概念来构造牛顿插值公式(计算量小,余项与拉格朗日余项相等),当节点之差为常数时,使用差分来代替差商构造牛顿向前插值公式 3)分段线性插值 -->高次插值存在震荡缺陷,采用低次分段函数(线性函数) y=interp1(x0,y0,x,'method') -->method可取n

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

插值:求过已知有限个数据点的近似函数 拟合:已知有限个数据点,求近似函数,不要求过已知数据点,只要求在某种意义下在这些点的误差最小 (一)插值方法 一.拉格朗日多项式插值 1.插值多项式 就是做出一个多项式函数,经过给出的n个节点,并尽可能的接近原函数,将点带入多项式函数得到一个线性方程组 当系数矩阵满秩时,有唯一解.而,系数矩阵的行列式为 这是一个范德蒙德行列式,只要各个节点不同时,行列式就不为0,因此可得,一定能够解出系数方程 还有些指标 2.拉格朗日插值多项式 3.MATLAB实现 fun

Matlab随笔之模拟退火算法

问题描述: 我方有一个基地,经度和纬度为( 70,40).假设我方飞机的速度为 1000 公里/小时. 我方派一架飞机从基地出发,侦察完敌方所有目标,再返回原来的基地.在敌方每一目 标点的侦察时间不计,求该架飞机所花费的时间(假设我方飞机巡航时间可以充分长). 这是一个旅行商问题.我们依次给基地编号为 1,敌方目标依次编号为 2, 3,…, 101, 最后我方基地再重复编号为 102(这样便于程序中计算). 距离矩阵 D = ( dij )102×102 , 其中 dij 表示表示 i, j 两

Matlab随笔(1)之画图函数总结

MATLAB函数画图 MATLAB不但擅长於矩阵相关的数值运算,也适合用在各种科学目视表示(Scientific visualization).本节将介绍MATLAB基本xy平面及xyz空间的各项绘图命令,包含一维曲线及二维曲面的绘制.列印及存档. 基本XY平面绘图命令 Plot 是绘制一维曲线的基本函数,但在使用此函数之前,我们需先定义曲线上每一点的x及y座标.下例可画出一条正弦曲线: close all; %linspace(5,100,20)和5:5:100的区别在于前者已知元素总个数而不

MATLAB实现牛顿插值

牛顿插值法的原理,在维基百科上不太全面,具体可以参考这篇文章.同样贴出,楼主作为初学者认为好理解的代码. function p=Newton1(x1,y,x2) %p为多项式估计出的插值 syms x n = length(x1); %差商的求法 for i=2:n f1(i,1)=(y(i)-y(i-1))/(x1(i)-x1(i-1)); end for i=2:n for j=i+1:n f1(j,i)=(f1(j,i-1)-f1(j-1,i-1))/(x1(j)-x1(j-i)); en

Matlab实现线性回归(直线拟合)

线性拟合:对于y=a*x+b的形式 a=(N*ΣXi*Yi-ΣXi*ΣYi)/(N*ΣXi*Xi-(ΣXi)2) b=((ΣXi*Xi)*(ΣYi)-(ΣXi)*(ΣXi*Yi))/(N*ΣXi*Xi-(ΣXi)2) 利用Matlab自带函数可实现: 拟合函数:pn=polyfit(x,y,n)  返回pn系数向量,降阶排列,n为阶数 函数:yy=polyval(pn,x) pn为降阶排列的多项式系数,x为向量或者矩阵,返回yy:将x带入pn后计算得到的向量或者矩阵. %初始数据 x=[1 2 

matlab 给某一列乘上一个系数

矩阵M是一个 mxn 的矩阵,现在要给M矩阵的第一列都要乘上10,使其第一列扩大10倍,那肿么做呢? 我第一时间用的是: M(:,1) = M(:,1)*10; //错误的 但是这个错了,结果是不对的 这里要用点乘才行,所以正确的写法是: M(:,1) = M(:,1).*10;  // 正确写法 这里要明白矩阵乘积与矩阵点乘(对应位相乘)的区别 矩阵乘法要求左矩阵的行数与右矩阵的列数相等,即MxN维矩阵乘以NxM维矩阵 例:A=[1 1 1;2 2 2]  B=[3 3;4 4;5 5] MA

matlab安装小坑----连接不上服务器

安装正版matlab时显示连接不上服务器,在贴吧看到有人说用阿里DNS.附上网址 http://www.alidns.com/setup/ 里面有一个一件快速设置小软件挺好用(下载时会报病毒,直接忽略这条信息),可以键设置DNS,不好用可以用管理员身份运行,再不好用就按照网址上的步骤手动设置吧. 但是我设置完了还是不好用,拔了网线,大功告成. -----------------------------------------------------------------------------