数值分析实验和程序

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

y=[0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.02 4.64].*1e-4 plot(t,y,'r*'),

legend('实验数据(t,yi)') xlabel('x'), ylabel('y'),

title('数据点(t,yi)的散点图')，显示结果为：

2.求参数a1、a2、a3的解析表达式

计算f(x)在(xi,yi)处的函数值，即输入程序 syms a1 a2 a3

t=[0 5 10 15 20 25 30 35 40 45 50 55]; fi=a1.*t+ a2.*t.^2+ a3.*t.^3

运行后屏幕显示关于a1,a2, a3的线性方程组：

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

构造误差平方和：

y=[0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.02 4.64].*1e-4;

fy=fi-y; fy2=fy.^2; J=sum(fy.^2); 运行后屏幕显示误差平方和如下:

J=12650*a1^2+1089000*a1*a2+49967500*a1*a3- (7819978335372091569501*a1)/28823037615171174400000

+

24983750*a2^2

+2386725000*a2*a3-(31331074233255294718193*a2)/2882303761517117440000+58593218750*a3^2-(274377591928296252150123*a3)/576460752303423488000+520374483464852566590953249225508026224249/332306998946228968225951765070086144000000000000

?J?0?a为求a1,a2, a3使J达到最小，只需利用极值的必要条件k(k=1,2,3)，

得到关于a1,a2, a3的线性方程组，这可以由下面的MATLAB程序完成，即输入程序

Ja1=diff(J,a1); Ja2=diff(J,a2); Ja3=diff(J,a3); Ja11=simple(Ja1), Ja21=simple(Ja2), Ja31=simple(Ja3), 运行后屏幕显示J分别对a1, a2 ,a3的偏导数如下

Ja11 =25300*a1 + 1089000*a2 + 49967500*a3 - 27131/100000 Ja21 =1089000*a1 + 49967500*a2 + 2386725000*a3 - 217403/20000 Ja31 =49967500*a1 + 2386725000*a2 + 117186437500*a3 - 1903877/4000解线性方程组Ja11 =0，Ja21 =0，Ja31 =0，Ja41 =0，输入下列程序

A=[25300, 1089000, 49967500; 1089000, 49967500, 2386725000; 49967500, 2386725000, 117186437500];

B=[27131/100000, 217403/20000, 1903877/4000]; C=B/A

运行后求得系数C如下: C =

2.6569e-005 -5.2948e-007 3.5168e-009 故所求的拟合曲线为

f(x)=2.6569×10-5t-5.2948×10-7t2+3.5168×10-9t3

3．估计其误差，并作出拟合曲线和数据的图形

输入程序：

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

y=[0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.02 4.64].*1e-4; n=length(t);

f=2.6569*1e-5*t-5.2948*1e-7*t.^2+3.5168*1e-9*t.^3; tx=0:1:100;

F=2.6569*1e-5*tx-5.2948*1e-7*tx.^2+3.5168*1e-9*tx.^3; fy=abs(f-y); fy2=fy.^2; Ew=max(fy), E1=sum(fy)/n, E2=sqrt((sum(fy2))/n)

plot(t,y,'r*'), hold on, plot(tx,F,'b-'), hold off legend('数据点(t,y)','拟合曲线y=f(x)'), xlabel('x'), ylabel('y'),

title('数据点(t,y)和拟合曲线y=f(x)的图形')

运行后屏幕显示数据(t,y)与拟合函数f的最大误差Ew，平均误差E1和均方根误差E2及其数据点(xi,yi)和拟合曲线y=f(x)的图形

结果为：

最大误差：Ew = 4.2350e-005 平均误差：E1 = 9.0785e-006 均方根误差：E2 = 1.4684e-005

4.用软件自带函数进行三次多项式拟合，比较拟合效果

编写程序如下：

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

y=[0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.02 4.64].*1e-4; tx=0:1:100;

F=2.6569*1e-5*tx-5.2948*1e-7*tx.^2+3.5168*1e-9*tx.^3; plot(t,y,'r*'), hold on, plot(tx,F,'b-') F1=polyfit(t,y,3)

F1=F1(1)*tx.^3+F1(2)*tx.^2+F1(3)*tx+F1(4); plot(tx,F1,'+'),hold off

legend('数据点(t,y)','拟合曲线y=f(x)','拟合曲线y=f1(x)'), xlabel('x'), ylabel('y'),

title('数据点(t,y)、拟合曲线y=f(x)、y=f1(x)图形') 结果如下图：

用软件拟合的三次多项式参数为