计算方法讲义：七 数值积分

第七章 数值积分如果函数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 的两点求积公式。

本题的解法很多，结果也不一定相同。