MATLAB实现非线性曲线拟合最小二乘法

1、非线性曲线拟合最小二乘法、问题提出设数据(Xj,yJ ,(i=0,1,2,3,4).由表3-1给出，表中第四行为In y二« ,可以看出数学模型为y二aebx,用最小二乘法确定a及b。i01234Xi1.001.251.501.752.00yi5.105.796.537.458.46yi1.6291.7561.8762.0082.135、理论基础根据最小二乘拟合的定义：在函数的最佳平方逼近中f(x). Ca,b，如果f(x) 只在一组离散点集Xi,i=0,1,m，上给定，这就是科学实验中经常见到的实验 数据 ( Xj, % ) , i=0,1,m的曲线拟合，这里yi二f(xj，i=0

2、,1,m，要求一个 函数y = S*(x)与所给数据 ( Xi, yi ) ,i=0,1,，m拟合，若记误差i 二 S*(xJ-% ,i=0,1，m,、=(0,1,，、m)T ,设 (x), (x)/ , :n(x)是 Ca,b上线性无关函数族，在二spar( (x), ：l(x)，,；(x)中找一函数S*(x)，使误差 平方和mmm2、2 八、i2 八S*(Xi)-y2 =min' S(Xi)-yi2,i=0i =0S(x)邯 im这里S(x)二 a。0(x) 4 !(x)an n(x)(n<m)这就是一般的最小二乘逼近，用几何语言来说，就称为曲线拟合的最小二乘 法。在建模的过

3、程中应用到了求和命令(sum)、求偏导命令(diff)、化简函数命 令(simple) >用迭代方法解二元非线性方程组的命令(fsolve)，画图命令(plot) 等。三、实验内容用最小二乘法求拟合曲线时，首先要确定S(x)的形式。这不单纯是数学问题，还与所研究问题的运动规律及所得观测数据(Xi,% )有关；通常要从问题的运动规律及给定数据描图，确定S(x)的形式，并通过实际计算选出较好的结果。S(x) 的一般表达式为线性形式，若 (x)是k次多项式，S(x)就是n次多项式，为了使问题的提法更有一般性，通常在最小二乘法中2都考虑为加权平方和m:2 八(Xi)S(Xj) - f(xj2.i

4、=0这里(X)_0是a,b上的权函数，它表示不同点(Xi，f(xj)处的数据比重不 同。用最小二乘法求拟合曲线的问题，就是求形如S(x)的一个函数y=S*(x)，2 m使I|2 =7(Xi)S(Xi) - f(Xi)2取得最小。它转化为求多元函数i z0mnl(a°,ai, ,an)八，(x)' a：j(xj - f (k)2i=0j=0的极小点(a；,a；,，a；)问题。再由求多元函数极值的必要条件，有mn=2、 ，(Xi)a(Xi) - f(xj MX) =0(k=0,1，,n)-aki =0j =0此题中假设(x)=1，由已知所给数据点(xi, yi)画出图形，根据离散

5、点的 位置观察出它们所拟合的曲线图形应类似于指数函数的曲线图形，故设拟合曲线的函数为y二aebx。本题编程过程中，令f=y,z1=a,z2=b,令拟合曲线中对应Xi的函数值与y的差的平方和为J,即J=sum(fy.A2)分别求 J关于z1,z2的偏导，简化后并令其分别为0得一关于z1,z2的二元非线性方程组， 最后利用fsolve命令求得z1,z2的值分别为z1=3.0751z2=0.5052故得到拟合曲线为c cr 厂 a 0.5052x目二 3.0751e为证明曲线拟合的正确性，我们将离散点(人，yQ与所得的拟合曲线c cr 厂 a 0.5052xy = 3.0751e画于同一图形中，图形

6、如下:例8的数据点(x(i),y(i)和拟合曲线y=f(x)的图形8.5 数值分析，李庆扬，王能超，易大义，2001，清华大学出版社(第四版)。 数值方法，关治，陆金甫，2006,清华大学出版社。 数值分析与实验学习指导，蔡大用，2001，清华大学出版社。 数值分析与实验，薛毅，2005,北京工业大学出版社.r1 r1r 11：-+ 数据点(x(i),y(i) 拟合曲线y=f(x)8- -7.5 -7-6.5 -6-5.5 -c |L|5数据点(x(i),y(i) 拟合曲线y=f(x)11.11.21.31.41.51.61.71.81.92x轴四、结果分析根据实验内容求得拟合曲线y = ae

7、bx中未知数a,b分别为a=3.0751b=0.5052即拟合曲线为y = 3.0751e0.5052x。由图形知拟合成功！参考文献附录程序 1: syms z1 z2x=1.00:0.25:2.00;y=5.10,5.79,6.53,7.45,8.46;f=z1*exp(z2.*x)fy=f-y;J=sum(fy.A2);Ja=diff(J,z1);Jb=diff(J,z2);Ja1=simple(Ja),Jb1=simple(Jb),程序 2:function y1,y2=fun(z)y1=2*z(1)*exp(2*z(2)-51/5*exp(z(2)+2*z(1)*exp(5/2*z(2

8、)-579/50*exp( 5/4*z(2)+2*z(1)*exp(3*z(2)-653/50*exp(3/2*z(2)+2*z(1)*exp(7/2*z(2) -149/10*exp(7/4*z (2) )+2*exp(4*z(2)*z(1)-423/25*exp(2*z (2);y2=-1/200*z(1)*(-400*z(1)*exp(2*z (2) )+2040*exp(z(2)-500*z(1)*exp(5/2 *z(2)+2895*exp(5/4*z(2)-600*z(1)*exp(3*z (2) )+3918*exp(3/2*z(2)-700 *z(1)*exp(7/2*z (2) )+5215*exp(7/4*z(2)-800*z(1)*exp(4*z (2) )+6768*exp (2*z (2);作图程序：x=1.00:0.25:2.00;y=5.10,5.79,6.53,7.45,8.46; f=3.0751*exp(0.5052*x); plot(x,y,'r*',x,f,'b-');xlabel(