1、 数值分析上机实验报告 选 题：曲线拟合的最小二乘法 指导老师： 专 业： 学 号： 姓 名： 昆明理工大学 数值分析上机实验报告 1 课题八 曲线拟合的最小二乘法 一、 问题提出 从随机的数据中找出其规律性，给出其近似表达式的问题，在生产实践和科学实验中大量存在，通常利用数据的最小二乘法求得拟合曲线。 在某冶炼过程中，根据统计数据的含碳量与时间关系，试求含碳量 y与时间t 的拟合曲线。 二、 要求 1、用最小二乘法进行曲线拟合； 2、近似解析表达式为 33221 tatatat ； 3、打印出拟合函数 t ，并打印出 jt 与 jty 的误差， 12,2,1 j ； 4、另外选取一个近似表达
2、式，尝试拟合效果的比较； 5、 *绘制出曲线拟合图 *。 三、 目的和意义 1、掌握曲线拟合的最小二乘法； 2、最小二乘法亦可 用于解超定线代数方程组； 3、探索拟合函数的选择与拟合精度间的关系。 四 、 计算 公式 对于给定的测量数据 (xi,fi)(i=1,2,， n)，设函数分布为 mj jj xaxy 0 )()( 特别的，取 )(xj 为多项式 jj xx )( (j=0, 1,， m) t(分 ) 0 5 10 15 20 25 30 35 40 45 50 55 y(10-4) 0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.02
3、 4.64 昆明理工大学 数值分析上机实验报告 2 则根据最小二乘法原理，可以构造泛函 ni mj ijjim xafaaaH 1 010 )(),( 令 0kaH (k=0, 1,， m) 则可以得到法方程 ),(),(),(),(),(),(),(),(),(),(),(),(1010101111000100mmmmmmmmfffaaa求该解方程组，则可以得到解 maaa , 10 ，因此可得到数据的最小二乘解 mj jj xaxf 0 )()( 曲线拟合：实际工作中，变量间未必都有线性关系，如服药后血药浓度与时间的关系；疾病疗效与疗程长短的关系；毒物剂量与致死率的关系等常呈曲线关系。曲线
4、拟合是指选择适当的曲线类型来拟合观测数据，并用拟合的曲线方程分析两变量间的关系。 五 、 结构程序设计 在程序 结构方面主要是按照顺序结构进行 设计 ， 在进行曲线的拟合时，为了进行比较，在程序设计中，直接调用了最小二乘法的拟合函数 polyfit， 并且依次调用了 plot、 figure、 hold on 函数进行图象的绘制， 最后调用了一个绝对值函数abs 用于计算拟合函数与原有数据的误差， 进行拟合效果的比较。 昆明理工大学 数值分析上机实验报告 3 5.1 用一元三次多项式 33221 tatatat 进行拟合 计算解析表达式系数： a1, a2, a3 t=0 5 10 15 20
5、 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)
6、/n) plot(xi,y,t*), hold on, plot(t,F,b-), hold off 所得函数为 4 3 3 2( t ) 0 . 3 4 3 6 4 1 0 5 . 2 1 5 6 1 0 0 . 2 6 3 4 0 . 0 1 3 8 3 9t t t 运行后屏幕显示数据 ),( ii yx 与拟合函数 f 的最大误差 Ew，平均误差 E1 和均方根误差 E2 及其数据点 ),( ii yx 和拟合曲线 y=f(x)的图形如图 5.1. Ew =0.4243 E1 =0.0911 E2 =0.1467 昆明理工大学 数值分析上机实验报告 4 图 5.1 一元三次多项式拟合曲
7、线误差图 5.2 用一元四次多项式 4433221 tatatatat 进行拟合： 计算多项式系数 ： a1, a2, a3, a4 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); 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; 昆明理工大学 数值分析上机实验报
8、告 5 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 所得函数为6 4 4 3 3 2( t ) 0 . 6 0 2 6 1 0 0 . 3 1 9 1 8 1 0 2 . 9 3 2 3 1 0 0 . 2 3 8 0 7 0 . 0 6 0 449t t t
9、t 运行后屏幕显示数据 ),( ii yx 与拟合函数 f 的最大误差 Ew，平均误差 E1 和均方根误差 E2 及其数据点 ),( ii yx 和拟合曲线 y=f(x)的图形如图 5.2。 Ew = 0.3897 E1 = 0.1034、 E2 =0.1429 图 5.2 一元四次多项式拟合曲线误差图 昆明理工大学 数值分析上机实验报告 6 5.3 用 一元二次多项式 2321 tataat 进行拟合： 计算多项式系数： a1, a2, a3 输入程序： syms a1 a2 a3 x=0 5 10 15 20 25 30 35 40 45 50 55; fi=a1.*x.2+ a2.*x+
10、 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
11、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; f
12、y=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
13、 + 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 为求 4321 , aaaa 使 J 达到最小，只需利用极值的必要条件 0kaJ)4,3,2,1( k ，得到关于 4321 , aaaa 的线性方程组，这可以由下面的 MATLAB 程序完成，即输入程序 ： 昆明理工大学 数值分析上机实验报告 7 syms a1 a2 a3 J =(100*a1 + 10*a2 + a3 - 54/25)
14、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
15、 - 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/2 Ja21 = 1089000*a1 + 25300*a2 + 660*a3 - 27131/10 Ja31 = 25300*
16、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 如下： C =-0.0024 0.2037 0.2305 F = (7338734818964133*x)/36028797018963968 - (5489104202452799*x2)/23058
17、43009213693952 + 8303449950332545/36028797018963968 故所求的拟合曲线为： 0 . 2 3 0 4 70 . 2 0 3 6 90 . 0 0 2 3 8 0 5)( 22 xxxf 编写下面的 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); 昆明理工大学 数值分析上机实验报告 8 f=-0.0023
18、805.*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)的图形 ) 运行后屏幕显示数据 ),( ii yx 与拟合
19、函数 f的最大误差 Ew，平均误差 E1 和均方根误差 E2 及其数据点 ),( ii yx 和拟合曲线 y=f(x)的图形如图 5.3 所示： Ew =0.4437 E1 = 0.1426 E2 =0.1983。 图 5.3 一元 二 次多项式拟合曲线误差图 昆明理工大学 数值分析上机实验报告 9 六 、结果讨论和分析： 由以上结果可知，拟合方程的选取至关重要，它决定了最大误差、平均误差以及均方根误差的大小，即拟合曲线的接近程度。本次实验，最初所选取的拟合解析方程 33221 tatatat 获得较好的拟合 ， 选用解析方程为 4433221 tatatatat 的曲线拟合时， 精确度有所下
20、降。由此 ，拟合函数的选择和拟合精度致密相关 ， 最小二乘法如果想将曲线拟合的比较完美，必须应用适当的模拟曲线，如果模拟曲线选择不够适当，那么用最小二乘法计算完后，会发现拟合曲线误差比较大，均方误差也比较大，而如果拟合曲线选择适当，那么效果较好 ，且根据本次结果可见 ，当采用更高次的多项式拟合数据，其结果的误差会更小 。 因此，需要对已知点根据分布规律选取多个可能的近似拟合曲线，算出后比较误差与均方误差，得到最佳拟合曲线。但是如果已知点分布非常不规律，无法观察或是无法正确观察出其近似曲线，那么根本无法使用最小二乘法进行曲线拟合，我们只能使用其它方法进行逼近。 通过 这次实验，我 学习并实践了 最小二乘法进行的曲线拟合 的知识，认识到数值分析这一分析方法在实际应用中的重要作用 。而且通过在实际 操作发现了各种问题，并寻找到解决问题的方法，使我获益良多。
