matlab拟合最小二乘法

最小二乘法拟合

一、 绘制离散点图

>> x=[0 5 10 15 20 25 30 35 40 45 50 55];

y=[0 1.27*10^-4 2.16*10^-4 2.86*10^-4 3.44*10^-4 3.87*10^-4 4.15*10^-4 4.37*10^-4

4.51*10^-4 4.58*10^-4 4.02*10^-4 4.64*10^-4];

plot(x,y,'r*');

legend('实验数据(xi,yi)');

xlabel('x'), ylabel('y'), title('数据点(xi,yi)的散点图')

运行如图

二、编写程序计算在x 处的函数值

>> syms a1 a2 a3 a4

x=[0 5 10 15 20 25 30 35 40 45 50 55];

fi=a1.*x.^3+ a2.*x.^2+ a3.*x+ a4

运行后显示a1、a2、a3、a4的方程组

fi =

[ a4, 125*a1 + 25*a2 + 5*a3 + a4, 1000*a1 + 100*a2 + 10*a3 + a4, 3375*a1 + 225*a2 + 15*a3 + a4, 8000*a1 + 400*a2 + 20*a3 + a4, 15625*a1 + 625*a2 + 25*a3 + a4, 27000*a1 + 900*a2 +

30*a3 + a4, 42875*a1 + 1225*a2 + 35*a3 + a4, 64000*a1 + 1600*a2 + 40*a3 + a4, 91125*a1 +

2025*a2 + 45*a3 + a4, 125000*a1 + 2500*a2 + 50*a3 + a4, 166375*a1 + 3025*a2 + 55*a3 + a4]

三、编写构造平方和的matlab程序

>> y=[0 1.27*10^-4 2.16*10^-4 2.86*10^-4 3.44*10^-4 3.87*10^-4 4.15*10^-4 4.37*10^-4

4.51*10^-4 4.58*10^-4 4.02*10^-4 4.64*10^-4];

fi=[ a4, 125*a1 + 25*a2 + 5*a3 + a4, 1000*a1 + 100*a2 + 10*a3 + a4, 3375*a1 + 225*a2 + 15*a3 + a4, 8000*a1 + 400*a2 + 20*a3 + a4, 15625*a1 + 625*a2 + 25*a3 + a4, 27000*a1 + 900*a2 + 30*a3 + a4, 42875*a1 + 1225*a2 + 35*a3 + a4, 64000*a1 + 1600*a2 + 40*a3 + a4, 91125*a1 + 2025*a2 + 45*a3 + a4, 125000*a1 + 2500*a2 + 50*a3 + a4, 166375*a1 + 3025*a2 + 55*a3 + a4] fy=fi-y;

fy2=fy.^2;

J=sum(fy.^2)

运行的如下结果

fi =

[ a4, 125*a1 + 25*a2 + 5*a3 + a4, 1000*a1 + 100*a2 + 10*a3 + a4, 3375*a1 + 225*a2 + 15*a3 + a4, 8000*a1 + 400*a2 + 20*a3 + a4, 15625*a1 + 625*a2 + 25*a3 + a4, 27000*a1 + 900*a2 + 30*a3 + a4, 42875*a1 + 1225*a2 + 35*a3 + a4, 64000*a1 + 1600*a2 + 40*a3 + a4, 91125*a1 + 2025*a2 + 45*a3 + a4, 125000*a1 + 2500*a2 + 50*a3 + a4, 166375*a1 + 3025*a2 + 55*a3 + a4]

J =

58593218750*a1^2 + 2386725000*a1*a2 + 49967500*a1*a3 + 1089000*a1*a4 - (274377591928296252150123*a1)/576460752303423488000 + 24983750*a2^2 + 1089000*a2*a3 + 25300*a2*a4 - (31331074233255294718193*a2)/2882303761517117440000 + 12650*a3^2 + 660*a3*a4 - (7819978335372091569501*a3)/28823037615171174400000 + 12*a4^2 - (2298349019433749545307*a4)/288230376151711744000000 + 520374483464852566590953249225508026224249/332306998946228968225951765070086144000000000000

四、求a1、a2、a3、a4使J达到最小，分别对a1、a2、a3、a4求偏导数，使之等于0 程序如下

syms a1 a2 a3 a4

J=58593218750*a1^2+2386725000*a1*a2+49967500*a1*a3+1089000*a1*a4-(274377591928296252150123*a1)/576460752303423488000+24983750*a2^2+1089000*a2*a3 + 25300*a2*a4 - (31331074233255294718193*a2)/2882303761517117440000+12650*a3^2+660*a3*a4-(7819978335372091569501*a3)/28823037615171174400000+12*a4^2-(2298349019433749545307*a4)/288230376151711744000000+520374483464852566590953249225508026224249/332306998946228968225951765070086144000000000000

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)

运行得

Ja11 =

117186437500*a1 + 2386725000*a2 + 49967500*a3 + 1089000*a4 - 1903877/4000

Ja21 =

2386725000*a1 + 49967500*a2 + 1089000*a3 + 25300*a4 - 217403/20000

Ja31 =

49967500*a1 + 1089000*a2 + 25300*a3 + 660*a4 - 27131/100000

Ja41 =

1089000*a1 + 25300*a2 + 660*a3 + 24*a4 - 3987/500000

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

>> A=[117186437500, 2386725000,49967500,1089000;

2386725000, 49967500, 1089000,25300;

49967500,1089000,25300, 660;

1089000, 25300, 660,24];

B=[274377591928296252150123/576460752303423488000,31331074233255294718193/2882303761517117440000,7819978335372091569501/28823037615171174400000,

2298349019433749545307/288230376151711744000000];

C=B/A, f=poly2sym(C)

运行即可得

C =

1.0e-004 *

0.0000 -0.0052 0.2634 0.0178

f =

(4154371348005165*x^3)/1208925819614629174706176 - (307872452451165*x^2)/590295810358705651712 + (971769044862535*x)/36893488147419103232 + 4212074133985165/2361183241434822606848

五、编写下面的MATLAB程序估计其误差，并作出拟合曲线和数据的图形.输入程序 >> xi=[0 5 10 15 20 25 30 35 40 45 50 55];

y=[0 1.27*10^-4 2.16*10^-4 2.86*10^-4 3.44*10^-4 3.87*10^-4 4.15*10^-4 4.37*10^-4

4.51*10^-4 4.58*10^-4 4.02*10^-4 4.64*10^-4];

n=length(xi);

f=4154371348005165./1208925819614629174706176.*xi.^3 307872452451165./590295810358705651712.*xi.^2

971769044862535./36893488147419103232.*xi

4212074133985165./2361183241434822606848.;

x=0:1: 55;

F=4154371348005165./1208925819614629174706176.*x.^3 307872452451165./590295810358705651712.*x.^2

971769044862535./36893488147419103232.*x

4212074133985165./2361183241434822606848.;

fy=abs(f-y);

fy2=fy.^2;

Ew=max(fy), E1=sum(fy)/n, E2=sqrt((sum(fy2))/n)

plot(xi,y,'r*'), hold , plot(x,F,'b-'), hold ;

legend('数据点(xi,yi)','拟合曲线y=f(x)');

xlabel('x'), ylabel('y'), title('数据点(xi,yi)和拟合曲线y=f(x)的图形') 如图 - + + - + +

最大误差Ew =4.2438e-005

平均误差E1 =9.1046e-006

均方根误差

E2=1.4672e-005