数值分析上机作业(8页).doc
-数值分析上机作业-第 6 页数值分析上机实验报告选 题：曲线拟合的最小二乘法指导老师： 专 业： 学 号： 姓 名： 课题八 曲线拟合的最小二乘法一、问题提出从随机的数据中找出其规律性，给出其近似表达式的问题，在生产实践和科学实验中大量存在，通常利用数据的最小二乘法求得拟合曲线。在某冶炼过程中，根据统计数据的含碳量与时间关系，试求含碳量y与时间t的拟合曲线。t(分)0 5 10 15 20 25 30 35 40 45 50 55y(×10-4)二、要求1、用最小二乘法进行曲线拟合；2、近似解析表达式为；3、打印出拟合函数，并打印出与的误差，；4、另外选取一个近似表达式，尝试拟合效果的比较；5、*绘制出曲线拟合图*。三、目的和意义1、掌握曲线拟合的最小二乘法；2、最小二乘法亦可用于解超定线代数方程组；3、探索拟合函数的选择与拟合精度间的关系。四、计算公式对于给定的测量数据(xi,fi)(i=1,2,，n)，设函数分布为特别的，取为多项式 (j=0, 1,，m)则根据最小二乘法原理，可以构造泛函令 (k=0, 1,，m)则可以得到法方程求该解方程组，则可以得到解，因此可得到数据的最小二乘解曲线拟合：实际工作中，变量间未必都有线性关系，如服药后血药浓度与时间的关系；疾病疗效与疗程长短的关系；毒物剂量与致死率的关系等常呈曲线关系。曲线拟合是指选择适当的曲线类型来拟合观测数据，并用拟合的曲线方程分析两变量间的关系。五、结构程序设计在程序结构方面主要是按照顺序结构进行设计，在进行曲线的拟合时，为了进行比较，在程序设计中，直接调用了最小二乘法的拟合函数polyfit，并且依次调用了plot、figure、hold on函数进行图象的绘制，最后调用了一个绝对值函数abs用于计算拟合函数与原有数据的误差，进行拟合效果的比较。进行拟合计算解析表达式系数： a1, a2, a3t=0 5 10 15 20 25 30 35 40 45 50 55;y=0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.02 4.64;>> n=length(xi);f=0.34364.*10.(-4)*x.3-5.2156.*10.(-3)*x.2+0.26340.*x+0.017839; x=0:0.01:55; F=0.34364.*10.(-4)*x.3-5.2156.*10.(-3)*x.2+0.26340.*x+0.017839;fy=abs(f-y);fy2=fy.2;Ew=max(fy),E1=sum(fy)/n,E2=sqrt(sum(fy2)/n)plot(xi,y,'t*'), hold on, plot(t,F,'b-'), hold off 所得函数为运行后屏幕显示数据与拟合函数f的最大误差Ew，平均误差E1和均方根误差E2及其数据点和拟合曲线y=f(x)的图形如图.Ew图一元三次多项式拟合曲线误差图进行拟合：计算多项式系数：a1, a2, a3, a4xi=0 5 10 15 20 25 30 35 40 45 50 55;y=0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.02 4.64;n=length(xi);x=0:0.01:55;f=0.6026.*10.(-6)*x.4-0.31918.*10.(-4)*x.3-0.0029323.*x.2+0.23807.*x+0.060449; x=0:0.01:55;F=0.6026.*10.(-6)*x.4-0.31918.*10.(-4)*x.3-0.0029323.*x.2+0.23807.*x+0.060449;fy=abs(f-y);fy2=fy.2;Ew=max(fy),E1=sum(fy)/n,E2=sqrt(sum(fy2)/n)plot(xi,y,'r*'), hold on, plot(x,F,'b-'), hold off 所得函数为运行后屏幕显示数据与拟合函数f的最大误差Ew，平均误差E1和均方根误差E2及其数据点和拟合曲线y=f(x)的图形如图5.2。Ew图一元四次多项式拟合曲线误差图用一元二次多项式进行拟合：计算多项式系数：a1, a2, a3输入程序：>> syms a1 a2 a3x=0 5 10 15 20 25 30 35 40 45 50 55;fi=a1.*x.2+ a2.*x+ a3运行后屏幕显示关于a1,a2和a3的线性方程组：fi= a3, 25*a1 + 5*a2 + a3, 100*a1 + 10*a2 + a3, 225*a1 + 15*a2 + a3, 400*a1 + 20*a2 + a3, 625*a1 + 25*a2 + a3, 900*a1 + 30*a2 + a3, 1225*a1 + 35*a2 + a3, 1600*a1 + 40*a2 + a3, 2025*a1 + 45*a2 + a3, 2500*a1 + 50*a2 + a3, 3025*a1 + 55*a2 + a3编写构造误差平方和的MATLAB程序：y=0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.02 4.64;fi = a3, 25*a1 + 5*a2 + a3, 100*a1 + 10*a2 + a3, 225*a1 + 15*a2 + a3, 400*a1 + 20*a2 + a3, 625*a1 + 25*a2 + a3, 900*a1 + 30*a2 + a3, 1225*a1 + 35*a2 + a3, 1600*a1 + 40*a2 + a3, 2025*a1 + 45*a2 + a3, 2500*a1 + 50*a2 + a3, 3025*a1 + 55*a2 + a3;fy=fi-y;fy2=fy.2;J=sum(fy.2) 运行后屏幕显示误差平方和如下：J =(100*a1 + 10*a2 + a3 - 54/25)2 + (25*a1 + 5*a2 + a3 - 127/100)2 + (225*a1 + 15*a2 + a3 - 143/50)2 + (400*a1 + 20*a2 + a3 - 86/25)2 + (900*a1 + 30*a2 + a3 - 83/20)2 + (625*a1 + 25*a2 + a3 - 387/100)2 + (1225*a1 + 35*a2 + a3 - 437/100)2 + (1600*a1 + 40*a2 + a3 - 451/100)2 + (2025*a1 + 45*a2 + a3 - 229/50)2 + (2500*a1 + 50*a2 + a3 - 201/50)2 + (3025*a1 + 55*a2 + a3 - 116/25)2 + a32为求使达到最小，只需利用极值的必要条件 ，得到关于的线性方程组，这可以由下面的MATLAB程序完成，即输入程序 ：>> syms a1 a2 a3J =(100*a1 + 10*a2 + a3 - 54/25)2 + (25*a1 + 5*a2 + a3 - 127/100)2 + (225*a1 + 15*a2 + a3 - 143/50)2 + (400*a1 + 20*a2 + a3 - 86/25)2 + (900*a1 + 30*a2 + a3 - 83/20)2 + (625*a1 + 25*a2 + a3 - 387/100)2 + (1225*a1 + 35*a2 + a3 - 437/100)2 + (1600*a1 + 40*a2 + a3 - 451/100)2 + (2025*a1 + 45*a2 + a3 - 229/50)2 + (2500*a1 + 50*a2 + a3 - 201/50)2 + (3025*a1 + 55*a2 + a3 - 116/25)2 + a32;Ja1=diff(J,a1);Ja2=diff(J,a2);Ja3=diff(J,a3);Ja11=simple(Ja1),Ja21=simple(Ja2),Ja31=simple(Ja3),运行后屏幕显示J分别对a1, a2 ,a3的偏导数如下：Ja11 =49967500*a1 + 1089000*a2 + 25300*a3 - 217403/2Ja21 = 1089000*a1 + 25300*a2 + 660*a3 - 27131/10Ja31 = 25300*a1 + 660*a2 + 24*a3 - 3987/50解线性方程组Ja11 =0，Ja21 =0，Ja31 =0输入下列程序：>>A=49967500,1089000,25300;1089000, 25300,660;25300,660,24;B=217403/2,27131/10,3987/50;C=B/A, F=poly2sym(C)运行后屏幕显示拟合函数f及其系数C如下：故所求的拟合曲线为：编写下面的MATLAB程序估计其误差，并作出拟合曲线和数据的图形。输入程序：>> xi=0 5 10 15 20 25 30 35 40 45 50 55;y=0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.02 4.64;n=length(xi);f=-0.0023805.*x.2+0.20369.*x+0.23047;x=0:0.01:55;F=-0.0023805.*x.2+0.20369.*x+0.23047;fy=abs(f-y);fy2=fy.2;Ew=max(fy),E1=sum(fy)/n,E2=sqrt(sum(fy2)/n)plot(xi,y,'r*'), hold on, plot(x,F,'b-'), hold off legend('数据点(xi,yi)','拟合曲线y=f(x)'),xlabel('x'),ylabel('y'),title('数据点(xi,yi)和拟合曲线y=f(x)的图形')运行后屏幕显示数据与拟合函数f的最大误差Ew，平均误差E1和均方根误差E2及其数据点和拟合曲线y=f(x)的图形如图所示：Ew图5.3一元二次多项式拟合曲线误差图六、结果讨论和分析： 由以上结果可知，拟合方程的选取至关重要，它决定了最大误差、平均误差以及均方根误差的大小，即拟合曲线的接近程度。本次实验，最初所选取的拟合解析方程获得较好的拟合，选用解析方程为的曲线拟合时，精确度有所下降。由此，拟合函数的选择和拟合精度致密相关，最小二乘法如果想将曲线拟合的比较完美，必须应用适当的模拟曲线，如果模拟曲线选择不够适当，那么用最小二乘法计算完后，会发现拟合曲线误差比较大，均方误差也比较大，而如果拟合曲线选择适当，那么效果较好，且根据本次结果可见，当采用更高次的多项式拟合数据，其结果的误差会更小。因此，需要对已知点根据分布规律选取多个可能的近似拟合曲线，算出后比较误差与均方误差，得到最佳拟合曲线。但是如果已知点分布非常不规律，无法观察或是无法正确观察出其近似曲线，那么根本无法使用最小二乘法进行曲线拟合，我们只能使用其它方法进行逼近。通过这次实验，我学习并实践了最小二乘法进行的曲线拟合的知识，认识到数值分析这一分析方法在实际应用中的重要作用。而且通过在实际操作发现了各种问题，并寻找到解决问题的方法，使我获益良多。
