第七章 数值积分如果函数f(x)在区间[a,b]上连续,且原函数为F(x),则可用牛顿―莱布尼兹公式:)()()(a F b F dx x f b a-=⎰来求得定积分。
然而很多函数无法用牛顿―莱布尼兹公式求积分。
一个简单被积函数,例如,其不定积分可能很复杂,见下面的MA TLAB 实例: >> syms a b c x>> int(sqrt(a+b*x+c*x*x),x)ans=1/4*(2*c*x+b)/c*(a+b*x+c*x^2)^(1/2)+1/2/c^(1/2)*log((1/2*b+c*x )/c^(1/2)+(a+b*x+c*x^2)^(1/2))*a-1/8/c^(3/2)*log((1/2*b+c*x)/c^(1/2)+(a+b*x+c*x^2)^(1/2))*b^2所以有必要研究简单、高效的计算定积分的方法(即数值积分方法)。
数值积分的基本思想是构造一个简单函数P n (x )来近似代替被积分函数f (x ),然后通过求⎰ba n dx x P )(得⎰ba dx x f )(的近似值。
7.1 插值型求积公式设⎰=ba dx x f I )(*,插值型求积公式就是构造插值多项式P n (x ),使⎰=≈ba n dx x P I I )(*。
构造以a ,b 为结点的线性插值多项式)()()(1b f ab ax a f b a b x x P --+--=,[])()()(21)()()(1b f a f a b dx b f a b a x a f b a b x dx x P T ba ba +-=⎥⎦⎤⎢⎣⎡--+--==⎰⎰称为梯形公式。
以a , 2ba c +=,b 为三个插值节点,构造二次插值多项式)())(())(( )())(())(()())(())(()(2b f c b a b c x a x c f b c a c b x a x a f b a c a b x c x x P ----+----+----=,则可以推出)()()()(2102b f c f a f dx x P S baλλλ++===⎰,)(61))(())((0a b dx b a c a b x c x ba-=----=⎰λ,)(64))(())((1a b dx b c a c b x a x ba-=----=⎰λ,)(61))(())((2a b dx c b a b c x a x b a -=----=⎰λ。
由此得公式:[])()(4)(6b fc f a f ab S ++-=,称为辛卜生(Sinpson )求积公式。
根据经典拉格朗日插值公式)()()(0k nk k n x f x l x P ∑==,代入求定积分则有)()()()(0k nk b a k k nk ba k x f dx x l dx x f x l I ⋅==∑⎰∑⎰==,引入记号dx x l ba k k )(⎰=λ,)(0k nk k x f I ∑==λ,λk 为求积系数,x k 为求积节点。
注意:一积分结果为函数值的一个代数和,二是ab dx x l nk ba k -=∑⎰=0)(。
如果积分区间比较大,直接使用上述求积公式精度难以保证。
可对f (x )用分段抛物插值。
通常采取的办法是复化求积方法: (1)等分求积区间,比如取步长nab h -=,分[a, b]为n 等分,分点为kh x x k +=0,k = 0, 1, 2,…, n 。
(2)在区间 [x k , x k+1]上使用以上求积公式求得I k 。
(3)取和值∑-==10n k k I I 作为整个区间上的积分值。
将梯型公式和辛卜生公式应用于各子区间[]1,+k k x x ()1,,0-=n k 上得到子区间的定积分,再将子区间的定积分加起来得到整个区间的定积分近似值,相关公式称为复化梯型公式和复化辛卜生公式。
相对于复化梯型公式,复化辛卜生公式是一种精度较高的求积公式。
例如对于复化梯型公式,令[])()(21++=k k k x f x f hI ,则∑-==10n k k n I T [])()(2110+-=+=∑k k n k x f x f h ⎥⎦⎤⎢⎣⎡++=∑-=)()(2)(211b f x f a f h n k kx k 0 1/8 1/4 3/8 1/2 5/8 3/4 7/8 1 f (x k ) 4 3.93846 3.7647 3.50685 3.2 2.8764 2.46 2.26549 2 例 利用数据表计算积分 dx xI ⎰+=12*14( 1415926.3|arctg 410*===πx I )。
解:取n = 8用复化梯形公式:()13899.31872432852212832412812)0(21818=⎥⎦⎤+⎪⎭⎫⎝⎛+⎪⎭⎫ ⎝⎛+⎪⎭⎫ ⎝⎛+⎢⎣⎡⎪⎭⎫ ⎝⎛+⎪⎭⎫ ⎝⎛+⎪⎭⎫ ⎝⎛+⎪⎭⎫ ⎝⎛+⨯=f f f f f f f f f T 7.2 变步长梯形方法使用复化求积公式须给出合适的步长,步长太大精度难保证,步长太小会增加计算量,事先给出一个合适的步长是十分困难的。
递推公式避免了老节点的重复计算,使计算量减少了一半。
变步长积分法思想是将区间逐次对分,比较前后两次计算结果,若满足精度要求就停止,否则再次对分,直到到达精度要求为止。
设将区间[a, b] n 等分,共有n+1个分点,按复化梯形公式计算T n ,需要计算n+1个f (x )的值。
T 2n 的全部分点中有n+1个是原有的点。
小区间[x k , x k +1]经过二分增加分点21+k x 后,用复化梯形公式得积分为:[])()(2)(4121++++=k k k k x f x f x f hI ,因此有[][]∑∑∑∑-=+-=+-=+++-=+=++=++=110101112)(221)(2)()(4)()(2)(4212121n k k n n k k n k k k k k k n k n x f h T x f h x f x f h x f x f x f hT例 计算。
解:根据梯形公式和复化梯形公式2/)]1()0([1f f T +=,∑-=++=15.02)(221n k k n n x f h T T ,于是有n 1 2 4 8 16 32 T n 0.9397 0.9445 0.9456 0.9459 0.9461 0.94617.3 求积公式的误差P n (x )是f (x )的n 次插值多项式,当)(x f 本身就是次数不超过n 的多项式时)()(x P x f n ≡,求积公式)()()(*0*k nk k ba n ba x f dx x P dx x f I ∑⎰⎰====λ是精确的。
由于a b nk k-=∑=0λ,若f (x k )的舍入误差小于ε ,则)()()()()(*0**a b x f x f x f x f I I k k nk k k nk k k nk k-<-≤-=-∑∑∑===ελλλ。
所以舍入误差对数值积分的影响不大。
应用插值多项式余项定理)()!1()()(1)1(x n f x R n n +++=ωξ,对于插值多项式次数为1的情况有:))((2)()(b x a x f x R --''=ξ,可以证明梯形公式的截断误差为:3*)(12)())((2))((a b f dx b x a x x f T I ba-''=--''=-⎰ξξ-(注意:这里用到了积分中值定理:设)(x f 在区间[a,b]上连续,)(x g 在[a,b]区间上可积且不变号,则在[a,b]区间上至少有一个ξ满足⎰⎰=babadx x g f dx x g x f )()()()(ξ)。
将[a, b]区间n 等分,取nab h -=考虑复化梯形积分公式的误差,这个误差是n 个等分区间段上得误差之和,即12)()(-12)(-)(12-12)(-2310313*h f a b h f n f h h f T I n i i n i i n ξξξξ''-=''=''=''=-∑∑-=-=(注意:这里用到了介值定理,即对于连续函数)(x f 和自然数n ,存在ξ使。
7.4 收敛条件及收敛加速梯形法简单,但精度低,收敛的速度慢。
如何提高收敛速度?设I 是精确积分值,根据复化梯形公式的余项表达式可知:),(12)()(-2*b a h f a b T I n ∈''-=-ηη,,),(12)2/)(()(-*2*2*b a h f a b T I n ∈''-=-ηη,。
假定)()(*''≈''ηηf f , 则有41*2*=--n n T I T I 。
整理得:)(3122*n n n T T T I -=-。
可见只要二分前后T n 与T 2n 相当接近,就可以保证T 2n 的误差很小。
T 2n 的误差大致等于)(312n n T T -,用误差值作为T 2n 的补偿,可期望所得到的)(3122*n n n T T T I -+=,可能是更好的结果。
也可以这样考虑,将所有T n 看做构成一个函数T ,变量是h 2。
当h 趋近0时,T(h 2)接近I *,即T(0)=I *,连接点),2n n T h (和),121++n n T h (得一直线,其方程为:1221221221++++--+--=n nn n n n n n T h h h x T h h h x y ,延伸该直线与Y 轴相交(见下图),3-4T ~11221221221n n n n n n n n n n T T T h h h T h h h +++++=--+-=。
这就是一种迭代的加速。
图 一种迭代加速7.5 高斯型求积公式在插值型求积公式中,插值节点是事先固定的,有时还进一步限定是等距的。
是否可以在[a, b]上自由选择节点的位置,使精度提高?称:)()(0i ni k ba x f A dx x f ∑⎰=≈为一般求积公式,这里A k 为不依赖f (x )的常数。
若对任意不高于m 次的多项式精确成立,而对于x m +1不能精确成立,就说具有m 次代数精确度。
下面讨论最高代数精确度的求积公式(叫做高斯型求积公式)。
例 求形如 ⎰-+≈111100)()()(x f A x f A dx x f 的两点求积公式。
本题的解法很多,结果也不一定相同。