matlab数值分析程序--高等数学,数值代数的matlab实现-文字版, matlab电子书, 和matlab 有关的电子书:

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

7.2  曲线拟合的线性最小二乘法及其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.8     0      0.1    1.5    2.7    3.6

yi

-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 -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 a4

x=[-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 =[ -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.^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程序完成,即输入程序

>> syms a1 a2 a3 a4

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.1905    6.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.8  0  0.1  1.5  2.7  3.6];

    y=[-192.9 -85.50 -36.15 -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 off

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.105 4         0.903 4           1.240 9

欢迎转载,转载请注明来自一手册:http://yishouce.com/book/3/3070200.html
友情链接It题库(ittiku.com)| 版权归yishouce.com所有| 友链等可联系 admin#yishouce.com|粤ICP备16001685号-1