zoukankan      html  css  js  c++  java
  • Matlab随笔之插值与拟合(下)

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

    已知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:

    image

    (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('三次线条插值'); 
    %可以观察出,三次线条插值的图像更平滑 

    运行结果:

    image

    (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('三次线条插值');

    运行结果:

    image

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

    运行结果:

    image

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

    %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=ry 
    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);

    运行结果:

    image

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

    image

    1、lsqlin函数

    image

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

    运行结果:

    image

    5、曲线拟合与函数逼近

    image

    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
    2%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 
      
    >>
    

      

    所以y的最佳平方逼近多项式为y=0.9996-0.4964x^2+0.0372x^4

  • 相关阅读:
    配置navigation bar外观
    加急审核
    UIButton中setTitleEdgeInsets和setImageEdgeInsets的使用
    打开某个好友的聊天界面
    ALAssets的两种用法
    更改appstore开发商名字
    回调函数的使用
    相应字体所占的位置大小
    IOS的变量前加extern和static字段
    iOS开发代码规范(通用)
  • 原文地址:https://www.cnblogs.com/vanker/p/4268292.html
Copyright © 2011-2022 走看看