数值计算方法上机实验题

1、精选优质文档-倾情为你奉上精选优质文档-倾情为你奉上专心-专注-专业专心-专注-专业精选优质文档-倾情为你奉上专心-专注-专业数值计算方法上机实验实验内容：1从地面发射一枚火箭,在最初80 s 内记录起加速度如下表,试求火箭在第50s,80s时的速度.t(s)01020304050607080a(m/s)30.0031.6333.4435.4737.7540.3342.3946.6950.67要求:分别用复化梯形法,复化Simpson 法和 Romberg公式计算.2给定积分 和 ,分别用下列方法计算积分值要求准确到 ,并比较分析计算时间.1）变步长梯形法; 2）变步长 Simpson 法;

2、3) Romberg 方法.算法描述：1、复合梯形法： 输入 被积函数数据点t,a.输出 积分值复合Simpson法： 输入 被积函数f(x),积分区间a,b和n输出 复合Simpson积分值步1 步2 对执行步3 步4 输出Romberg积分法： 根据已知数据对其进行多项式拟合得出p(x);f(x)p(x);输入 被积函数f(x)，积分区间端点a,b,允许误差输出 Romberg积分值步1 步2 反复执行步3步9.步3 步4 反复执行步5步6.步5 步6 若xb,则退出本层循环.步7 执行步8 执行步9 若e且k5,则退出循环.步10 步11 输出2、变步长梯形算法：功能 求积分，允许误差为

3、。输入 被积函数f(x),a,b, 。输出 复合梯形积分值。步1 hb-a.步2 。步3 反复执行步4步10.步4 S0；xa+h/2.步5 反复执行步6步7.步6 SS+f(x);xx+h.步7 若xb,则退出本层循环。步8 .步9 .步10 若e,退出循环.步11 步12 输出复合Simpson法算法：功能 用复合Simpson公式求积分。输入 被积函数f(x),积分区间a,b和n。输出 复合Simpson积分值。步1 .步2 对k=1，2，n执行。步3 .步4 输出Romberg积分法算法：功能 计算积分，允许误差为。输入 被积函数f(x), 积分区间端点a,b, 允许误差为。输出 Ro

4、mberg积分值。步1 .步2 反复执行步3步9.步3 S0；xa+h/2.步4 反复执行步5步6.步5 SS+f(x);xx+h.步6 若xb,则退出本层循环。步7 .步8 执行.步10 若e且,则退出循环.步11 步12 输出源程序清单：1、复合梯形法源程序清单： t=0 10 20 30 40 50 60 70 80; a=30.00 31.63 33.44 35.47 37.75 40.33 42.39 46.69 50.67; h=10;v0=0; v50=v0+(h/2)*(a(1)+2*(a(2)+a(3)+a(4)+a(5)+a(6) v80=v0+(h/2)*(a(1)+2*

5、(a(2)+a(3)+a(4)+a(5)+a(6)+a(7)+a(8)+a(9)复合Simpson法源程序清单： t=0 10 20 30 40 50 60 70 80; a=30.00 31.63 33.44 35.47 37.75 40.33 42.39 46.69 50.67; h=20;v0=0; v80=v0+(h/6)*(a(1)+4*(a(2)+a(4)+a(6)+a(8)+2*(a(3)+a(5)+a(7)+a(9) x=0 10 20 30 40 50;y=30.00 31.63 33.44 35.47 37.75 40.33; p=polyfit(x,y,2);poly2s

6、ym(p);x=5 15 25 35 45 ;q=polyval(p,x); a=30.00 q(1) 31.63 q(2) 33.44 q(3) 35.47 q(4) 37.75 q(5) 40.33; h=10; v50=v0+(h/6)*(a(1)+4*(a(2)+a(4)+a(6)+a(8)+a(10)+2*(a(3)+a(5)+a(7)+a(9)+a(11)Romberg积分法源程序清单：function R2n=Romberg(f,a,b,tol)h=b-a;T1=(h/2)*(feval(f,a)+feval(f,b);S1=0;C1=0;R1=0;k=0;while 1 S=0

7、;x=a+h/2; while 1 S=S+feval(f,x);x=x+h; if x=b break end end T2=T1/2+(h/2)*S;S2=(4/3)*T2-(1/3)*T1; C2=(16/15)*S2-(1/15)*S1;R2=(64/63)*C2-(1/63)*C1; e=abs(R2-R1);h=h/2; T1=T2;S1=S2;C1=C2;R1=R2;k=k+1; if e=5 break endendR2n=R2; x=0 10 20 30 40 50 60 70 80; y=30.00 31.63 33.44 35.47 37.75 40.33 42.39 4

8、6.69 50.67; p=polyfit(x,y,3);y=poly2sym(p) function y=f(x)y=(26077*x3)/ - (50087*x2)/ + (12137*x)/ + 74033/2475 a=0;b=50;tol=0.; v50=Romberg(f,a,b,tol) a=0;b=80;tol=0.;80=Romberg(f,a,b,tol)2、变步长梯形法源程序清单：function T2n=Vsm(f,a,b,tol)h=b-a;T1=h/2*(feval(f,a)+feval(f,b);while 1 S=0;x=a+h/2; while 1 S=S+f

9、eval(f,x);x=x+h; if x=b break endendT2=T1/2+h*S/2;e=abs(T2-T1);h=h/2;T1=T2;if e a=1;b=3;tol=0.; I=Vsm(f,a,b,tol)function y=f(x)y=1./x; a=1;b=3;tol=0.; I=Vsm(f1,a,b,tol)变步长 Simpson 法源程序清单：function I=gauss2(gfun,a,b,n)h=(b-a)/n;gp1=1/sqrt(3);gp2=-gp1;k=1:n;g1=feval(gfun,a+h/2*(gp2+2*k-1);g2=feval(gfun

10、,a+h/2*(gp1+2*k-1);I=h/2*sum(g1+g2);function f=gfun(x)y=exp（x）; I=gauss2(gfun,1,3,8)function f=gfun(x)y=1./； I=gauss2(gfun,1,3,8)Romberg 方法源程序清单：function R2n=Vsm(f,a,b,tol)h=b-a;T1=h/2*(feval(f,a)+feval(f,b);S1=0;C1=0;R1=0;k=0;while 1 S=0;x=a+h/2; while 1 S=S+feval(f,x);x=x+h; if x=b break endendT2=T1/2+h*S/2;S2=4*T2/3-T1/3;C2=(16/15)*S2-(1/16)*S1;R2=(64/63)*C2-(1/63)*C1;e=abs(R2-R1);h=h/2;T1=T2;S1