用Romberg方法求解积分.docx

用Romberg方法求解积分.docx

《用Romberg方法求解积分.docx》由会员分享，可在线阅读，更多相关《用Romberg方法求解积分.docx（17页珍藏版）》请在冰豆网上搜索。

用Romberg方法求解积分

1.用Romberg方法求解积分

，要求误差不超过

解：

Romberg.m文件：

function[I,step]=Romberg（f,a,b,EPS）

%Romberg.m是用龙贝格公式求积分

%f为被积函数

%EPS为积分结果精度

%a,b为积分区间的上下限

%I为积分结果；step为积分的子区间数

m=1

k=0

Er=0.1

H=b-a

S=zeros（1,1）

S（1,1）=（H/2）*（subs（sym（f）,findsym（sym（f））,a）+subs（sym（f）,findsym（sym（f））,b））

whileEr>EPS

k=k+1

f1=0

H=H/2

fori=1:

m

x=a+H*（2*i-1）

f1=f1+subs（sym（f）,findsym（sym（f））,x）

end

S（k+1,1）=S（k,1）/2+H*f1

m=2*m

forn=1:

k

S（k+1,n+1）=S（k+1,n）+（S（k+1,n）-S（k,n））/（4^n-1）

end

Er=abs（S（k+1,n+1）-S（k,n））

end

I=S（k+1,k+1）

step=k

命令：

clear

clc

formatshort

a=0;b=0.8;EPS=1e-2;

[I,step]=Romberg（'x^（1/2）',a,b,EPS）

计算结果：

m=

1

k=

0

Er=

0.1000

H=

0.8000

S=

0

S=

0.3578

k=

1

f1=

0

H=

0.4000

x=

0.4000

f1=

0.6325

S=

0.3578

0.4319

m=

2

S=

0.35780

0.43190.4566

Er=

0.0988

k=

2

f1=

0

H=

0.2000

x=

0.2000

f1=

0.4472

x=

0.6000

f1=

1.2218

S=

0.35780

0.43190.4566

0.46030

m=

4

S=

0.35780

0.43190.4566

0.46030.4698

S=

0.357800

0.43190.45660

0.46030.46980.4707

Er=

0.0141

k=

3

f1=

0

H=

0.1000

x=

0.1000

f1=

0.3162

x=

0.3000

f1=

0.8640

x=

0.5000

f1=

1.5711

x=

0.7000

f1=

2.4077

S=

0.357800

0.43190.45660

0.46030.46980.4707

0.470900

m=

8

S=

0.357800

0.43190.45660

0.46030.46980.4707

0.47090.47450

S=

0.357800

0.43190.45660

0.46030.46980.4707

0.47090.47450.4748

S=

0.3578000

0.43190.456600

0.46030.46980.47070

0.47090.47450.47480.4748

Er=

0.0042

I=

0.4748

step=

3

I=

0.4748

step=

3

2.设方程组

试用Jacobi迭代法求解此方程，

，当

时终止迭代。

解：

Jacobi.m文件：

functionJacobi（A,b,max,eps）%max为最大迭代次数，eps为容许误差

n=length（A）;x=zeros（n,1）;x1=zeros（n,1）;k=0;

while1

x1

（1）=（b

（1）-A（1,2:

n）*x（2:

n,1））/A（1,1）

fori=2:

n-1

x1（i）=（b（i）-A（i,1:

i-1）*x（1:

i-1,1）-A（i,i+1:

n）*x（i+1:

n,1））/A（i,i）

end

x1（n）=（b（n）-A（n,1:

n-1）*x（1:

n-1,1））/A（n,n）

k=k+1

ifsum（abs（x1-x））

fprintf（'number=%d\n',k）

break

end

ifk>=max

fprintf（'TheMethodisdisconvergent\n'）

break

end

x=x1

end

ifk

fori=1:

n

fprintf（'x[%d]=%f\n',i,x1（i））

end

end

命令：

clear

clc

formatshort

A=[521;-142;2-310];

b=[-12203]';

max=100;

eps=1e-5

Jacobi（A,b,max,eps）

计算结果：

i=

1

A=

521

-142

2-310

b=

-12

20

3

D=

500

040

0010

L=

000

100

-230

U=

0-2-1

00-2

000

D0=

0.200000

00.25000

000.1000

x0=

0

0

0

B=

0-0.4000-0.2000

0.25000-0.5000

-0.20000.30000

f=

-2.4000

5.0000

0.3000

x=

-2.4000

5.0000

0.3000

x0=

-2.4000

5.0000

0.3000

i=

2

x=

-4.4600

4.2500

2.2800

x0=

-4.4600

4.2500

2.2800

i=

3

x=

-4.5560

2.7450

2.4670

x0=

-4.5560

2.7450

2.4670

i=

4

x=

-3.9914

2.6275

2.0347

x0=

-3.9914

2.6275

2.0347

i=

5

x=

-3.8579

2.9848

1.8865

x0=

-3.8579

2.9848

1.8865

i=

6

x=

-3.9712

3.0922

1.9670

x0=

-3.9712

3.0922

1.9670

i=

7

x=

-4.0303

3.0237

2.0219

x0=

-4.0303

3.0237

2.0219

i=

8

x=

-4.0139

2.9815

2.0132

x0=

-4.0139

2.9815

2.0132

i=

9

x=

-3.9952

2.9900

1.9972

x0=

-3.9952

2.9900

1.9972

i=

10

x=

-3.9954

3.0026

1.9960

x0=

-3.9954

3.0026

1.9960

i=

11

x=

-4.0002

3.0031

1.9999

x0=

-4.0002

3.0031

1.9999

i=

12

x=

-4.0012

3.0000

2.0010

x0=

-4.0012

3.0000

2.0010

i=

13

x=

-4.0002

2.9992

2.0002

x0=

-4.0002

2.9992

2.0002

i=

14

x=

-3.9997

2.9998

1.9998

x0=

-3.9997

2.9998

1.9998

i=

15

x=

-3.9999

3.0002

1.9999

x0=

-3.9999

3.0002

1.9999

i=

16

x=

-4.0000

3.0001

2.0000

x0=

-4.0000

3.0001

2.0000

i=

17

x=

-4.0000

3.0000

2.0000

x0=

-4.0000

3.0000

2.0000

i=

18

x=

-4.0000

3.0000

2.0000

x0=

-4.0000

3.0000

2.0000

i=

19

x=

-4.0000

3.0000

2.0000

x0=

-4.0000

3.0000

2.0000

i=

20

x=

-4.0000

3.0000

2.000

x=

-4.0000

3.0000

2.0000

i=

20