曲线拟合的线性最小二乘法及其MATLAB程序.docx
《曲线拟合的线性最小二乘法及其MATLAB程序.docx》由会员分享，可在线阅读，更多相关《曲线拟合的线性最小二乘法及其MATLAB程序.docx（24页珍藏版）》请在冰豆网上搜索。
曲线拟合的线性最小二乘法及其MATLAB程序
1曲线拟合的线性最小二乘法及其MATLAB程序
例7.2.1给出一组数据点
列入表7–2中，试用线性最小二乘法求拟合曲线，并用（7.2），（7.3）和（7.4）式估计其误差，作出拟合曲线.
表7–2例7.2.1的一组数据
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.5068.04];
plot（x,y,'r*'）,
legend（'实验数据（xi,yi）'）
xlabel（'x'）,ylabel（'y'）,
title（'例7.2.1的数据点（xi,yi）的散点图'）
运行后屏幕显示数据的散点图（略）.
（3）编写下列MATLAB程序计算
在
处的函数值，即输入程序
>>symsa1a2a3a4
x=[-2.5-1.7-1.1-0.800.11.52.73.6];
fi=a1.*x.^3+a2.*x.^2+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.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,5832/125*a1+324/25*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*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
为求
使
达到最小，只需利用极值的必要条件
，得到关于
的线性方程组，这可以由下面的MATLAB程序完成，即输入程序
>>symsa1a2a3a4
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*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）;
Ja11=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/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,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/125,-14859/25];
C=B/A,f=poly2sym（C）
运行后屏幕显示拟合函数f及其系数C如下
C=5.0911-14.19056.4102-8.2574
f=716503695845759/140737488355328*x^3
-7988544102557579/562949953421312*x^2
+1804307491277693/281474976710656*x
-4648521160813215/562949953421312
故所求的拟合曲线为
.
（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.^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*'）,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）的图形'）
运行后屏幕显示数据
与拟合函数f的最大误差Ew，平均误差E1和均方根误差E2及其数据点
和拟合曲线y=f（x）的图形（略）.
Ew=E1=E2=
3.10540.90341.2409
7.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.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程序计算
在
处的函数值，即输入程序
>>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];
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）];
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*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/50）^2
为求
使
达到最小，只需利用极值的必要条件，得到关于
的线性方程组，这可以由下面的MATLAB程序完成，即输入程序
>>symsab
J=（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+（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/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）^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*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）
用解二元非线性方程组的牛顿法的MATLAB程序求解线性方程组Ja1=0，Jb1=0，得
a=b=
2.81100.5816
故所求的拟合曲线（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.2652.81198.27165.6059.1741.6625.9222.3713.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.^2;Ew=max（fy）,
E1=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）的图形'）
运行后屏幕显示数据
与拟合函数f的最大误差Ew=390.1415，平均误差E1=36.9422和均方根误差E2=106.0317及其数据点
和拟合曲线y=f（x）的图形（略）.
7.4多项式拟合及其MATLAB程序
例7.4.1给出一组数据点
列入表7–3中，试用线性最小二乘法求拟合曲线，并用（7.2），（7.3）和（7.4）式估计其误差，作出拟合曲线.
表7–3例7.4.1的一组数据
xi
-2.9-1.9-1.1-0.800.11.52.73.6
yi
53.9433.6820.8816.928.798.984.179.1219.88
解
（1）首先根据表7–3给出的数据点
，用下列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程序求待定系数
.输入程序
>>a=polyfit（x,y,2）
运行后输出（7.16）式的系数
a=
2.8302-7.37219.1382
故拟合多项式为
.
（4）编写下面的MATLAB程序估计其误差，并做出拟合曲线和数据的图形.输入程序
>>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.^2-7.3721.*xi+9.1382
x=-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,
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）的图形'）
运行后屏幕显示数据
与拟合函数f的最大误差Ew，平均误差E1和均方根误差E2及其数据点（xi,yi）和拟合曲线y=f（x）的图形（略）.
Ew=E1=E2=
0.7457,0.3892,0.4363
7.5拟合曲线的线性变换及其MATLAB程序
例7.5.1给出一组实验数据点
的横坐标向量为x=（7.56.85.104.53.63.42.62.52.11.52.73.6），纵横坐标向量为y=（359.26165.6059.1741.6625.9222.3713.4712.8711.876.6914.8724.22），试用线性变换和线性最小二乘法求拟合曲线，并用（7.2），（7.3）和（7.4）式估计其误差，作出拟合曲线.
解
（1）首先根据给出的数据点
，用下列MATLAB程序画出散点图.
在MATLAB工作窗口输入程序
>>x=[7.56.85.104.53.63.42.62.52.11.52.73.6];
y=[359.26165.6059.1741.6625.9222.3713.4712.8711.876.6914.8724.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.56.85.104.53.63.42.62.52.11.52.73.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.
升级会员 