Matlab实现插值与拟合

最新推荐文章于 2024-08-31 10:02:07 发布

键盘里的青春

最新推荐文章于 2024-08-31 10:02:07 发布

阅读量7.9k

点赞数 10

转自：点击打开链接

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（或y'n）

（2）已知插值节点xi 和相应的导数值 y'i (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) = y'n 。由这种边界条件建立的样条插值函数称为 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) = N'a (a), s' (b) = N'b (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
x 0 3 5 7 9 11 12 13 14 15
y 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('一阶','二阶');

运行结果：

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

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

最后，将上述所有情况封装起来，变成下列函数：

function y = showAllInterp( x0,y0,s,t)
%显示x0,y0之间所有不同类型的插值情况
%字符串s选择要输出的插值类型：
%all：全部类型           lagrange：拉格朗日插值函数
%linear：分段线性插值    spline：三次线条插值A
%csape：三次线条插值B
if(nargin<4)
    t=linspace(x0(1),x0(length(x0)),500);%默认
end
switch s
    case 'lagrange' %拉格朗日插值函数
        y1=lagrange(x0,y0,t);%调用编写的lagrange函数
        plot(x0,y0,'ro',t,y1);%画出曲线
        title('拉格朗日插值函数');
        if(nargout==1)
            y=y1;
        end
    case 'linear' %分段线性插值
        y2=interp1(x0,y0,t);%默认linear 
        plot(x0,y0,'ro',t,y2,'b');%画出曲线
        title('分段线性插值');
        if(nargout==1)
            y=y2;
        end
    case 'spline' %三次线条插值A
        y3=spline(x0,y0,t); %等价于interp1(x0,y0,t,'spline');
        plot(x0,y0,'ro',t,y3);%画出曲线 
        title('三次线条插值A'); 
        if(nargout==1)
            y=y3;
        end
    case 'csape' %三次线条插值B
        pp1=csape(x0,y0);%默认的边界条件,即给定边界一阶导数 
        pp2=csape(x0,y0,'second');%给定边界二阶导数 
        y4=ppval(pp1,t); 
        Y4=ppval(pp2,t); 
        plot(t,y4,'b',t,Y4,'r',x0,y0,'ro');%画出曲线
        title('三次线条插值B');
        legend('一阶','二阶');
        if(nargout==1)
            y=y4;
        end
    case 'all' %显示全部
        y1=lagrange(x0,y0,t);%调用编写的lagrange函数
        subplot(2,2,1);
        plot(x0,y0,'ro',t,y1);%画出曲线
        title('拉格朗日插值函数');
        
        y2=interp1(x0,y0,t);%默认linear 
        subplot(2,2,2);
        plot(x0,y0,'ro',t,y2);%画出曲线
        title('分段线性插值');
            
        y3=spline(x0,y0,t); %等价于interp1(x0,y0,t,'spline');
        subplot(2,2,3);
        plot(x0,y0,'ro',t,y3);%画出曲线 
        title('三次线条插值A');

        pp1=csape(x0,y0);%默认的边界条件,即给定边界一阶导数 
        pp2=csape(x0,y0,'second');%给定边界二阶导数 
        y4=ppval(pp1,t); 
        Y4=ppval(pp2,t);
        subplot(2,2,4);
        plot(t,y4,'b',t,Y4,'r',x0,y0,'ro');%画出曲线
        title('三次线条插值B');
        legend('一阶','二阶');
end

二维插值之插值节点为网格节点

已知m x n个节点：（xi，yj，zij）(i=1…m,j=1…n)，且xi,yi递增。求（x，y）处的插值z。

Matlab可以直接调用interp2(x0,y0,z0,x,y,`method`)

其中 x0,y0 分别为 m 维和 n 维向量，表示节点， z0 为 n × m 维矩阵，表示节点值， x,y
为一维数组，表示插值点， x 与 y 应是方向不同的向量，即一个是行向量，另一个是列
向量， z 为矩阵，它的行数为 y 的维数，列数为 x 的维数，表示得到的插值， 'method'
的用法同上面的一维插值。

如果是三次样条插值，可以使用命令
pp=csape({x0,y0},z0,conds,valconds), z=fnval(pp,{x,y})
其中 x0,y0 分别为 m 维和 n 维向量， z0 为 m × n 维矩阵， z 为矩阵，它的行数为 x 的维
数，列数为 y 的维数，表示得到的插值，具体使用方法同一维插值。

eg:

(1)、用interp2函数插值：

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]; 
p=100:1:500; 
q=100:1:400; 
q=q';%须为列向量 
z0=interp2(x,y,z,p,q);%分段线性插值 
z1=interp2(x,y,z,p,q,'spline');%三次线条插值 
subplot(2,1,1); 
mesh(p,q,z0); 
title('分段线性插值'); 
subplot(2,1,2); 
mesh(p,q,z1); 
title('三次线条插值'); 
%可以观察出，三次线条插值的图像更平滑

运行结果：

