用Romberg方法求解积分.doc

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))
while Er > EPS
k = k + 1
f1 = 0
H = H/2
for i = 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
for n = 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
format short
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.3578 0
0.4319 0.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.3578 0
0.4319 0.4566
0.4603 0

m =

4

S =

0.3578 0
0.4319 0.4566
0.4603 0.4698

S =

0.3578 0 0
0.4319 0.4566 0
0.4603 0.4698 0.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.3578 0 0
0.4319 0.4566 0
0.4603 0.4698 0.4707
0.4709 0 0

m =

8

S =

0.3578 0 0
0.4319 0.4566 0
0.4603 0.4698 0.4707
0.4709 0.4745 0

S =

0.3578 0 0
0.4319 0.4566 0
0.4603 0.4698 0.4707
0.4709 0.4745 0.4748

S =

0.3578 0 0 0
0.4319 0.4566 0 0
0.4603 0.4698 0.4707 0
0.4709 0.4745 0.4748 0.4748

Er =

0.0042

I =

0.4748

step =

3

I =

0.4748

step =

3

2. 设方程组试用Jacobi迭代法求解此方程，，当时终止迭代。
解：Jacobi.m文件：
function Jacobi(A, b, max, eps) %max为最大迭代次数，eps为容许误差
n = length(A); x = zeros(n, 1); x1 = zeros(n, 1); k = 0;
while 1
x1(1) = ( b(1) - A(1,2:n) * x(2:n,1) )/A(1,1)
for i = 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
if sum(abs(x1 - x)) < eps
fprintf('number = %d\n',k)
break
end
if k >= max
fprintf('The Method is disconvergent\n')
break
end
x = x1
end
if k < max
for i = 1:n
fprintf( 'x[ %d ] = %f\n',i,x1(i) )
end
end

命令：
clear
clc
format short
A = [5 2 1; -1 4 2; 2 -3 10];
b = [-12 20 3]';
max = 100;
eps = 1e-5
Jacobi(A, b, max, eps)

计算结果：
i =
1
A =
5 2 1
-1 4 2
2 -3 10
b =
-12
20
3
D =
5 0 0
0 4 0
0 0 10
L =
0 0 0
1 0 0
-2 3 0
U =
0 -2 -1
0 0 -2
0 0 0
D0 =
0.2000 0 0
0 0.2500 0
0 0 0.1000
x0 =
0
0
0
B =
0 -0.4000 -0.2000
0.2500 0 -0.5000
-0.2000 0.3000 0
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