数值计算方法上机实验题

数值计算方法上机实验题

数值计算方法上机实验

实验内容：

1

.

要求:分别用复化梯形法,复化Simpson 法和 Romberg 公式计算.

2．给定积分

dx e x

3

1

和

dx x ?3

1

1 ,分别用下列方法计算积分值要求准确到510- ,并比较分析计算时间. 1）变步长梯形法; 2）变步长 Simpson 法; 3) Romberg 方法.

算法描述：

1、复合梯形法：?=t

dt t a t V 0)()( ))()(2)((21

1

∑-=++=n k k n b f x f a f h

T

输入 被积函数数据点t,a. 输出 积分值.n T

复合Simpson 法：?

=t

dt t a t V 0)()( ))()(2)(4)((6101

12

1∑∑---=++++=n k n k k k n b f x f x f a f h S

输入 被积函数f(x),积分区间[a,b]和n 输出 复合Simpson 积分值n S

步1 .);()(;a x b f a f S n

a

b h n ?-?-? 步2 对n k ,,2,1 =执行).(2;2

);(4;2x f S S h

x x x f S S h x x n n n n +?+?+?+?

步3 n n S h

S ??6

步4 输出n S

Romberg 积分法：

根据已知数据对其进行多项式拟合得出p(x);f(x)?p(x); 输入 被积函数f(x)，积分区间端点a,b,允许误差ε 输出 Romberg 积分值n R 2 步1 .0;0;0;0));()((2;1111?===+?

-?k R C S b f a f h

T a b h 步2 反复执行步3→步9. 步3 .2

;0h a x S +

步4 反复执行步5→步6. 步5 ;);(h x x x f S S +?+?

步6 若x ≥b,则退出本层循环. 步7 执行.63

1

6364;1511516;3134;2212212212212C C R S S C T T S S h T T

-?-?-?+?

步8 执行.1;;;;;2

;2121212112+

-?k k R R C C S S T T h

h R R e 步9 若e ≤ε且k ≥5,则退出循环. 步10 .22R R n ? 步11

输出.2n R

2、变步长梯形算法：

功能 求积分

b

a

)(dx x f ，允许误差为ε。

输入 被积函数f(x),a,b, ε。 输出 复合梯形积分值T n 2。 步1

h ?b-a. 步2

T 1?

))()((2

b f a f h

-。 步3 反复执行步4——步10. 步4 S ?0；x ?a+h/2. 步5 反复执行步6——步7. 步6 S ?S+f(x);x ?x+h.

步7 若x ≥b,则退出本层循环。 步8

S T T

+?2

h

212

.

步9 T T T T h h e 2

112;2

;??-?. 步10 若e ε≤,退出循环. 步11

.22T T

n

步12 输出T n 2

复合Simpson 法算法：

功能 用复合Simpson 公式求积分

b

a

)(dx x f 。 输入 被积函数f(x),积分区间[a,b]和n 。 输出 复合Simpson 积分值S n 。 步1 .);()(;h a x b f a f n a

b S n

-?-?

. 步2 对k=1，2，…，n 执行)(2;2

);(4;

2

x x f h x x x f h x S S S S n n n n +?+?+?+?。

步3

S S n n h

2

. 步4 输出S n

Romberg 积分法算法：

功能 计算积分

b

a

)(dx x f ，允许误差为ε。

输入 被积函数f(x), 积分区间端点a,b, 允许误差为ε。 输出

Romberg 积分值R n 2。

步1 0;0;0;0));()((2

;

h 1111?===+?-?k b f a f h a b R C S T .

步2 反复执行步3——步9. 步3 S ?0；x ?a+h/2. 步4 反复执行步5——步6. 步5 S ?S+f(x);x ?x+h.

步6 若x ≥b,则退出本层循环。 步7

C C R S S C T T S T T

S 1

2212212212 63

16364;1511516;31342h 2-?-?-??+?；. 步8 执行1k k ;;;;;2;2

1212121

12

+

-?

R

R C C S S T T R R

h

h e .

步10 若e ε≤且5≥k ,则退出循环. 步11

.22R R

n

步12 输出R n 2

源程序清单：

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*(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=************];y=[30.0031.6333.4435.4737.7540.33]; >> p=polyfit(x,y,2);poly2sym(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;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<=tol&k>=5

break

end

end

R2n=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 46.69 50.67];

>> p=polyfit(x,y,3);y=poly2sym(p) function y=f(x)

y=(5930286613326077*x^3)/295147905179352825856 -

(3295392375650087*x^2)/4611686018427387904 +

(3387527922312137*x)/18014398509481984 + 74033/2475

>> a=0;b=50;tol=0.000005;

>> v50=Romberg('f',a,b,tol)

>> a=0;b=80;tol=0.000005;

>>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+feval(f,x);x=x+h;

if x>=b

break

end

end

T2=T1/2+h*S/2;

e=abs(T2-T1);h=h/2;T1=T2;

if e<=tol

break

end

end

T2n=T2;

function y=f(x)

y=exp(x);

>> a=1;b=3;tol=0.000005;