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

1、曲线拟合的线性最小二乘法及其MATLAB程序1 曲线拟合的线性最小二乘法及其MATLAB程序例7.2.1 给出一组数据点列入表72中，试用线性最小二乘法求拟合曲线，并用（7.2），（7.3）和（7.4）式估计其误差，作出拟合曲线.表72 例7.2.1的一组数据xi-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.6; y=-192.9 -

2、85.50 -36.15 -26.52 -9.10 -8.43 -13.12 6.50 68.04;plot(x,y,r*),legend(实验数据(xi,yi)xlabel(x), ylabel(y),title(例7.2.1的数据点(xi,yi)的散点图)运行后屏幕显示数据的散点图（略）. （3）编写下列MATLAB程序计算在处的函数值，即输入程序 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.3+ a2.*x.2+ a3.*x+ a4运行后屏幕显示关于a1,a2, a3和a4的线性方程组fi = -12

3、5/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编写构造误差平方和的MATLAB程序 y=-192.9 -85.50 -

4、36.15 -26.52 -9.10 -8.43 -13.12 6.50 68.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, 5832/125*a1+324/25*

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

6、2+3/2*a3+a4+328/25)2+(19683/1000*a1+729/100*a2+27/10*a3+a4-13/2)2+(5832/125*a1+324/25*a2+18/5*a3+a4-1701/25)2为求使达到最小，只需利用极值的必要条件 ，得到关于的线性方程组，这可以由下面的MATLAB程序完成，即输入程序 syms a1 a2 a3 a4J=(-125/8*a1+25/4*a2-5/2*a3+a4+1929/10)2+(-4913/1000*a1+289/100*a2-17/10*a3+a4.+171/2)2+(-1331/1000*a1+121/100*a2-11/10

7、*a3+a4+723/20)2+(-64/125*a1+16/25*a2-4/5*a3+a4+663/25)2+(a4+91/10)2+(1/1000*a1+1/100*a2+1/10*a3+a4+843/100)2+(27/8*a1+9/4*a2+3/2*a3+a4+328/25)2+(19683/1000*a1+729/100*a2+27/10*a3+a4-13/2)2+(5832/125*a1+324/25*a2+18/5*a3+a4-1701/25)2; Ja1=diff(J,a1); Ja2=diff(J,a2); Ja3=diff(J,a3); Ja4=diff(J,a4);Ja1

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

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

10、5, -14859/25; C=B/A, f=poly2sym(C)运行后屏幕显示拟合函数f及其系数C如下C = 5.0911 -14.1905 6.4102 -8.2574f=7165/1428*x3-7988544102557579/562949953421312*x2+77693/2856*x-4648521160813215/562949953421312 故所求的拟合曲线为.（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

11、-26.52 -9.10 -8.43 -13.12 6.50 68.04;n=length(xi); f=5.0911.*xi.3-14.1905.*xi.2+6.4102.*xi -8.2574;x=-2.5:0.01: 3.6; F=5.0911.*x.3-14.1905.*x.2+6.4102.*x -8.2574;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 offlegend(数据点(xi,yi),拟合曲线y=

12、f(x), xlabel(x), ylabel(y),title(例7.2.1的数据点(xi,yi)和拟合曲线y=f(x)的图形)运行后屏幕显示数据与拟合函数f的最大误差Ew，平均误差E1和均方根误差E2及其数据点和拟合曲线y=f(x)的图形（略）.Ew = E1 = E2 =3.105 4 0.903 4 1.240 97.3 函数的选取及其MATLAB程序例7.3.1 给出一组实验数据点的横坐标向量为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

13、.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)，试用线性最小二乘法求拟合曲线，并用（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.

