1.变步长Romberg 积分法的原理
复化求积方法对于提高精度是行之有效的方法,但复化公式的一个主要缺点在于要事先估计出部长。若步长过大,则精度难于保证;若步长过小,则计算量又不会太大。而用复化公式的截断误差来估计步长,其结果是步长往往过小,而且''()f x 和(4)()f x 在区间[,]a b 上的上界M 的估计是较为困难的。在实际计算中通常采用变步长的方法,即把步长逐次分半(也就是把步长二等分),直到达到某种精度为止,这种方法就是Romberg 积分法的思想。
在步长的逐步分半过程中,要解决两个问题: 1. 在计算出N T 后,如何计算2N T ,即导出2N T 和N T 之间的递推公式; 2. 在计算出N T 后,如何估计其误差,即算法的终止的准则是什么。
首先推导梯形值的递推公式,在计算N T 时,需要计算1N +个点处的函数值在计算出N T 后,在计算2N T 时,需将每个子区间再做二等分,共新增N 个节点。为了避免重复计算,计算2N T 时,将已计算的1N +个点的数值保留下来,只计算新增N 个节点处的值。为此,把2N T 表示成两部分之和,即
21
2221
122211
122211
1
[()()2()]21
[()()2()2((21))]21[()()2()]((21))]22N N
N N k N N N N N k k N N
N N N N k k T h f a f b f a kh h f a f b f a kh f a k h h f a f b f a kh h f a k h -=-==-===+++=+++++-=+++++-∑∑∑∑∑ 由此得到梯形值递推公式
2221
1
((21))12N
N
N N N k T T h f a k h ==++-∑
因此
1
11212122,[()()],211
,()
22h h b a T f a f b h h T T h f a h =-=
+==++
由复化梯形公式的截断误差有
2''
112''2222(),12(),12
N N N
N b a I T h f a b b a I T h f a b
ηηηη--=-≤≤--=-≤≤ 若''()f x 变化不大时,即''''12()()f f ηη≈,则有
22241
()
2413
N N N N N T T I T T T -≈
=+--
式(2)表明,用2N T 作为定积分I 的近似值,其误差大致为21
()3
N N T T -,
因此其终止条件为
2N N T T ε-≤
其中ε是预先给定的精度。
积分公式
将上述方法不断推广下去,可以得到一个求积分的序列,而且这个序列很快收敛到所求的定积分。记
(0)N N T T =,将区间N 等分的梯形值。(1)N N T S =,将区间N 等分的Simpson (2)
N N T C =,将区间N 等分的Cotes 。(3)N
N T R =,将区间N 等分的Romberg 。 由其可构造一个序列(){}k N T ,次序列称为Romberg 序列,并满足如下递推关系:
(0)
(0)(0)1
211[()()],((21)),2222N N N k b a b a b a T
f a f b T T f a k N N
=---=+=++-∑ (1)(1)
24,1,2,41k k k k N N
N
k
T T T k ---==
-
以上递推公式就是Romberg 积分递推公式。
积分程序
1. 置1N =,精度要求ε,1h b a =-;
2. 计算(0)1[()()]2
b a
T f a f b -=+; 3. 置22
N N
h h =,并计算(0)
(0)211((21))222N N N k b a b a T T f a k N N =--=++-∑;
4. 置,2,1;M N N N K ===
5. 计算(1)(1)
2441
k k k k M M
M
k T T T ---=-;
6. 若 1M =,则转(7);否则置2
M
M =
,1k k =+转(5); 7. 若()(1)11k k T T ε--≤,则停止计算(输出()1k T ),否则转(3)。
积分法的应用
function [T,n] = romb(f,a,b,eps) double R ;
if nargin<4,eps=1e-8; end
h=b-a;R(1,1)=(h/2)*(feval(f,a)+feval(f,b)); n=1;J=0;err=1; while (err>eps)
J=J+1;h=h/2;S=0; for i=1:n
x=a+h*(2*i-1); S=S+feval(f,x); end
R(J+1,1)=R(J,1)/2+h*S; for k=1:J
R(J+1,k+1)=(4^k*R(J+1,k)-R(J,k))/(4^k-1); end
err=abs(R(J+1,J+1)-R(J+1,J)); n=2*n; end R;
T=R(J+1,J+1) End
其中输入项:f 为被积函数,ab 为积分区间的端点值,ep 为积分精度;输出项:T 是逐次积分表值,n 是迭代次数,R 是最后积分值。 程序调用
可以将被积分函数编成函数文件,也可以直接使用内联函数来表示被积分函数,示例如下:
>>f=inline('1/(1+x.^2)','x'); >> [T,n,R]=romb(f,2,9,1e-9)
运行后得出其迭代次数,最终积分结果以及龙贝格积分矩阵如表2-1所示, 迭代次数N=64,最终的积分值R=.
表2-1 龙贝格积分矩阵
3.课本例题求解
2
1
1
1
2
0000
1ln(1)ln(1)sin()(1),(2),(3),(4)1+1x x x dx dx dx dx x x x x π
+++????
1 当迭代精度ep=1e-9的条件下,迭代次数N=32,迭代结果R=
表2-2 式1对应的龙贝格积分矩阵
2 当迭代精度ep=1e-9的条件下,迭代次数N=32,迭代结果R=.
表2-3 式2对应的龙贝格积分矩阵
3对于积分1
ln(1)
x dx x +?
,由于积分下限0为其奇点,理论上无法进行数值积分,本题中近似取下限为1*10-9来进行计算。当迭代精度ep=1e-9的条件下,迭代次
数N=16,迭代结果R=.
表2-4 式3对应的龙贝格积分矩阵
4.对于积分2
0sin()x
dx x
π
?,同样积分下限0为积分函数的奇点,理论上无法进行数
值积分运算,本题中仍取积分下限近似为1*10-9进行计算。当迭代精度ep=1e-9的条件下,迭代次数N=16,迭代结果R=.
表2-5 式4对应的龙贝格积分矩阵