曲线拟合的线性最小二乘法及其MATLAB程序.docx

1、曲线拟合的线性最小二乘法及其MATLAB程序1曲线拟合的线性最小二乘法及其 MATLA程序例 7.2.1给出一组数据点（人，yj列入表72中，试用线性最小二乘法求拟合曲线,并用（7.2）, （ 7.3 ）和（7.4）式估计其误差，作出拟合曲线表7 - 2例7.2.1的一组数据（xyjxi-2.5 -1.7 -1.1 -0.8 0 0.1 1.5 2.7 3.6yi-192.9 -85.50 -36.15 -26.52 -9.10 -8.43 -13.12 6.50 68.04解 (1)在MATLAB 工作窗口输入程序 x=-2.5 -1.7 -1.1 -0.8 0 0.1 1.5 2.7 3.

2、6;y=-192.9 -85.50 -36.15 -26.52 -9.10 -8.43 -13.12 6.5068.04;plot(x,y, r*),legend(实验数据(xi,yi) )xlabel( x ), ylabel( y),title(例7.2.1的数据点(xi,yi) 的散点图)运行后屏幕显示数据的散点图(略)(3)编写下列MATLAB程序计算f(x)在(xi, yi)处的函数值，即输入程序 syms a1 a2 a3 a4x=-2.5 -1.7 -1.1 -0.8 0 0.1 1.5 2.7 3.6;fi=a1.*x.A3+ a2.*x.A2+ a3.*x+ a4运行后屏幕显

3、示关于a1,a2, a3和a4的线性方程组fi = -125/8*a1+25/4*a2-5/2*a3+a4,-4913/1000*a1+289/100*a2-17/10*a3+a4,-1331/1000*a1+121/100*a2-11/10*a3+a4, -64/125*a1+16/25*a2-4/5*a3+a4, a4, 1/1000*a1+1/100*a2+1/10*a3+a4,27/8*a1+9/4*a2+3/2*a3+a4, 19683/1000*a1+729/100*a2+27/10*a3+a4, 5832/125*a1+324/25*a2+18/5*a3+a4编写构造误差平方和的

4、 MATLAB程序 y=-192.9 -85.50 -36.15 -26.52 -9.10 -8.43 -13.12 6.5068.04;fi=-125/8*a1+25/4*a2-5/2*a3+a4,-4913/1000*a1+289/100*a2-17/10*a3+a4,-1331/1000*a1+121/100*a2-11/10*a3+a4,-64/125*a1+16/25*a2-4/5*a3+a4, a4,1/1000*a1+1/100*a2+1/10*a3+a4,27/8*a1+9/4*a2+3/2*a3+a4,19683/1000*a1+729/100*a2+27/10*a3+a4,