14、69,14.87,24.22;plot(x,y,r*),legend(实验数据(xi,yi)xlabel(x), ylabel(y),title(例7.3.1的数据点(xi,yi)的散点图)运行后屏幕显示数据的散点图（略）.（3）编写下列MATLAB程序计算在处的函数值，即输入程序 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(

15、71/10*b), a*exp(34/5*b), a*exp(51/10*b), a*exp(9/2*b), a*exp(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;fi = a*exp(17/2*b

16、), 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(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);fy=fi-y;fy2=fy.2;J=sum(fy.2)运行后屏幕显示误差平方和如下J = (a*exp(17/2*b)-22963/50)2+(a*exp(87/10*b)-5281/100)2+(a

17、*exp(71/10*b)-19827/100)2+(a*exp(34/5*b)-828/5)2+(a*exp(51/10*b)-5917/100)2+(a*exp(9/2*b)-2083/50)2+(a*exp(18/5*b)-648/25)2+(a*exp(17/5*b)-2237/100)2+(a*exp(13/5*b)-1347/100)2+(a*exp(5/2*b)-1287/100)2+(a*exp(21/10*b)-1187/100)2+(a*exp(3/2*b)-669/100)2+(a*exp(27/10*b)-1487/100)2+(a*exp(18/5*b)-1211/5

18、0)2为求使达到最小，只需利用极值的必要条件，得到关于的线性方程组，这可以由下面的MATLAB程序完成，即输入程序 syms a bJ=(a*exp(17/2*b)-22963/50)2+(a*exp(87/10*b)-5281/100)2+(a*exp(71/10*b)-19827/100)2+(a*exp(34/5*b)-828/5)2+(a*exp(51/10*b)-5917/100)2+(a*exp(9/2*b)-2083/50)2+(a*exp(18/5*b)-648/25)2+(a*exp(17/5*b)-2237/100)2+(a*exp(13/5*b)-1347/100)2+(

19、a*exp(5/2*b)-1287/100)2+(a*exp(21/10*b)-1187/100)2+(a*exp(3/2*b)-669/100)2+(a*exp(27/10*b)-1487/100)2+(a*exp(18/5*b)-1211/50)2;Ja=diff(J,a); Jb=diff(J,b); Ja1=simple(Ja), Jb1=simple(Jb),运行后屏幕显示J分别对的偏导数如下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/

20、50*exp(3/2*b)-1487/50*exp(27/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*exp(26/5*b)+2*a*exp(71

21、/5*b)+2*a*exp(51/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)2+8500*a*exp(17/2*b)2+6800*a*exp(34/5*b)2-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

22、*exp(13/5*b)-301767*exp(51/10*b)-32175*exp(5/2*b)-187470*exp(9/2*b)-24927*exp(21/10*b)+7100*a*exp(71/10*b)2+5100*a*exp(51/10*b)2+4500*a*exp(9/2*b)2+7200*a*exp(18/5*b)2+3400*a*exp(17/5*b)2+2600*a*exp(13/5*b)2+2500*a*exp(5/2*b)2+1500*a*exp(3/2*b)2+2700*a*exp(27/10*b)2+8700*a*exp(87/10*b)2)用解二元非线性方程组的牛

23、顿法的MATLAB程序求解线性方程组Ja1 =0，Jb1 =0，得a = b=2.811 0 0.581 6故所求的拟合曲线（7.13）为e. （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.37 13.47 12.87 11.87 6.69 14.8

24、7 24.22;n=length(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.2; Ew=max(fy),E1=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)的图形)运行后屏幕显示

25、数据与拟合函数f的最大误差Ew = 390.141 5，平均误差E1=36.942 2和均方根误差E2=106.031 7及其数据点和拟合曲线y=f(x)的图形（略）.7.4 多项式拟合及其MATLAB程序例7.4.1 给出一组数据点列入表73中，试用线性最小二乘法求拟合曲线，并用（7.2），（7.3）和（7.4）式估计其误差，作出拟合曲线.表73 例7.4.1的一组数据xi-2.9 -1.9 -1.1 -0.8 0 0.1 1.5 2.7 3.6yi53.94 33.68 20.88 16.92 8.79 8.98 4.17 9.12 19.88解 （1）首先根据表73给出的数据点，用下列M

26、ATLAB程序画出散点图.在MATLAB工作窗口输入程序 x=-2.9 -1.9 -1.1 -0.8 0 0.1 1.5 2.7 3.6; y=53.94 33.68 20.88 16.92 8.79 8.98 4.17 9.12 19.88;plot(x,y,r*), legend(数据点(xi,yi)xlabel(x), ylabel(y), title(例7.4.1的数据点(xi,yi)的散点图)运行后屏幕显示数据的散点图（略）.（3）用作线性最小二乘拟合的多项式拟合的MATLAB程序求待定系数 .输入程序 a=polyfit(x,y,2)运行后输出（7.16）式的系数a =2.8302

27、 -7.3721 9.1382故拟合多项式为 .（4）编写下面的MATLAB程序估计其误差，并做出拟合曲线和数据的图形.输入程序 xi=-2.9 -1.9 -1.1 -0.8 0 0.1 1.5 2.7 3.6; y=53.94 33.68 20.88 16.92 8.79 8.98 4.17 9.12 19.88;n=length(xi); f=2.8302.*xi.2-7.3721.*xi+9.1382x=-2.9:0.001:3.6;F=2.8302.*x.2-7.3721.*x+8.79;fy=abs(f-y); fy2=fy.2; Ew=max(fy), E1=sum(fy)/n,E

28、2=sqrt(sum(fy2)/n), plot(xi,y,r*, x,F,b-),legend(数据点(xi,yi),拟合曲线y=f(x)xlabel(x), ylabel(y), title(例7.4.1 的数据点(xi,yi)和拟合曲线y=f(x)的图形)运行后屏幕显示数据与拟合函数f的最大误差Ew，平均误差E1和均方根误差E2及其数据点(xi,yi)和拟合曲线y=f(x)的图形（略）.Ew = E1 = E2 =0.745 7, 0.389 2, 0.436 37.5 拟合曲线的线性变换及其MATLAB程序 例7.5.1 给出一组实验数据点的横坐标向量为x=（7.5 6.8 5.10

29、4.5 3.6 3.4 2.6 2.5 2.1 1.5 2.7 3.6），纵横坐标向量为y=（359.26 165.60 59.17 41.66 25.92 22.37 13.47 12.87 11.87 6.69 14.87 24.22），试用线性变换和线性最小二乘法求拟合曲线，并用（7.2），（7.3）和（7.4）式估计其误差，作出拟合曲线.解 （1）首先根据给出的数据点，用下列MATLAB程序画出散点图.在MATLAB工作窗口输入程序 x=7.5 6.8 5.10 4.5 3.6 3.4 2.6 2.5 2.1 1.5 2.7 3.6;y=359.26 165.60 59.17 41.6

30、6 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),title(例7.5.1的数据点(xi,yi)的散点图)运行后屏幕显示数据的散点图（略）.（2）根据数据散点图，取拟合曲线为 e , （7.19） 其中是待定系数.令，则（7.19）化为.在MATLAB工作窗口输入程序 x=7.5 6.8 5.10 4.5 3.6 3.4 2.6 2.5 2.1 1.5 2.7 3.6;y=359.26 165.60 59.17 41.66 25.92 22.37 13.47 12.87 11.87 6.69 14.87 24.22; Y=log(y); a=polyfit(x,Y,1); B=a(1);A=a(2); b=B,a=exp(A)n=length(x); X=8:-0.01:1; 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(例7.5.1 的数据点(xi,yi)和拟合曲线y=f(x)的图形)fy=abs(f-y)