第七章 数值积分
( Numerical Integration)
武汉大学数学与统计学院内容提纲
数值积分的必要性
求积公式及其代数精度
插值型求积公式
Newton-Cotes公式及数值稳定性
复化求积公式及误差估计数值积分的必要性本章主要讨论如下形式的一元函数积分在微积分里,按 Newton-Leibniz公式 求定积分要求被积函数 f(x)
有解析表达式 ;
f(x)的原函数 F(x)为初等函数,
( ) ( ) ( ) ( )baI f f x dx F b F a

b
a
dxxffI )()(
实际问题
1.2 f(x)的原函数 F(x)不能用初等函数表示例如函数,
2,1,
ln
1,s i n,c o s,s i n 322 xex
xx
xxx
考虑一个实际问题,
建筑上用的一种铝制波纹瓦是用一种机器将一块平整的铝板压制而成的,
假若要求波纹瓦长 4英尺,每个波纹的高度 (从中心线 )为 1英寸,且每个波纹以近似 2π 英寸为一个周期,求制做一块波纹瓦所需铝板的长度 L.
这个问题就是要求由函数 f(x)=sin x给定的曲线,从 x=0到 x=48英寸间的 弧长 L.
由微积分学我们知道,所求的弧长可表示为,
dxxdxxfL 480 2480 2' )( co s1))((1
上述积分称为第二类椭圆积分,它不能用普通方法来计算,
2,有些被积函数其原函数虽然可以用初等函数表示成有限形式,但表达式相当复杂,计算极不方便,例如函数
32 22?xx
并不复杂,但它的原函数却 十分复杂,
)322l n (216 93216 33241 2222 xxxxxx
3.f(x)没有解析表达式,只有数表形式,
x 1 2 3 4 5
f(x) 4 4.5 6 8 8.5
这些都说明,通过原函数来计算积分有它的局限性,因而,研究关于积分的数值方法具有很重要的实际意义,
求积公式及其代数精度求积公式的概念积分值在几何上可解释为由 x=a,x=b,y=0和 y=f(x)
所围成的 曲边梯形的面积,积分计算之所以有困难,就是因为这个曲边梯形有一条边 y=f(x)
是曲的,

b
a
dxxffI )()(
依据 积分中值定理,对于连续函数 f(x),在
[a,b]内存在一点 ξ,使得
)()()()(?fabdxxffI
b
a

称 f(ξ)为区间 [a,b]的平均高度,问题在于点 ξ的具体位置一般是不知道的,这样,只要对平均高度 f(ξ)提供一种 算法,相应地便获得一种 数值求积方法,
如果简单地 选取区间 [a,b]的一个端点或区间中点的高度作为平均高度,这样建立的求积公式分别是,
左矩形公式,I(f)≈(b-a)f(a)
右矩形公式,I(f)≈(b-a)f(b)
中矩形公式,I(f)≈(b-a)f[(a+b)/2]
此外,众所周知的 梯形公式,
I(f)≈(b-a)[f(a)+f(b)]/2
和 Simpson公式,
I(f)≈(b-a)[f(a)+4f((a+b)/2)+f(b)]/6
则分别可以看作用 a,b,c=(a+b)/2,三点高度的加权平均值 [f(a)+f(b)]/2

[f(a)+4f(c)+f(b)]/6
作为 平均高度 f(ξ)的近似值,
更一般地,取区间 [a,b]内 n+1个点 {xi},(i=0,1,
2,…n) 处的高度 {f(xi)} (i=0,1,…,n) 通过 加权平均 的方法近似地得出平均高度 f(ξ),这类求积方法称为 机械求积,
)()()(
0
i
b
a
n
i
i xfabdxxf

或写成,
数值积分公式 求积系数求积节点
)()(
0
k
b
a
n
k
k xfAdxxf
(1)

)2()()(
0
k
n
k
kn xfAfI?
)3(,)()()()()(
0

b
a
n
k
kkn xfAdxxffIfIfR
称 (2)为数值求积公式,(3)为求积公式余项 (误差 ).
构造或确定一个求积公式,要讨论解决的问题有
(i) 确定求积系数 Ak和求积节点 xk ;
(ii) 求积公式的误差估计和收敛性为了构造形如式 (2)的求积公式,需要提供一种 判定求积方法精度高低准则求积公式的代数精度定义 1 称求积公式 (2)具有 m次代数精度,如果它满足如下两个条件,
(i)对所有次数 ≤ m次的多项式,有
(ii)存在 m+1次多项式,使得
)(xPm
0)()()( mnmm PIPIPR
)(1 xPm?
0)()()( 111 mnmm PIPIPR
定义 1中的条件 (i),(ii)等价于,
0)()(
)0(,0)()()()(
1?

m
k
n
kk
xRii
mkxIxIxRi
插值型求积公式在积分区间 [a,b] 上 取 n+1个节点 xi,i=0,1,
2,…,n,作 f(x)的 n次代数插值多项式 (拉格朗日插值公式),
则有为插值余项于是有
n
j
jjn xfxlxL
0
)()()(
)()()( xRxLxf nn
)()!1( )()( 1
)1(
xwnfxR n
n

取称 (4)式为插值型求积公式,其中 求积系 Ak由 (5)
式确定,






b
a
j
n
j
b
a
j
b
a
n
b
a
n
b
a
dxxRxfdxxl
dxxRdxxLdxxf
)()()(
)()()(
0
(4)
(5)

b
a
b
a k
n
k
k dxxlxfdxxf )()()(
0
Ak
0
()
()
nb
i
k a
i ki
ik
xxA d x
xx?


由 节点 决定,
与 f(x) 无关。
推论 1 求积系数满足,
0
( 1 )
0
[ ] ( ) ( ) [ ( ) ( ) ]
()
()
( 1 ) !
nbb
k k naa
k
n n
b
x
ka
k
R f f x d x A f x f x L x d x
f
x x d x
n




误 差定理 1形如 的求积公式至少有 n 次代数精度
该 公式为 插值型 (即,)
n
k
kk xfA
0
)(
ba kk dxxlA )(
abA
n
j
j
0
取 节点为 等距分布,
,,0,1,.,,,i bax a i h h i nn
dxxx xxA nx
x ij ji
j
i


0 )(
)(

n
ji
inn
ji
dtjtinin abdthhji hjt
00
)()!(! )1)(()( )(
令 htax
Cotes系数 ()n
kC
注,Cotes 系数仅取决于 n 和
k,可查表得到。与 f (x)
及区间 [a,b]均无关。
Newton-Cotes 公式由此构造的插值型求积公式称为 Newton-Cotes公式,此时 求积系数记
dtkt
njnj
C
n n
jkk
jn
n
j


0,0
)( )(
)!(!
)1(

njCabA njj,,2,1,0,)( )(
)()()(
0
)(
j
b
a
n
j
n
j xfCabdxxf

求积公式 (4)变为
(6)
(7)
(8)
称 (8)式为 n阶 闭型 Newton-Cotes求积公式,
),(,)()(
)!1(
)(
)!1(
)(
)(
0
0
)1(
2
1
)1(
badtjtf
n
h
dxxw
n
f
fR
n
n
j
n
n
n
b
a
n



Newton-Cotes公式的误差为,
与 x有关注意,由 (6)式确定的 Cotes系数只与 j和 n有关,
与 f(x)和积分区间 [a,b]无关,且 满足,
1
0
)(
n
j
n
jC
(9)
定理 2 当阶数 n为偶数时,Newton-Cotes
公式 (8)至少具有 n+1次代数精度,
证明 只需验证当 n为偶数时,Newton-Cotes
公式对 f(x)=xn+1的余项为零,
由于 f(x)=xn+1,所以 f(n+1)(x)=(n+1)!,由式 (9)得

n
n
j
n dtjthfR
0 0
2 )()(
引进变换 t=u+n/2,因为 n为偶数,故 n/2为整数,
于是有

2
2 0
2 )
2
()(
n
n
n
j
n dujnuhfR
据此可断定 R(f)=0,因为上述被积函数是个奇函数,
Newton-Cotes公式的数值稳定性现在讨论 舍入误差 对计算结果产生的影响,设用公式近似计算积分时,其中计算函数值 f(xj)有误差 εj(j=0,1,2,…,n),设计算 Cj(n)没有误差,中间计算过程中的舍入误差也不考虑,则在式 (10 )的计算中,由 εj引起的误差为

n
j
j
n
jn xfCabfI
0
)( )()()(

b
a
dxxffI )()(
(10)




n
j
j
n
j
n
j
jj
n
j
n
j
j
n
jn
Cab
xfCabxfCabe
0
)(
0
)(
0
)(
)(
))(()()()(
如果 Cj(n)都是正数,并设 ||m a x
0 jnj


)(||)(||
0
)( abCabe
n
j
n
jn

故 en是有界的,即由 εj引起的误差受到控制,不超过
ε的 (b-a)倍,保证了 数值计算的稳定性,而当 n>7
时,Cj(n)将出现 负数,?
n
j
n
jC
0
)( ||
保证数值稳定性,因此高阶公式不宜采用,有实用价值的仅仅是几种 低阶的求积公式,
将随 n增大,因而 不能则有
2
1,
2
1 )1(
1)1(0 CC
n = 1:
)]()([2)( bfafabdxxfba
Trapezoidal Rule
dxbxaxffR ba x ))((!2 )(][
/* 令 x = a+th,h = b?a,用中值定理 */
1,],[,)(12
1 3 abhbafh
代数精度 = 1
n = 2:
6
1,
3
2,
6
1 )2(
2
)2(
1
)2(
0 CCC
)]()(4)([6)( 2 bffafabdxxf baba
Simpson’s Rule
代数精度 = 3
2,),(,)(90
1][ )4(5 abhbafhfR
n = 4,Cotes Rule,代数精度 = 5,)(
9458][ )6(7?fhfR
复化型求积公式高次插值有 Runge 现象,高阶 Newton-
Cotes公式会出现 数值不稳定,低阶
Newton-Cotes公式有时又不能满足精度要求,解决这个矛盾的办法是将积分区间
[a,b]分成若干小区间,在每个小区间上用低阶求积公式计算,然后将它们加起来,
这就是 复化求积方法,
复化梯形公式:
,( 0,.,,,)ibah x a i h i nn
在每个 上用梯形公式:],[
1?ii xx
1 1
1( ) [ ( ) ( ) ],0,.,,,12
i
i
x ii
iix
xxf x d x f x f x i n

1
1
( ) 2 ( ) ( )2
n
k
i
h f a f x f b?


1
1
0
( ) [ ( ) ( ) ]2
nb
iia
i
hf x dx f x f x?

= Tn
1
321
0
0
2
()
[ ] [ ( ) ] ( )
1 2 1 2
( ) ( ),(,)
12
n
in
i
i
i
f
hh
R f f b a
n
h
b a f a b




/*中值定理 */
复化梯形公式积分法收敛性由上述的误差估计式可知,当 f(x)?C2[a,b] 时,只要
h?0时 数列 Tn(f)? I(f),且收敛速度为二阶 O(h2).
但是 f(x)?C2[a,b] 条件相对苛刻,现假定 f(x)在 [a,b]
上 Riemann可积,讨论复化求积公式的收敛性
1
1
0
( ) [ ( ) ( ) ]2
n
n i i
i
hT f f x f x?

1
01
1 ( ( ) ( ) )
2
nn
i i i i
ii
f x x f x x


0
1l i m ( ) l i m ( ) [ ( ) ( ) ] ( )
2nn nT f T f I f I f I f
,( 0,1,.,,,2 )2 jbah x a j h j nn
2 2 2 1 2
1
( ) [ ( ) 4 ( ) ( ) ]3
n
j j j
j
hI f f x f x f x


22jx? 21jx? 2jx
2 2 2 1 2
1
( ) [ ( ) 4 ( ) ( ) ]
3
n
n j j j
j
h
S f f x f x f x

复化 Simpson 公式:
1
2 1 2
11
( ) [ ( ) 4 ( ) 2 ( ) ( ) ]3
nnb
jja
jj
hf x dx f a f x f x f b?


复化 Simpson公式积分法误差估计
5
( 4 )
2 2 2( 1,) ( ),,9 0 2j j j j j
h b aR f f h x x
n

每个子区间上的误差估计式为
( 4 ) 4( ) ( ) ( ),
180n
baI f S f f h a b
将 n个子区间的误差相加得
5
( 4 )
1
( ) ( ) ( )90
n
nj
j
hI f S f f?

由闭区间上连续函数的介值性质可知在 [a,b]上至少存在一点?,使
( 4 ) ( 4 )
1
1( ) ( )n
j
j
ff n

可见,当 f(x)有四阶导数时,复化 Simpson公式具有 4阶收敛,
例 对于函数,
x
xxf s in)(? 试用数据表计算积分
10 s in)( dxx xfI
x f (x)
0 1
1/8 0.9973978
2/8 0.9896158
3/8 0.9767267
4/8 0.9588510
5/8 0.9361556
6/8 0.9088516
7/8 0.8771925
1 0.8414709
解 将区间 [0,1]划分为 8等分,应用复化梯形法 求得
9 4 5 6 9 0 9.0
)1(
8
2)0(
28
1
)()(2)(
2
7
1
7
1
8


k
k
k
f
k
ff
bfxfaf
h
T
应用 复化 Simpson法 计算,得
9 4 6 0 8 3 2.0
)1(
8
2
2
8
12
4)0(
38
1
)()(2)(4)(
3
4
1
3
1
4
1
3
1
2124








k k
j j
jj
f
k
f
k
ff
bfxfxfaf
h
S
比较上面两个结果 T8和 S4,它们都需要提供 9个点上的函数值工作量基本相同,然而精度却差别很大,
同积分的 准确值 I(f)=0.9460831比较,复化梯形法的结果 T8=0.9456909只有 两位有效数字,而复化
Simpson法的结果 S4=0.9460832却有 六位有效数字,
习 题
p388,1,5;
p389,8.