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

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

《曲线拟合的线性最小二乘法及其MATLAB程序.docx》由会员分享，可在线阅读，更多相关《曲线拟合的线性最小二乘法及其MATLAB程序.docx（19页珍藏版）》请在冰豆网上搜索。

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

1

曲线拟合的线性最小二乘法及其MATLA程序

例7.2.1

给出一组数据点（人，yj列入表7^2中，试用线性最小二乘法求拟合曲线,

并用（7.2）,（7.3）和（7.4）式估计其误差，作出拟合曲线

表7-2例7.2.1的一组数据（x^yj

xi

-2.5-1.7-1.1-0.800.11.52.73.6

yi

-192.9-85.50-36.15-‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎26.52-9.10-8.43-13.126.5068.04

解

（1）在MATLAB工作窗口输入程序

>>x=[-2.5-1.7-1.1-0.800.11.52.73.6];

y=[-192.9-85.50-36.15-26.52-9.10-8.43-13.126.50

68.04];

plot（x,y,'r*'）,

legend（'实验数据（x‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎i,yi）'）

xlabel（'x'）,ylabel（'y'）,

title（'例7.2.1的数据点（xi,yi）的散点图’）

运行后屏幕显示数据的散点图（略）

（3）编写下列MATLAB程序计算f（x）在（xi,yi）处的函数值，即输入程序

>>symsa1a2a3a4

x=[-2.5-1.7-1.1-0.800.11.52.73.6];

fi=a1.*x‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎.A3+a2.*x.A2+a3.*x+a4

运行后屏幕显示关于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]

编写构造误差平方和的MATLAB程序

>>y=[-192.9-85.50-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*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+2

89/1‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎00*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+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+18/5*a3+a4-1701/25）A2

QI

为求a1,a2,a3,a4使J达到最小，只需利用极值的必要条件——-0（k=1,2,3,4），

得到关于a1,a2,a3,a4的线性方程组，这可以由下面的MATLAB程序完成，即‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎输入程序

>>symsala2a3a4

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+6‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎63/25）A2+（a4+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+18/5*a3+a4-1701/25）A2;

Ja1=d‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎iff（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/625

Ja21=

32097579/25000*a1+1377283/2500*a2+23667/250*a3+67*a4+767319/625

Ja31=1377283/2500*a1+23667/250*a2+67*a3+18/5*a4-232638/125

Ja41=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;

1377283/2500,2‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎3667/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.19056.4102-8.2574

f=716503695845759/140737‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎488355328*xA3-7988544102557579/562949953421312*xA2+1804307491277693/281474976710656*x-4648521160813215/562949953421312

故所求的拟合曲线为

32

f（x）=5.0911x-14.1905x6.4102x-8.2574.

（4）编写下面的MATLAB程序估计‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎其误差，并作出拟合曲线和数据的图形.输入程

序

>>xi=[-2.5-1.7-1.1-0.800.11.52.73.6];

y=[-192.9-85.50-36.15-26.52-9.10-8.43-13.126.5068.04];

n=length（xi）;

f=5.0911.*xi.A3-14.1905.*xi.A2+6.4102.*xi-8.2574;‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎

x=-2.5:

0.01:

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*'）,holdon,plot（x,F,'b-'）,holdoff

legend（'‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎数据点（xi,yi）'，拟合曲线y=f（x）'）,

xlabel（'x'）,ylabel（'y'）,

title（'例7.2.1的数据点（xi,yi）和拟合曲线y=f（x）的图形'）

运行后屏幕显示数据（Xj,%）与拟合函数f的最大误差Ev,平均误差E和均方根误差E及其数据点（Xj,yj和拟合曲线y=f（x）的图形（略）.

Ew=E1=E2=

3.10540.903‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎41.2409

7.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,2‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎2.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.69,14.87,24.22];

plot（x,y,'r*'）,legend（'实验数据（xi,yi）'）

xlabel（'x'）,ylabel（'y'）,

title（'例7.3.1的数据点（xi,yi）的散点图’）

运行后屏幕‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎显示数据的散点图（略）

（3）编写下列MATLAB程序计算f（x）在（xi,yi）处的函数值，即输入程序

>>symsab

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];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（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/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）-1‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎9827/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/25F2+（a*exp（17/5*b）-2237/100）A2+（a*exp（13/5*b）-1347/100）A2+（a*exp（5/2*b）-1287/100）A2+（a*e‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎xp（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）A

2

为求a,b使J达到最小，只需利用极值的必要条件，得到关于这可以由下面的MATLAB程序完成，即输入程序

>>symsab

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+（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*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仁simp‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎le（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/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（3‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎6/5*b）+2*a*exp（26/5*b）+2*a*exp（71/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）A2+8500*a*exp（17/2*b）A2+6800*a*exp（34/5*b）A2-10035*exp（3/2*b）-4‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎0149*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/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）A2+5100*a*exp（51/10*b）A2+4500*a*exp（9/2*b）A2+7200*a*exp（18/5*b）A2+3400*a*exp（17/5*b）A2+2600*a*exp（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）

用解二元非线性方程组的牛顿法的MATLAB^序求解线性方程组Jai=0，Jbi=0,得

a=b=

2.81100.5816

故所求的拟合曲线（7.13）为

f（X）=2.8110e^.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.2652.81198.27165.6059.1741.6625.9222.37

13.4712.8711.876.6914.872‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎4.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.A2;

Ew=max（fy）,

E仁sum（fy）/n,E2=sqrt（（sum（fy2））/n）,plot（xi,y,'r*'）,holdon

plot（x,F‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎,'b-'）,holdoff,

legend（'数据点（xi,yi）'，'拟合曲线y=f（x）'）

xlabel（'x'）,ylabel（'y'）,

title（'例7.3.1的数据点（xi,yi）和拟合曲线y=f（x）的图形'）

运行后屏幕显示数据（xi,yi）与拟合函数f的最大误差Ew=390.1415，平均误差

E1=36.9422和均方根误差‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎‎