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

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

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

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

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

例7.2.1给出一组数据点（xj,yi）列入表7-中，试用线性最小二乘法求拟合曲线,

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

表7-2例7.2.1的一组数据（xi,yi）

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（'实验数据（xi,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.126.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/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+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

ak

为求a1,a2,a3,a4使J达到最小，只需利用极值的必要条件一丄0（k1,2,3,4）,

运行后屏幕显示J分别对a1,a2,a3,a4的偏导数如下

Ja11=

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

解线性方程组Ja11=0，Ja21=0，Ja31=0，Ja41=0，输入下列程序

>>A=[56918107/10000,32097579/25000,

23667/250;32097579/25000,1377283/2500,

1377283/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.19056.4102-8.2574

f=716503695845759/140737488355328*xA3-7988544102557579/562949953421312*xA2+1804307491277693/281474976710656*x-4648521160813215/562949953421312

故所求的拟合曲线为

1.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）,

E1=sum（fy）/n,E2=sqrt（（sum（fy2））/n）plot（xi,y,'r*'）,holdon,plot（x,F,'b-'）,holdofflegend（'数据点（xi,yi）','拟合曲线y=f（x）'）,xlabel（'x'）,ylabel（'y'）,

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

运行后屏幕显示数据（x^yj与拟合函数f的最大误差平均误差E和均方根误差E

及其数据点（Xj,yj和拟合曲线y=f（x）的图形（略）.

Ew=E1=E2=

3.10540.90341.2409

7.3函数m（x）的选取及其MATLA程序

例7.3.1给出一组实验数据点（xi,yi）的横坐标向量为

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），试用线性最小二乘法求拟合曲线，并用（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）处的函数值，即输入程序

>>symsabx=[-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）,

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）-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/100）A2+（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/100）A2+（a*exp（18/5*b）-1211/50）A

2

为求a,b使J达到最小，只需利用极值的必要条件，得到关于a,b的线性方程组,

这可以由下面的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/100）A2+（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/100）A2+（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/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/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）-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/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）为

0.5816x

f（x）2.8110e.（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.8724.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,yj与拟合函数f的最大误差Ew=390.1415，平均误差

E1=36.9422和均方根误差巳=106.0317及其数据点（Xi,yj和拟合曲线y=f（x）的图

形（略）.

7.4多项式拟合及其MATLA程序

例7.4.1给出一组数据点（Xi,yi）列入表743中，试用线性最小二乘法求拟合曲线,并用（7.2）,（7.3）和（7.4）式估计其误差，作出拟合曲线.

表7-3例7.4.1的一组数据&,yi）

Xi

-2.9

■1.9-1.1-0.80

0.1

1.5

2.7

3.6

53.94

33.6820.8816.92

8.79

8.984.179.1219.88

解

（1）首先根据表7-3给出的数据点（Xi,yj，用下列MATLAB程序画出散点图

在MATLAB工作窗口输入程序

>>x=[-2.9-1.9-1.1-0.800.11.52.73.6];y=[53.9433.6820.8816.928.798.984.179.1219.88];

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

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

title（'例7.4.1的数据点（xi,yi）的散点图'）运行后屏幕显示数据的散点图（略）.

（3）用作线性最小二乘拟合的多项式拟合的MATLAB程序求待定系数ak（k1,2,3）.输入程序

>>a=polyfit（x,y,2）运行后输出（7.16）式的系数

a=

2.8302-7.37219.1382

故拟合多项式为

f（x）2.8302x27.3721x9.1382.

（4）编写下面的MATLA龍序估计其误差，并做出拟合曲线和数据的图形•输入程序

>>xi=[-2.9-1.9-1.1-0.800.11.52.73.6];y=[53.9433.6820.8816.928.798.984.179.1219.88];

n=length（xi）;f=2.8302,*xi.A2-7.3721.*xi+9.1382

x=-2.9:

0.001:

3.6;F=2.8302.*x.A2-7.3721.*x+8.79;

fy=abs（f-y）;fy2=fy.A2;Ew=max（fy）,E1=sum（fy）/n,E2=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）的图形'）

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

Ew=E1=E2=

0.7457,0.3892,0.4363

7.5拟合曲线的线性变换及其MATLA程序

例7.5.1给出一组实验数据点（xi,yi）的横坐标向量为x=（7.56.85.104.5

3.63.42.62.52.11.52.73.6），纵横坐标向量为y=（359.26165.6059.17

41.6625.9222.3713.4712.8711.876.6914.8724.22），试用线性变换和线性

最小二乘法求拟合曲线，并用（7.2），（7.3）和（7.4）式估计其误差，作出拟合曲线.

解

（1）首先根据给出的数据点（Xi，yj，用下列MATLAB^序画出散点图.

在MATLAB工作窗口输入程序

>>x=[7.56.85.104.53.63.42.62.52.11.52.7

3.6];

y=[359.26165.6059.1741.6625.9222.3713.4712.87

11.876.6914.8724.22];

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

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

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

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

（2）根据数据散点图，取拟合曲线为

bx

yaebx（a0,b0）,（7.19）

其中a,b是待定系数.令YIny,AIna,Bb，则（7.19）化为丫ABx.在

MATLAB工作窗口输入程序

>>x=[7.56.85.104.53.63.42.62.52.11.52.7

3.6];

y=[359.26165.6059.1741.6625.9222.3713.4712.8711.876.6914.8724.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）;fy2=fy.A2;Ew=max（fy）,E仁sum（fy）/n,

E2=sqrt（（sum（fy2））/n）

数据（xi,yi）与拟合函数f

运行后屏幕显示yaebx的系数b=0.6241,a=2.7039,

的最大误差Ew=67.6419，平均误差E1