Matlab实验四拟合

实验四.Matlab拟合一、实验目的直观了解拟合基本内容;掌握最小二乘拟合的基本原理；掌握用matlab命令求解拟合问题；掌握matlab拟合工具箱的基本使用方法二.实验原理与方法在科学计算中经常要建立实验数据的数学模型。给定函数的实验数据，需要用比较简单和合适的函数来逼近（或拟合）实验数据。这种逼近的特点是：(a)适度的精度是需要的；(b)实验数据有小的误差；(c)对于某些问题，可能有某些特殊的信息能够用来选择实验数据的数学模型。曲线拟合的数学描述描述：已知一组（二维）数据，即平面上n个点),...,2,1(),,(niyxii，ix互不相同。寻求一个函数（曲线）)(xfy使)(xf在某种准则下与所有数据点最为接近，即拟合的最好，不要求曲线（面）通过所有数据点，而是要求它反映对象整体的变化趋势。逼近离散数据的基本方法就是曲线拟合，常采用最小二乘拟合。最小二乘拟合分为线性最小二乘拟合与非线性最小二乘拟合1．线性最小二乘法（一）线性最小二乘法的基本思路：第一步:先选定一组基函数（线性无关）)(),(),(21xrxrxrm（mn）,构造函数)()()()(2211xraxraxraxfmm其中maaa,,,21为待定系数。第二步:确定系数maaa,...,,21，使n个点),(iiyx与曲线)(xf的距离i的平方和最小，该准则称为最小二乘准则，求系数maaa,...,,21的方法称为线性最小二乘拟合法。问题归结为，求a1,a2,…am使J(a1,a2,…am)最小。（1）记211211221])([])([),,(iiknimkkininiiimyxrayxfaaaJ则J为maaa,...,21的n元二次函数。由数学分析（或高等数学）的极值理论，J达到最小的充分必要条件是maaa,...,,21满足0...21maJaJaJ于是得到求使J达到最小的maaa,...,21的方法是求解线性方程组：mjyxraxrnimkiikkij,...,2,10])()[(11也即求解线性方程组：mjyxrxrxrainiijnkikijmkk,...,2,1)()]()([111若记mnnmnmxrxrxrxrR)(...)(.........)(...)(1111，TnTmyyyaaa],...,[,],...,[11则线性方程组可表示为：yRRaRTT由于当)(),(),(21xrxrxrm线性无关时，R列满秩，RRT可逆所以上述方程组有唯一解：yRRRaTT*)(1（称为该线性方程组的最小二乘意义下的解）。此即线性最小二乘拟合所得系数。用作线性最小二乘拟合的误差也可以考虑以下两种形式：1．最小平方误差：niiimxfyaaaJ1221)))(((),...,(2．最大偏差：)(),...,(max121iinimxfyaaaJ（二）基函数组的选取线性最小二乘拟合f(x)=a1r1(x)+…+amrm(x)中函数{r1(x),…rm(x)}的选取1.通过机理分析建立数学模型来确定f(x)，如果通过机理分析能知道y和x之间的关系，则12(),(),...,()mrxrxrx很容易可以确定。2.若无法知道y和x之间的关系，则可以通过对数据作图，直观的判断该组数据应用什么样的曲线去拟和。f=a1+a2xf=a1+a2x+a3x2常用曲线有如下几种：+++++++++++++++f=a1+a2/x+++++f=aebx（1）直线：21axay（2）多项式：1121...mmmmaxaxaxay（一般m=2,3,不宜太高）（3）双曲线（一支）：21axay，拟合前需做变量替换xt1转换成多项式，这样做比较简单。（4）指数函数：xaeay21，拟合前需做变量替换yzln，转换成xaaz21ln，这样做比较简单在实际计算过程中，面对一组数据，到底用什么曲线拟合好，可以在直观判断的基础上，选几种曲线分别拟合，然后作比较，选择拟合效果好的曲线。用MATLAB作线性最小二乘拟合1.作多项式f(x)=a1xm+…+amx+am+1拟合,也即选取基函数1121(),(),...,(),()1mmmmrxxrxxrxxrx可利用命令:a=polyfit(x,y,m)说明:yx,为原始样本点构成的向量，m为选定的多项式的阶数，a为拟合得到的多项式的系数，从高次项到低次项降幂排列，为长度为1m的数组。注：多项式在x处的值y可用以下命令计算：y=polyval（a，x）2．选取基函数：1121(),(),...,(),()1mmmmrxxrxxrxxrx建立超定方程组:yRRaRTT其中1............1...)(...)(.........)(...)(11111mnmmnnmnmxxxrxrxrxrRTnTmyyyaaa],...,[,],...,[11使用matlab命令求解：yRRRaTT*)(1或者yRa\例对下面一组数据作二次多项式拟合xi00.10.20.30.40.50.60.70.80.91yi-0.4471.9783.286.167.087.347.669.569.489.3011.2即要求出二次多项式:3221f(x)axaxa中的系数：321,,aaa,使得1012))((iiyxf达到最小解法1．用解方程组的方法1111211121xxxxR此时1）输入以下命令：x=0:0.1:1;y=[-0.4471.9783.286.167.087.347.669.569.489.3011.2];R=[(x.^2)'x'ones(11,1)]；A=R\y'计算结果：Ａ=-9.810820.1293-0.0317所求拟合多项式2()9.810820.12930.0317fxxx解法2．用多项式拟合的命令1）输入以下命令：x=0:0.1:1;y=[-0.4471.9783.286.167.087.347.669.569.489.3011.2];A=polyfit(x,y,2)z=polyval(A,x);plot(x,y,'k+',x,z,'r')%作出数据点和拟合曲线的图形2）计算结果：Ａ=-9.810820.1293-0.0317所求拟合多项式2()9.810820.12930.0317fxxx两种方法是等价的，所求拟合系数也是相同的。用MATLAB作非线性最小二乘拟合非线性最小二乘是假设)(xf是待定系数maaa,,,21的非线性函数，在最小二乘意义下求解其系数maaa,,,21。例如，前面所述常用的双曲线和指数曲线就是非线性最小二乘拟合中最常用的非线性函数，只不过在上面使用中通过变量代换将它们转变成线性最小二乘拟合方法。对于给定的实验数据，通常应根据实验数据的走向、趋势选择合适的数学模型，即拟合函数。例如，当实验数据具有单调性和凸性时，可选择下述适当的数学模型)(xfy来拟合实验数据：xbaxfaxxfaexfaexfbxbbx)(,)(,)(,)(其中ba,为待定参数。Matlab的提供了两个求非线性最小二乘拟合的函数：lsqcurvefit和lsqnonlin。两个命令都要先建立M-文件，在其中定义函数F(x)，但两者定义F(x)的方式是不同的.命令一：lsqcurvefit已知数据点：xdata=（xdata1，xdata2，…，xdatan），ydata=（ydata1，ydata2，…，ydatan）lsqcurvefit用以求含参量x（向量）的向量值函数F(x,xdata)=（F（x，xdata1），…，F（x，xdatan））T中的参变量x(向量),使得niiiydataxdataxF12)),((各种使用格式为:(1)x=lsqcurvefit(‘fun’,x0,xdata,ydata);(2)x=lsqcurvefit(‘fun’,x0,xdata,ydata,options);(3)x=lsqcurvefit(‘fun’,x0,xdata,ydata,options,’grad’);(4)[x,options]=lsqcurvefit(‘fun’,x0,xdata,ydata,…);(5)[x,options,funval]=lsqcurvefit(‘fun’,x0,xdata,ydata,…);(6)[x,options,funval,Jacob]=lsqcurvefit(‘fun’,x0,xdata,xdata,ydata,…);例：用下面一组数据拟合ktbeatc02.0)(中的参数a，b，k,数据如下:jt1002003004005006007008009001000310jc4.544.995.355.655.906.106.266.396.506.59该问题即求解最优化问题：210102.0,,][),,(minjjktkbacbeakbaFj应用lsqcurvefit求解，此时：),...,(),(10102.002.0ktktbeabeatdataxF1）编写M-文件curvefun1.mfunctionf=curvefun1(x,tdata)%x为由待拟合参数%所组成的数组，其中kxbxax)3(,)2(,)1(f=x(1)+x(2)*exp(-0.02*x(3)*tdata)%f为长度与tdata%相同的数组2)建立脚本文件qwe.mtdata=100:100:1000cdata=1e-03*[4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59];x0=[0.2,0.05,0.05];%x0可看作待拟合参数的近似值，或者算法%的初始值x=lsqcurvefit('curvefun1',x0,tdata,cdata)f=curvefun1(x,tdata)%计算在数据点上的拟合误差。命令二：lsqnonlin已知数据点：xdata=（xdata1，xdata2，…，xdatan）ydata=（ydata1，ydata2，…，ydatan）lsqnonlin用以求含参量x（向量）的向量值函数f(x)=(f1(x),f2(x),…,fn(x))T中的参量x，使得其中),,()(ydataixdataixfxfi输入格式为：1）x=lsqnonlin（‘fun’，x0）；2）x=lsqnonlin（‘fun’，x0，options）；3）x=lsqnonlin（‘fun’，x0，options，‘grad’）；4）[x，options]=lsqnonlin（‘fun’，x0，…）；5）[x，options，funval]=lsqnonlin（‘fun’，x0，…）说明:fun是一个事先建立的定义函数f(x)的M-文件，自变量为x,为由待拟和参数所构成的数组，x0为初值（或者看作待拟合参数的近似值）。例：用下面一组数据拟合ktbeatc02.0)(中的参数a，b，k,数据如下:jt1002003004005006007008009001000310jc4.544.995.355.655.906.106.266.396.506.59用命令lsqnonlin求解，此时),...,(),,(1002.0102.0101cbeacbeaydatatdataxFktkt1）编写M-文件curvefun2.mfunctionf=curvefun2(x)tdata=100:100:1000;cdata=1e-03*[4.54,4.99,5.35,5.65,5.90,…6.10,6.26,6.39,6.50,6.59];f=x(1)+x(2)*