(2)、用csape函数插值：

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]; 
p=100:1:500; 
q=100:1:400; 
q=q'; 
%三次线条插值 
pp=csape({x,y},z');%注意跟interp2的区别,有个转置 
z0=fnval(pp,{p,q}); 
mesh(p,q,z0');%注意跟interp2的区别,有个转置 
title('三次线条插值');

运行结果：

2、二维插值之插值节点为散乱节点

已知 n 个节点： ( xi , yi , zi )(i = 1,2,…, n) ，求点 (x, y) 处的插值 z 。
对上述问题， Matlab 中提供了插值函数 griddata，其格式为：
ZI = GRIDDATA(X,Y,Z,XI,YI)
其中 X、 Y、 Z 均为 n 维向量，指明所给数据点的横坐标、纵坐标和竖坐标。向量 XI、
YI 是给定的网格点的横坐标和纵坐标，返回值 ZI 为网格（ XI， YI）处的函数值。 XI
与 YI 应是方向不同的向量，即一个是行向量，另一个是列向量。

eg:

%散乱节点的二维插值 
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]; 
x0=[75:1:200]; 
y0=[-85:1:145]'; 
z0=griddata(x,y,z,x0,y0,'cubic');%保凹凸性3次插值 
%[xx,yy]=meshgrid(x0,y0);无需采样，故不需要该函数 
mesh(x0,y0,z0);

运行结果:

在上述问题中,补上寻找最大值的程序:

%max(z0)返回一个行向量，向量的第i个元素是矩阵A的第i列上的最大值 
%find(A) 寻找矩阵A非零元素下标,返回矩阵A中非零元素所在位置 
%[i,j,v]=find(A)返回矩阵A中非零元素所在的行i,列j,和元素的值v(按所在位置先后顺序输出) 
[p,q]=find(z0==max(max(z0))); 
zmax=z0(p,q)

3、最小二乘法实现曲线拟合

（1）用最小二乘法求一个形如 y = a + bx^ 2 的经验公式：

%等价于[1,x^2][a;b]=y,转换成解超定方程问题，超定方程的解是根据最小二乘法得来的

x=[19 25 31 38 44]'; 
y=[19.0 32.3 49.0 73.3 97.8]'; 
r=[ones(5,1),x.^2] 
ab=r\y 
x0=19:0.1:44; 
y0=ab(1)+ab(2)*x0.^2; 
plot(x,y,'o',x0,y0,'r')

（2）多项式拟合

%a=polyfit(x,y,n)用多项式求过已知点的表达式，其中x为源数据点对应的横坐标，可为行向量、矩阵，y为源数据点对应的纵坐标，可为行向量、矩阵，n为你要拟合的阶数，一阶直线拟合，二阶抛物线拟合，并非阶次越高越好，看拟合情况而定，a为m+1的行向量。polyfit函数的数学基础是最小二乘法曲线拟合原理，所得到的函数值在基点处的值与原来点的坐标偏差最小，常用于数据拟合，polyfit 做出来的值从左到右表示从高次到低次的多项式系数。

如果要求拟合函数在x`点的函数值，可以调用polyval(a,x`)函数

eg：

x0=[1990 1991 1992 1993 1994 1995 1996]; 
y0=[70 122 144 152 174 196 202]; 
%画出散点图 
plot(x0,y0,'ro'); 
hold on 
%用线性拟合 
p=polyfit(x0,y0,1); 
z0=polyval(p,x0); 
plot(x0,z0);

运行结果：

4、最小二乘优化 (最小二乘:least square)

1、lsqlin函数

eg：

%拟合形如y=a+bx^2的函数 
%采样点 
x=[19 25 31 38 44]'; 
y=[19 32.3 49 73.3 97.8]'; 
r=[ones(5,1),x.^2]; 
ab=lsqlin(r,y) 
x0=19:0.1:44; 
y0=ab(1)+ab(2)*x0.^2; 
plot(x,y,'o',x0,y0,'r')

运行结果：

5、曲线拟合与函数逼近

eg：

求 f(x) =cos x, （-pi/2<=x<=pi/2）在 H = Span{1, x^2 , x^4} 中的最佳平方逼近多项式。

程序如下：

syms x%定义符号数值 
base=[1,x^2,x^4]; 
y1=base.'*base 
y2=cos(x)*base.' 
r1=int(y1,-pi/2,pi/2) 
r2=int(y2,-pi/2,pi/2) 
a=r1\r2%a为符号数值 
xishu1=double(a)%化简符号数值 
digits(8)%设置符号数值的精度 
xishu2=vpa(a)%任意精度（符号类）数值

运行结果：

 
         y1 =  
        
         [   1, x^2, x^4]  
        
         [ x^2, x^4, x^6]  
        
         [ x^4, x^6, x^8]  
        
         y2 =  
        
              cos(x)  
        
         x^2*cos(x)  
        
         x^4*cos(x)  
        
         r1 =  
        
         [      pi,  pi^3/12,   pi^5/80]  
        
         [ pi^3/12,  pi^5/80,  pi^7/448]  
        
         [ pi^5/80, pi^7/448, pi^9/2304]  
        
         r2 =  
        
                             2  
        
                    pi^2/2 - 4  
        
         pi^4/8 - 6*pi^2 + 48  
        
         a =  
        
         (15*(pi^4 - 308*pi^2 + 3024))/(4*pi^5)  
        
            -(210*(pi^4 - 228*pi^2 + 2160))/pi^7  
        
            (1260*(pi^4 - 180*pi^2 + 1680))/pi^9  
        
         xishu1 = 
        
             0.9996  
        
            -0.4964  
        
             0.0372 
        
         xishu2 =  
        
           0.99957952  
        
         -0.49639233  
        
         0.037209327  
        
         >>