5、5832/125*a1+324/25*a2+18/5*a3+a4;fy=fi-y; fy2=fy.A2; J=sum(fy.A2)运行后屏幕显示误差平方和如下J=(-125/8*a1+25/4*a2-5/2*a3+a4+1929/10)A2+(-4913/1000*a1+289/100*a2-17/10*a3+a4+171/2)A2+(-1331/1000*a1+121/100*a2-11/10* a3+a4+723/20)A2+(-64/125*a1+16/25*a2-4/5*a3+a4+663/25)A2+(a4+91 /10)A2+(1/1000*a1+1/100*a2+1/10*a3+

6、a4+843/100)A2+(27/8*a1+9/4*a 2+3/2*a3+a4+328/25)A2+(19683/1000*a1+729/100*a2+27/10*a3+a4-13/ 2)A2+(5832/125*a1+324/25*a2+18/5*a3+a4-1701/25)A2Q I为求a1,a2,a3,a4使J达到最小，只需利用极值的必要条件 -0 (k =1,2,3,4)，得到关于a1,a2,a3,a4的线性方程组，这可以由下面的 MATLAB程序完成，即输入程序 syms al a2 a3 a4J=(-125/8*a1+25/4*a2-5/2*a3+a4+1929/10)A2+(-

7、4913/1000*a1 +289/100*a2-17/10*a3+a4.+171/2)A2+(-1331/1000*a1+121/100*a2-11 /10*a3+a4+723/20)A2+(-64/125*a1+16/25*a2-4/5*a3+a4+663/25)A2+(a 4+91/10)A2+(1/1000*a1+1/100*a2+1/10*a3+a4+843/100)A2+(27/8*a1+9 /4*a2+3/2*a3+a4+328/25)A2+(19683/1000*a1+729/100*a2+27/10*a3+a4 -13/2)A2+(5832/125*a1+324/25*a2+

8、18/5*a3+a4-1701/25)A2;Ja1=diff(J,a1); Ja2=diff(J,a2); Ja3=diff(J,a3);Ja4=diff(J,a4);Ja11=simple(Ja1), Ja21=simple(Ja2), Ja31=simple(Ja3),Ja41=simple(Ja4),运行后屏幕显示J分别对a1, a2 ,a3 ,a4的偏导数如下Ja1仁56918107/10000*a1+32097579/25000*a2+1377283/2500*a3+ 23667/250*a4-8442429/625Ja21 =32097579/25000*a1 + 1377283/

9、2500*a2+23667/250*a3+67*a4 +767319/625Ja31 = 1377283/2500*a1+23667/250*a2+67*a3+18/5*a4-232638/125Ja41 = 23667/250*a1+67*a2+18/5*a3+18*a4+14859/25解线性方程组 Jan =0， Ja21 =0， Ja31 =0， Ja41 =0，输入下列程序A=56918107/10000, 32097579/25000, 1377283/2500,23667/250; 32097579/25000, 1377283/2500, 23667/250, 67;13772

10、83/2500, 23667/250, 67, 18/5; 23667/250, 67, 18/5, 18;B=8442429/625, -767319/625, 232638/125, -14859/25;C=B/A, f=poly2sym(C)运行后屏幕显示拟合函数 f及其系数C如下C = 5.0911 -14.1905 6.4102 -8.2574f=716503695845759/140737488355328*xA3 -7988544102557579/562949953421312*xA2 +1804307491277693/281474976710656*x -464852116

11、0813215/562949953421312故所求的拟合曲线为3 2f(x) =5.0911 x -14.1905 x 6.4102 x - 8.2574 .(4)编写下面的MATLAB程序估计其误差，并作出拟合曲线和数据的图形 .输入程序 xi=-2.5 -1.7 -1.1 -0.8 0 0.1 1.5 2.7 3.6;y=-192.9 -85.50 -36.15 -26.52 -9.10 -8.43 -13.12 6.50 68.04;n=len gth(xi);f=5.0911.*xi.A3-14.1905.*xi.A2+6.4102.*xi -8.2574;x=-2.5:0.01:

12、3.6;F=5.0911.*x.A3-14.1905.*x.A2+6.4102.*x -8.2574;fy=abs(f-y); fy2=fy.A2; Ew=max(fy),E仁 sum(fy)/n, E2=sqrt(sum(fy2)/n)plot(xi,y, r*), hold on, plot(x,F, b-), hold offlegend(数据点(xi,yi)，拟合曲线 y=f(x),xlabel( x), ylabel( y),title(例7.2.1的数据点(xi,yi)和拟合曲线y=f(x)的图形)运行后屏幕显示数据(Xj,%)与拟合函数f的最大误差 Ev,平均误差 E和均方根误差

13、 E 及其数据点(Xj,yj和拟合曲线y=f(x)的图形(略).Ew = E1 = E2 =3.105 4 0.903 4 1.240 97.3函数m(x)的选取及其MATLA程序例7.3.1 给出一组实验数据点(x )的横坐标向量为x = (-8.5,-8.7,-7.1,-6.8,-5.10,-4.5,-3.6,-3.4,-2.6,-2.5, -2.1,-1.5, -2.7,-3.6 )，纵横坐标向量为 y=( 459.26,52.81,198.27,165.60,59.17,41.66,25.92, 22.37,13.47,12.87, 11.87,6.69,14.87,24.22 )，试

14、用线性最小二乘法求拟合曲线，并用( 7.2)，( 7.3)和(7.4)式估计其误差，作出拟合曲线 .解 (1)在MATLAB工作窗口输入程序 x=-8.5,-8.7,-7.1,-6.8,-5.10,-4.5,-3.6,-3.4,-2.6,-2.5, -2.1,-1.5, -2.7,-3.6;y=459.26,52.81,198.27,165.60,59.17,41.66,25.92,22.37,13.47, 12.87, 11.87,6.69,14.87,24.22;plot(x,y, r* ),legend( 实验数据(xi,yi) )xlabel( x ), ylabel( y),titl

15、e(例7.3.1的数据点(xi,yi) 的散点图)运行后屏幕显示数据的散点图(略)(3)编写下列MATLAB程序计算f(x)在(xi, yi)处的函数值，即输入程序 syms a bx=-8.5,-8.7,-7.1,-6.8,-5.10,-4.5,-3.6,-3.4,-2.6,-2.5,-2.1,-1.5,-2.7,-3.6; fi=a.*exp(-b.*x)运行后屏幕显示关于a和b的线性方程组fi = a*exp(17/2*b), a*exp(87/10*b), a*exp(71/10*b),a*exp(34/5*b), a*exp(51/10*b), a*exp(9/2*b), a*exp

16、(18/5*b), a*exp(17/5*b), a*exp(13/5*b), a*exp(5/2*b), a*exp(21/10*b), a*exp(3/2*b), a*exp(27/10*b), a*exp(18/5*b)编写构造误差平方和的 MATLAB程序如下y=459.26,52.81,198.27,165.60,59.17,41.66,25.92,22.37,13.47,12.87, 11.87, 6.69,14.87,24.22;a*exp(71/10*b),a*exp(9/2*b), a*exp(13/5*b),a*exp(3/2*b),a,b的线性方程组,fi = a*exp

17、(17/2*b), a*exp(87/10*b),a*exp(34/5*b), a*exp(51/10*b),a*exp(18/5*b), a*exp(17/5*b),a*exp(5/2*b), a*exp(21/10*b),a*exp(27/10*b), a*exp(18/5*b);fy=fi-y;fy2=fy.A2;J=sum(fy.A2)运行后屏幕显示误差平方和如下J =(a*exp(17/2*b)-22963/50)A2+(a*exp(87/10*b)-5281/100)A2+( a*exp(71/10*b)-19827/100)A2+(a*exp(34/5*b)-828/5)A2+(

18、a*exp(51/1 0*b)-5917/100)A2+(a*exp(9/2*b)-2083/50)A2+(a*exp(18/5*b)-648/25 F2+(a*exp(17/5*b)-2237/100)A2+(a*exp(13/5*b)-1347/100)A2+(a*ex p(5/2*b)-1287/100)A2+(a*exp(21/10*b)-1187/100)A2+(a*exp(3/2*b)- 669/100F2+(a*exp(27/10*b)-1487/100)A2+(a*exp(18/5*b)-1211/50)A2为求a,b使J达到最小，只需利用极值的必要条件，得到关于 这可以由下面

19、的 MATLAB程序完成，即输入程序 syms a bJ=(a*exp(17/2*b)-22963/50)A2+(a*exp(87/10*b)-5281/100)A2 +(a*exp(71/10*b)-19827/100)A2+(a*exp(34/5*b)-828/5)A2+(a*exp(51 /10*b)-5917/100)A2+(a*exp(9/2*b)-2083/50)A2+(a*exp(18/5*b)-648/ 25)A2+(a*exp(17/5*b)-2237/100)A2+(a*exp(13/5*b)-1347/100F2+(a* exp(5/2*b)-1287/100)A2+(a

20、*exp(21/10*b)-1187/100)A2+(a*exp(3/2*b )-669/100)A2+(a*exp(27/10*b)-1487/100F2+(a*exp(18/5*b)-1211/50)A2；Ja=diff(J,a); Jb=diff(J,b);Ja仁simple(Ja), Jb仁simple(Jb),运行后屏幕显示J分别对a, b的偏导数如下Ja1 =2*a*exp(3*b)+2*a*exp(17*b)+2*a*exp(87/5*b)+2*exp(68/5*b) *a+2*exp(9*b)*a+2*a*exp(34/5*b)-669/50*exp(3/2*b)-1487/5

21、0*exp(2 7/10*b)-2507/25*exp(18/5*b)-22963/25*exp(17/2*b)-5281/50*exp(87 /10*b)-19827/50*exp(71/10*b)-2237/50*exp(17/5*b)-1656/5*exp(34/ 5*b)-1347/50*exp(13/5*b)-5917/50*exp(51/10*b)-1287/50*exp(5/2*b )-2083/25*exp(9/2*b)-1187/50*exp(21/10*b)+4*a*exp(36/5*b)+2*a*e xp(26/5*b)+2*a*exp(71/5*b)+2*a*exp(5

22、1/5*b)+2*a*exp(5*b)+2*a*exp (21/5*b)+2*a*exp(27/5*b)Jb1 =1/500*a*(2100*a*exp(21/10*b)A2+8500*a*exp(17/2*b)A2+6800 *a*exp(34/5*b)A2-10035*exp(3/2*b)-40149*exp(27/10*b)-180504*exp (18/5*b)-3903710*exp(17/2*b)-459447*exp(87/10*b)-1407717*exp(71 /10*b)-76058*exp(17/5*b)-1126080*exp(34/5*b)-35022*exp(13/

23、5*b)- 301767*exp(51/10*b)-32175*exp(5/2*b)-187470*exp(9/2*b)-24927*ex p(21/10*b)+7100*a*exp(71/10*b)A2+5100*a*exp(51/10*b)A2+4500*a*e xp(9/2*b)A2+7200*a*exp(18/5*b)A2+3400*a*exp(17/5*b)A2+2600*a*ex p(13/5*b)A2+2500*a*exp(5/2*b)A2+1500*a*exp(3/2*b)A2+2700*a*exp( 27/10*b)A2+8700*a*exp(87/10*b)A2)用解二元非

24、线性方程组的牛顿法的 MATLAB序求解线性方程组 Jai =0，Jbi =0,得a = b=2.811 0 0.581 6故所求的拟合曲线(7.13)为f(X)= 2.811 0e.5816x. (7.14)(4)根据(7.2)，( 7.3)，( 7.4)和(7.14)式编写下面的 MATLAB程序估计其误差， 并做出拟合曲线和数据的图形.输入程序 xi=-8.5 -8.7 -7.1 -6.8 -5.10 -4.5 -3.6 -3.4 -2.6 -2.5-2.1 -1.5 -2.7 -3.6;y=459.26 52.81 198.27 165.60 59.17 41.66 25.92 22.

25、3713.47 12.87 11.87 6.69 14.87 24.22;n=le ngth(xi); f=2.8110.*exp(-0.5816.*xi); x=-9:0.01: -1;F=2.8110.*exp(-0.5816.*x); fy=abs(f-y); fy2=fy.A2;Ew=max(fy),E仁 sum(fy)/n, E2=sqrt(sum(fy2)/n), plot(xi,y, r* ), hold onplot(x,F, b- ), hold off,legend(数据点(xi,yi) ，拟合曲线 y=f(x) )xlabel( x ), ylabel( y),title( 例7.3.1的数据点(xi,yi) 和拟合曲线y=f(x) 的图形)运行后屏幕显示数据 (xi , yi)与拟合函数f的最大误差 Ew = 390.141 5 ，平均误差E1=36.942 2 和均方根误差