1、精選優(yōu)質(zhì)文檔-傾情為你奉上精選優(yōu)質(zhì)文檔-傾情為你奉上專心-專注-專業(yè)專心-專注-專業(yè)精選優(yōu)質(zhì)文檔-傾情為你奉上專心-專注-專業(yè)數(shù)值計算方法上機實驗實驗內(nèi)容：1從地面發(fā)射一枚火箭,在最初80 s 內(nèi)記錄起加速度如下表,試求火箭在第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、復合梯形法： 輸入 被積函數(shù)數(shù)據(jù)點t,a.輸出 積分值復合Simpson法： 輸入 被積函數(shù)f(x),積分區(qū)間a,b和n輸出 復合Simpson積分值步1 步2 對執(zhí)行步3 步4 輸出Romberg積分法： 根據(jù)已知數(shù)據(jù)對其進行多項式擬合得出p(x);f(x)p(x);輸入 被積函數(shù)f(x)，積分區(qū)間端點a,b,允許誤差輸出 Romberg積分值步1 步2 反復執(zhí)行步3步9.步3 步4 反復執(zhí)行步5步6.步5 步6 若xb,則退出本層循環(huán).步7 執(zhí)行步8 執(zhí)行步9 若e且k5,則退出循環(huán).步10 步11 輸出2、變步長梯形算法：功能 求積分，允許誤差為

3、。輸入 被積函數(shù)f(x),a,b, 。輸出 復合梯形積分值。步1 hb-a.步2 。步3 反復執(zhí)行步4步10.步4 S0；xa+h/2.步5 反復執(zhí)行步6步7.步6 SS+f(x);xx+h.步7 若xb,則退出本層循環(huán)。步8 .步9 .步10 若e,退出循環(huán).步11 步12 輸出復合Simpson法算法：功能 用復合Simpson公式求積分。輸入 被積函數(shù)f(x),積分區(qū)間a,b和n。輸出 復合Simpson積分值。步1 .步2 對k=1，2，n執(zhí)行。步3 .步4 輸出Romberg積分法算法：功能 計算積分，允許誤差為。輸入 被積函數(shù)f(x), 積分區(qū)間端點a,b, 允許誤差為。輸出 Ro

4、mberg積分值。步1 .步2 反復執(zhí)行步3步9.步3 S0；xa+h/2.步4 反復執(zhí)行步5步6.步5 SS+f(x);xx+h.步6 若xb,則退出本層循環(huán)。步7 .步8 執(zhí)行.步10 若e且,則退出循環(huán).步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