用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
ifkfori=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