基于matlab的最小二乘法应用

基于matlab非线性曲线最小二乘拟摘要：在工程计算与科学研究中，常常需要从一组测量数据出发，寻找变量x与y的函数关系式，有时很难找出他们之间精确地函数表达式)(xfy，这时就要观察所给数据值，利用最小二乘曲线拟合来构造一个近似的解析式)()(xxfy。利用这种方法拟合出的曲线虽然不能保证通过所有的数据点，但是很好的逼近它们，从而充分反映已知数据间内在的数量关系。因此这种方法在科学实验和生产实践中具有广泛的应用前景。一般构造)(xy的方法有很多，本文先介绍了最小二乘法的原理，并通过实例用matlab实现了曲线的拟合以得到函数关系的方法和步骤。通过不同的经验公式得到不同的的拟合结果，并分析其结果。关键字：最小二乘法matlab曲线拟合1前言在现代科学研究中,物理量之间的相互关系通常是用函数来描述的。有些函数关系是由经典理论分析推导得出的,这些函数关系不仅为我们进一步的分析研究工作提供了物理的理论基础,也使我们可以十分方便的运用丰富的数学知识来解决物理问题。在现实的物理研究过程中,有一些问题很难由经典物理理论推导出物理量的函数表达式,或者推导出的表达式十分复杂,不利于进一步的分析。但由于研究需要,又很希望能得到这些量之间的函数关系,这时就可以利用曲线拟合的方法,用实验数据结合数学方法得到物理量之间的近似函数表达式。2最小二乘法原理在函数的最佳平方逼近中),()(baxf，如果)(xf只在一组离散数据点集{ix,mi,.....,1,0}上给出，这就是科学实验中经常见到的实验数据{（iiyx,），mi,.....,1,0}的曲线拟合，这里)(iixfy（mi,.....,1,0），要求一个函数xys与所给数据{（iiyx,），mi,.....,1,0}拟合，若误差iiyxs（mi,.....,1,0），Tm,.....,,10，设baCxxxn,,......,10是上线性无关函数族，在xxxxspansn中找一函数,.....,,10，使得误差平方和,min20200222miiixSmiiimiiyxSyxs（2-1）这里xaxaxaxSnn...1100nm(2-2)这就是一般的最小二乘逼近，用几何语言说，就称为曲线拟合的最小二乘法[1]。用最小二乘法求拟合曲线时，首先要确定S(x)的形式。这不单纯是数学问题，还与所研究问题的运动规律及所得观测数据（iiyx,）有关；通常要从问题的运动规律及给定数据描图，确定S(x)的形式，并通过实际计算选出较好的结果。S(x)的一般表达式为线性形式，若)(xk是k次多项式，S(x)就是n次多项式，为了使问题的提法更有一般性，通常在最小二乘法中22都考虑为加权平方和miiiixfxSx0222)]()()[(（2-3）这里0)(x是[a,b]上的权函数，它表示不同点))(,(iixfx处的数据比重不同一般取值为1.用最小二乘法求拟合曲线的问题，就是求形如S(x)的一个函数)(*xSy，使miiiixfxSx0222)]()()[(取得最小。它转化为求多元函数20010)]()([)(),,,(injijjmiinxfxaxaaaI（2-4）的极小点),,,(**1*0naaa问题。再由求多元函数极值的必要条件，有miiknjiijjikxxfxaxaI000)(])()()[(2(k=0,1,…,n)若记)()()(),(0ikijmiikjxxx（2-5）kikimiikdxxfxf)()()(),(0上式可改写成),...,1,0(,),(nkdakjnojjk（2-6）这个方程成为法方程，可写成距阵形式dGa（2-7）其中,),...,,(,),...,,(1010TnTnddddaaaa),(),(),()(),(),(),(),(),(101110101000nnnnnnG（2-8）它的均方误差为：.)]()([)(||||2022iimiixfxSx(2-9)最大偏差为：yfEw（2-10）3问题描述在化学反应中，由实验测得分解物浓度与时间的关系如下表2所示表2浓度（y）与时间(x)的关系实验数据表x0510152025y01.272.162.863.443.87x303540455055y4.154.374.514.584.624.64（1）试画出以上数据的散点图；（2）根据数据的变化趋势，使用合适的经验公式)(tf拟合以上数据，求出偏差和均方误差，并说明优劣。4问题求解及分析过程由所给的实验数据点通过orgin软件做出如图1所示的散点图，从图中可以01020304050012345xy图1.分解物浓度与时间的散点图及问题的物理背景可以看出，拟合曲线xy应具有如下特点：①曲线随时间的增加而上升，但上升速度由快到慢。②当x=0时反应尚未开始，即y=0当x时，y趋于某一常数。故曲线应通过原点，且有一水平渐近线。具有这两种特点的曲线有很多，悬着不同的数学模型可以后的不同的拟合曲线和经验公式本文提供如下两种方案。方案一：设想xy具有指数形式)0,0)(/exp(baxbay（4-1）在求解参数a和b时，为了避免求解一个非线性方程组，对上式两边取对数方程变为：xbaylnln（4-2）引入新的变量xXyY/1,ln，并记bBaA,ln上式变为BXAY,此时的问题就转化为求形如BXAY的最小二乘解。运用matlab语言编写计算和画图程序，程序一见附录部分[2]。运算的结果：a=5.2151,b=-7.4962。最大偏差3044.0wE，均方误差5811.02E。故拟合的曲线为：xy/4962.7exp2151.5（4-3）拟合的曲线图形如图2所示：010203040506000.511.522.533.544.55xy数据点(xi,yi)和拟合曲线y=f(x)的图形数据点(xi,yi)拟合曲线y=f(x)图2.方案二：设想xy是双曲线型的，具有下面的形式：baxxy（4-4）此时如直接按最小二乘原则来确定参数a和b的值，问题则归结为求二元函数的极小点，这将导致要求解非线性方程组，从而给计算带来麻烦，我们可以通过变量替换将其转化为关于求解待定参数的线性函数。为此将（4-3）式转化为xbay1（4-5）引入新变量xXyY1,1，方程变为：bXaY(4-6)问题同样转变为求解形如bxay的最小二乘解。同样运用matlab语言来编写其计算和作图的程序，其程序见附录第二部分。计算结果1074.4,0110.3,2063.3,1413.02EEbaw。此种方案的拟合曲线方程为2063.31413.0xxy（4-7）拟合曲线图形见图3所示：01020304050600123456xy数据点(xi,yi)和拟合曲线y=f(x)的图形数据点(xi,yi)拟合曲线y=f(x)图3.把两种方案的所得到的经验方程（4-3）和（4-7）进行比较见表4-1，从最大偏差和均方误差这两个角度来看前者优于后者。表4-1经验公式最大偏差均方误差式（4-3）0.30440.5811式（4-7）3.01104.1074因此在解决实际问题时，常常要经过反复的分析，多次选择、计算与比较，最终才能得到较好的数学模型。5参考文献[1]李庆扬，王能超，易大义.数值分析[M].北京：清华大学出版社，2008.[2]郑阿奇，曹戈。Matlab实用教程[M].北京：电子工业出版社，2007.6附录程序1：x=[510152025303540455055];y=[1.272.162.863.443.874.154.374.514.584.624.64];Y=log(y);X=1./x;a=polyfit(X,Y,1);B=a(1);A=a(2);b=B,a=exp(A),X=0:1:55;Y=a*exp(b./X);f=a*exp(b./x);plot(x,y,'r*',X,Y,'b-'),xlabel('x'),ylabel('y')legend('数据点(xi,yi)','拟合曲线y=f(x)')title('数据点(xi,yi)和拟合曲线y=f(x)的图形')fy=abs(f-y);fy2=fy.^2;Ew=max(fy),E1=sum(fy),E2=sqrt(sum(fy2))程序2：x=[510152025303540455055];y=[1.272.162.863.443.874.154.374.514.584.624.64];Y=1./y;X=1./x;a=polyfit(X,Y,1);B=a(1);A=a(2);b=B,a=AX=0:1:55;Y=X./(a*X+b);f=x/(a*x+b);plot(x,y,'r*',X,Y,'b-'),xlabel('x'),ylabel('y')legend('数据点(xi,yi)','拟合曲线y=f(x)')title('数据点(xi,yi)和拟合曲线y=f(x)的图形')fy=abs(f-y);fy2=fy.^2;Ew=max(fy),E1=sum(fy),E2=sqrt(sum(fy2))