数值计算方法上机实验题

数值计算方法上机实验

实验内容：

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=[0 10 20 30 40 50];y=[30.00 31.63 33.44 35.47 37.75 40.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;

>> I=Vsm('f',a,b,tol)

function y=f(x)

y=1./x;

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

>> 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,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

end

end

T2=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=S2;C1=C2;R1=R2;k=k+1; if (e<=tol)&(k>=5)

break

end

end

R2n=R2;

function y=f(x)

y=exp(x);

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

>> I=Vsm('f',a,b,tol)

function y=f(x)

y=1./x;

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

>> I=Vsm('f1',a,b,tol)