曲线拟合的线性最小二乘法及其MATLAB程序47476Word文档格式.docx

曲线拟合的线性最小二乘法及其MATLAB程序47476Word文档格式.docx

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

title（'

2.1的数据点（xi，yi）的散点图'

）

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

（3）编写下列MATLAB程序计算

在

处的函数值,即输入程序

symsa1a2a3a4

x=［—2.5—1。

7-1。

1—0.800.11.52.73。

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。

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。

6］；

15-26。

52-9。

10—8。

43—13。

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。

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）’）,

），ylabel（'

y'

）,

title（’例7。

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）式估计其误差,作出拟合曲线。

〉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）的散点图'

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

处的函数值，即输入程序

〉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。

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。

66，25。

92,22。

37，13.47，12.87，11.87，6。

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=

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）;

^2;

Ew=max（fy），

E1=sum（fy）/n,E2=sqrt（（sum（fy2））/n），plot（xi,y,’r*'

），holdon

plot（x，F，’b—'

）,holdoff,

数据点（xi，yi）'

，’拟合曲线y=f（x）’）

xlabel（’x’），ylabel（’y'

），

title（’例7.3。

与拟合函数f的最大误差Ew=390。

1415，平均误差E1=36.9422和均方根误差E2=106。

0317及其数据点

和拟合曲线y=f（x）的图形（略）。

4多项式拟合及其MATLAB程序

4。

列入表7–3中，试用线性最小二乘法求拟合曲线,并用（7。

3）和（7。

表7–3例7.4.1的一组数据

-2.9—1。

9-1.1-0.800.11。

52.73.6

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。

y=[53。

9433.6820。

928。

798。

984.179.1219。

88］;

plot（x,y，’r＊’），legend（'

数据点（xi，yi）’）

xlabel（’x'

），ylabel（’y’），

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。

y=［53.9433.6820。

8816.928。

88]；

f=2.8302。

＊xi.^2—7。

3721。

＊xi+9.1382

9:

001：

3。

6;

8302。

*x.^2—7。

＊x+8。

79；

fy=abs（f-y）;

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）’）

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

5拟合曲线的线性变换及其MATLAB程序

例7.5。

的横坐标向量为x=（7。

56。

85。

104。

53.63.42.62。

52.11.52。

73.6），纵横坐标向量为y=（359。

26165。

6059。

1741.6625。

4712。

8711。

22），试用线性变换和线性最小二乘法求拟合曲线，并用（7。

2），（7.3）和（7.4）式估计其误差，作出拟合曲线。

解

（1）首先根据给出的数据点

，用下列MATLAB程序画出散点图。

〉〉x=[7.56.85.104.53.63。

42.62。

11.52.73.6］；

y=［359.26165.6059.1741。

9222.3713。

4712.8711.876。

plot（x，y,'

），legend（’数据点（xi,yi）'

），ylabel（’y'

title（’例7.5。

1的数据点（xi,yi）的散点图'

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

，（7。

19）

其中

是待定系数.令

，则（7。

19）化为

.在MATLAB工作窗口输入程序

〉〉x=［7.56。

85.104.53。

63.42.62.52。

y=［359.26165。

9222.3713.4712.8711。

Y=log（y）；

a=polyfit（x，Y,1）;

B=a

（1）;

A=a

（2）；

b=B,a=exp（A）

n=lengt