第1页 / 共8页
第2页 / 共8页
第3页 / 共8页
第4页 / 共8页
第5页 / 共8页
第6页 / 共8页
第7页 / 共8页
第8页 / 共8页
复合梯形及复合辛普森求积计算积分、龙贝格求积.docx
计算实习题
计算方法大作业
计算实习题
1. 用不同数值方法计算积分 1
0
x
ln d
x x
4
9
(1) 取不同的步长 h. 分别用复合梯形及复合辛普森求积计算积分, 给出误差中关于 h 的函数,
并与积分精确值比较两个公式的精度, 是否存在一个最小的 h, 使得精度不能再被改善?
(2) 用龙贝格求积计算完成问题(1).
(3) 用自适应辛普森积分, 使其精度达到 10−4.
解：(1)根据题意，可令 (
f x
)
=
x
ln
x
，a=0, b=1, 且有
lim ln
x
x
0
x
0
( )
f
f a
( )
f
f b
( )
f x
0
lim ( )
(0)
f x
x
(1) 0
ln
x
4
x x
a
, 等分节点
n
将积分区间[a,b]划分成 n 等分, 步长 b
h
kx
a
kh
k
，
0,1,
,
n
.
复合梯形公式为
T
n
h
2
其余项为
n
1
k
0
f x
k
f x
k
1
h
2
( ) 2
f a
n
1
k
1
f x
k
( )
f b
,
(
R f
n
)
I T
n
ha
2
b
12
f
)
(
,
,a b
.
复合辛普森公式为
f x
k
S
1
n
n
h
6
k
0
4
f x
k
1/2
+
f x
k
1
h
6
( ) 4
f a
n
1
k
0
f
x
k
+2
n
1
k
1
f x
k
( )
f b
,
其余项为
(
R f
n
)
I
S
n
4
h
h
180 2
(4)
( )
f
,a b
,
.
(1)
(2)
(3)
(4)
根据(1)式和(3)式, 使用 MATLAB 编程计算, 程序分别如附录 1 和 2. 两种算法都是给定最
大等分数 nmax, 并计算 n=1~nmax 对应的复合梯形积分值 T(n)和复合梯形积分值 S(n).
若最终结果(积分值)只需保留至小数点后第 4 位数字。取 nmax=1000 时，计算结果如表 1.
由表 1 可知对于同一步长 h，复合辛普森法求积分精度明显比复合梯形法求积的精度要高，且
当步长取不同值时即 h 越小时，积分精度越高。随着 n 的增大，复合辛普森法求积复合梯形积
分值复合梯形积分值的误差都逐渐减小，由于积分值的精度保留至小数点后第 4 位数字，故存
1
在一个最小的 h, 使得精度不能再被改善。对于复合梯形求积，当 n=694, 即 h=1.4409E-03 时，
计算精度不再改善。对于复合辛普森求积，当 n=202, 即 h=4.9505E-03 时，计算精度不再改善。
显然，当增加积分值的精度时，计算精度不再改善对于的步长 h 值将减小，但仍然存在。
表 1 复合梯形积分值复合梯形积分值和误差(部分数据)
辛普森积分值
(T(n))
-0.3268
-0.3958
-0.4157
-0.4248
-0.4298
-0.4329
-0.4351
-0.4366
-0.4377
-0.4386
⋮
-0.4443
-0.4443
-0.4443
-0.4443
-0.4443
-0.4443
-0.4443
-0.4444
-0.4444
-0.4444
-0.4444
⋮
辛普森积分值误差
1.1769E-01
4.8661E-02
2.8729E-02
1.9692E-02
1.4664E-02
1.1511E-02
9.3737E-03
7.8417E-03
6.6970E-03
5.8136E-03
⋮
1.4444E-04
1.4444E-04
1.4444E-04
1.4444E-04
1.4444E-04
1.4444E-04
1.4444E-04
4.4444E-05
4.4444E-05
4.4444E-05
4.4444E-05
⋮
等分数(n)
1
2
3
4
5
6
7
8
9
10
⋮
195
196
197
198
199
200
201
202
203
204
205
⋮
690
691
692
693
694
695
696
697
698
699
步长(h)
1.0000E+00
5.0000E-01
3.3333E-01
2.5000E-01
2.0000E-01
1.6667E-01
1.4286E-01
1.2500E-01
1.1111E-01
1.0000E-01
⋮
5.1282E-03
5.1020E-03
5.0761E-03
5.0505E-03
5.0251E-03
5.0000E-03
4.9751E-03
4.9505E-03
4.9261E-03
4.9020E-03
4.8780E-03
⋮
1.4493E-03
1.4472E-03
1.4451E-03
1.4430E-03
1.4409E-03
1.4388E-03
1.4368E-03
1.4347E-03
1.4327E-03
1.4306E-03
梯形积分值(T(n))
梯形积分值误差
0.0000
-0.2451
-0.3218
-0.3581
-0.3789
-0.3922
-0.4014
-0.4081
-0.4131
-0.4171
⋮
-0.4439
-0.4439
-0.4439
-0.4439
-0.4439
-0.4439
-0.4439
-0.4439
-0.4439
-0.4439
-0.4439
⋮
-0.4443
-0.4443
-0.4443
-0.4443
-0.4444
-0.4444
-0.4444
-0.4444
-0.4444
-0.4444
4.4444E-01
1.9938E-01
1.2266E-01
8.6340E-02
6.5535E-02
5.2212E-02
4.3028E-02
3.6354E-02
3.1311E-02
2.7382E-02
⋮
5.4444E-04
5.4444E-04
5.4444E-04
5.4444E-04
5.4444E-04
5.4444E-04
5.4444E-04
5.4444E-04
5.4444E-04
5.4444E-04
5.4444E-04
⋮
1.4444E-04
1.4444E-04
1.4444E-04
1.4444E-04
4.4444E-05
4.4444E-05
4.4444E-05
4.4444E-05
4.4444E-05
4.4444E-05
(2) 用龙贝格求积计算步骤如下, 取二分区间最大次数为 kmax.
①取 0,
k
h
, 求
b
a
(0)
T
0
h
2
(
f a
)
( )
f b
.
令 k=0(k 为区间[a,b]的二分次数)
②求梯形值 ( )
kT ，按照复合梯形公式(1)计算，并取
0
h
b a
2k
.
③求加速值，按(5)逐个求出 (
T
j
k j
)(
j
1,2,
, )
k
.
2
④若
(0)
T
k
(0)
T
1
k
(预定给定的精度), 则终止计算，并取 (0)
kT
I ，否则，令 k=k+1, 转而
继续从②计算.
( )
k
T
m
4
m
4
m
1
1)
(
k
T
1
m
1
4
m
1
( )
k
T
1
m
,
k
1,2,
(5)
计算时，若取 kmax=10, =10−4, 最终确定的二分次数 k=9, 输出的 T 表 2, 若=10−4 不变,
取 kmax=20, 输 出的 T 不 变, 显 然计 算时 正确 的 ， 最终 积分 近似 值为−0.444386. 计 算 的
MATLAB 程序见附录 3.
)
(kT
1
)
(kT
2
)
(kT
3
)
(kT
4
)
(kT
5
)
(kT
6
)
(kT
7
)
(kT
8
)
(kT
9
表 2 龙贝格求积 T 表
)
(
kT
0
0
-0.245065
-0.326753
-0.358104
-0.408090
-0.395784
-0.424752
-0.400386
-0.426683
-0.429475
-0.436603
-0.437393
-0.438389
-0.441361
-0.441678
-0.442031
-0.443244
-0.443370
-0.443494
-0.443981
-0.444030
-0.444074
-0.444301
-0.444267
-0.444377
-0.444286
-0.444384
-0.427101
-0.437563
-0.441746
-0.443397
-0.444041
-0.444290
-0.444386
-0.437604
-0.441763
-0.443403
-0.444043
-0.444291
-0.444386
-0.441767
-0.443405
-0.444044
-0.444291
-0.444386
-0.443405
-0.444044
-0.444044
-0.444291
-0.444386
-0.444291
-0.444386
-0.444291
-0.444386
-0.444386
k
0
1
2
3
4
5
6
7
8
9
(3) 用自适应辛普森积分, 使其精度达到 10−4.
设给定精度要求 0 , 计算积分的近似值。先取步长 h1=b−a，应用辛普森公式有
I
(
f
)
b
a
( )
f x dx
( , )
S a b
1
4
b a h
180
2
f
4
( ),
( , )
a b
其中
1( , )
S a b
把区间[a,b]对分，步长
h
2
h
1
2
b
2
h
6
a
( ) 4
f a
f a
，在每个小区间上用辛普森公式，得
( )
f b
h
1
2
I
(
f
)
( , )
S a b
2
b a
180
h
2
2
4
f
4
( ),
( , )
a b
(6)
(7)
(8)
其中，
( , )
S a b
2
S a
,
a b
2
h
2
6
h
2
6
,
S a
a b
2
a b
2
S
,
b
3
S
( ) 4
f a
,
b
,
a b
2
f a
h
2
2
(
f h
2
) ,
( + ) 4
f a h
2
f a
3
h
2
2
( ) .
f b
( )
I f
S
,在适当假设下, 可以推出必须有 1
( , )
2 15
S a b
S
，若满足此不等式，
要满足
2
此时 2( , )
S a b 作为 (
f 的近似，则可以达到给定的误差精度. 若不等式不成立, 则分别对子
)
I
a b
a b b
a b
a
b
区间 ,
,
,
2
2
a b
S a
2
, 得到 3
h
2
2
S
3
，
.
,
,
a b
2
进行辛普森公式, 此时步长
h
3
b
2
a b
2
I S
3
2
及
,
2
,
I S a
3
只要分别考察
是否成立. 对于满足要求的区间不再
细分，对于对不满足要求的还要继续上述过程，直到满足要求为止.计算程序如附录 4 所示.
最终，计算积分值为−0.444438, 实际误差为 6.3636−06.
4
附录 1 复合梯形求积 MATLAB 程序
function T=Tixing(nmax)
T=zeros(nmax,1);
for n=1:nmax
T(n)=Tn(n);
end
fid=fopen('tixing.txt','w');fprintf(fid,'%8.4f\n',T);fclose(fid);
function TT=Tn(n)
a=0;b=1;
T0=f(a)+f(b);
h=(b-a)/n;
A=0;
if n==1
TT=h/2*T0;
else
x=zeros(1,n-1);
for k=1:n-1
x(k)=a+k*h;
A=A+f(x(k));
end
TT=h/2*(T0+2*A);
end
end
function y=f(x)
if x==0
y=0;
else
end
y=sqrt(x)*log(x);
end
end
5
附录 2 复合辛普森求积 MATLAB 程序
function S=Simpson(nmax)
S=zeros(nmax,1);
for n=1:nmax
S(n)=Sn(n);
end
fid=fopen('Simpson.txt','w');fprintf(fid,'%8.4f\n',S);fclose(fid);
function [SS]=Sn(n)
a=0;b=1;
S0=f(a)+f(b);
h=(b-a)/n;
A=0;B=0;
if n==1
SS=h/6*(S0+4*f(a+1/2*h));
else
x=zeros(1,n-1);
x1=zeros(1,n);
for k=1:n-1
x(k)=a+k*h;
A=A+f(x(k));
end
for k=0:n-1
k1=k+1;
x1(k1)=a+(k+1/2)*h;
B=B+f(x1(k1));
end
SS=h/6*(S0+2*A+4*B);
end
function y=f(xx)
if xx==0
y=0;
else
end
y=sqrt(xx)*log(xx);
end
end
end
6
附录 3 龙贝格求积 MATLAB 程序
function [T1,T,k]=Romberg (kmax,esp)
T=zeros(kmax+1);
k=0;
n=2^k;
T(k+1,1)=T_0(n);
delta=1;
while ((delta>=esp)&(k<=kmax))
k=k+1;
n=2^k;
T(k+1,1)=T_0(n);
for m=1:k
T(k+1,m+1)=4^m/(4^m-1)*T(k+1,m)-1/(4^m-1)*T(k,m);
end
delta=abs(T(k+1,k+1)-T(k,k));
end
T1=T(1:1+k,1:1+k); %最终输出的 T 表
function [Tn]=T_0(n)
a=0;b=1;
T0=f(a)+f(b);
h=(b-a)/n;
A=0;
if n==1
Tn=h/2*T0;
else
x=zeros(1,n-1);
for kk=1:n-1
x(kk)=a+kk*h;
A=A+f(x(kk));
end
Tn=h/2*(T0+2*A);
end
function y=f(x)
if x==0
y=0;
else
end
y=sqrt(x)*log(x);
end
end
end
7
附录 4 自适应辛普森求积 MATLAB 程序
function [S,eps1]=Adaptive_Simpson(eps)
a=0;b=1;
h1=b-a;
S=h1/6*(f(a)+4*f(a+h1/2)+f(b)); %第一次 Simpson 公式积分值
S=Simpson(a,b,eps,S);
eps1=abs(-4/9-S);
function SS=Simpson(aa,bb,eps,S)
%最终求得的自适应方法积分值
%实际误差
h=(bb-aa)/2;
S1=h/6*(f(aa)+4*f(aa+h/2)+f(aa+h));
S2=h/6*(f(aa+h)+4*f(aa+3/2*h)+f(bb));
if abs(S-S1-S2)
相关推荐
2023年江西萍乡中考道德与法治真题及答案.doc
2012年重庆南川中考生物真题及答案.doc
2013年江西师范大学地理学综合及文艺理论基础考研真题.doc
2020年四川甘孜小升初语文真题及答案I卷.doc
2020年注册岩土工程师专业基础考试真题及答案.doc
2023-2024学年福建省厦门市九年级上学期数学月考试题及答案.doc
2021-2022学年辽宁省沈阳市大东区九年级上学期语文期末试题及答案.doc
2022-2023学年北京东城区初三第一学期物理期末试卷及答案.doc
2018上半年江西教师资格初中地理学科知识与教学能力真题及答案.doc
2012年河北国家公务员申论考试真题及答案-省级.doc
2020-2021学年江苏省扬州市江都区邵樊片九年级上学期数学第一次质量检测试题及答案.doc
2022下半年黑龙江教师资格证中学综合素质真题及答案.doc
