数值计算方法第七章

) Rn(f ) = 0
能源学院，航空航天中心
18
Newton-Cotes求积─误差估计
定理7.4
设函数 ( ) 2 C2[ ] ，则梯形积分公式 f x a, b Zb
a
f
( )d xx
⇡
b
a 2
( )+ ( ) fa fb
的截断误差是 R1(f ) =
(b
a)3 12
f
00(⌘),
⌘
2 [a, b]
)
+
f
(xk+1)
k=0
k=0
余项：
R=
(b a)5 2880n4
f
(4)(⌘)
能源学院，航空航天中心
23
变步长梯形积分法
能源学院，航空航天中心
24
变步长梯形积分公式
前面的方法都是固定步长的积分方法。实际应用时往往给定 一一个精度要求，要获得指定精度要求的积分，一一种方法是自 动选择积分步长。即在积分过程中，将步长逐步折半，反复 应用复合积分公式，直到满足精度要求。
Akf (xk) 的系数 Ak > 0
a
k=0
则求积公式稳定。
能源学院，航空航天中心
13
Newton-Cotes 公式
n 等分节点的插值型求积公式
能源学院，航空航天中心
14
Newton-Cotes 求积公式
将区间 [a, b]
n 等分，步长为 h = b a ，节点为 n
xk = a + kh
(k = 0, 1, · · · , n) ，构造插值型求积公式
A1
+
8A2
+
27A3
=
81 4
.
A0 =
3 8,
A1
=
9 8,
A2
=
9 8,
A3 =
3 8
Z3
0
f
( )d xx
⇡
3 8
(0) + 3 (1)
f
f
+ 3 (2) + (3) ff
能源学院，航空航天中心
8
插值型求积公式
根据节点处的函数值，构造插值多项式
pn
() x
wk.baidu.com
，用插值多项式近似
被积函数
() fx
能源学院，航空航天中心
3
基本思想
Zb
= ( )d
I
fx x
a
定积分的几何意义是求函数
( ) 对应的曲线在区间 fx
[a,
b]
所包络
的曲边梯形的面积。数值方法就是用别的曲线或者折线代替原曲
线，计算其对应的曲边形的面积。
最简单的数值方法就是用矩形来近似曲边梯形，计算公式如下
Zb
= ( )d ⇡ ( ) ( + (1 ) ) 2 [0 1]
22
复合 Simpson 求积公式
将子区间 [xk, xk+1]
对分，中点记为
xk+0.5 =
xk + xk+1 2
。
由 Simpson 求积公式，
⇣
⌘
Ik
=
h 6
f
(xk
)
+
4 f
(xk+0.5
)
+
f
(xk+1)
对子区间的积分求和：
nX1
nX1 ⇣
⌘
I⇡
Ik
=
h 6
f
(xk
)
+
4 f
(xk+0.5
I
f x x b a f ✓a
✓b, ✓ ,
a
✓ = 0 时，右矩形公式， ✓ = 0.5 时，中矩形公式，✓ = 1 时，左 矩形公式。
能源学院，航空航天中心
4
基本思想
若用梯形来近似曲边梯形，得到梯形积分公式：
Zb
= I
a
( )d ⇡ b fx x
a 2
( )+ ( ) fa fb
用过三点
A(a, f (a)), B(b, f (b)), C( a
点互异，则有唯一一解。
Zb
Xn
= I
( )d ⇡ fx x
Akf (xk)
a
k=0
• 待定系数法用解线性方程组的方法来构造求积公式。 • 精度给定，节点给定之后，求积公式就确定了。
能源学院，航空航天中心
7
代数精度─待定系数法例题
例7.1
用待定系数法构造具有三次代数精度的求积公式：
Z3
f
( )d xx
Ck(n)
=
1 n
Zn
0
0 Yn
@
i=0,i6=k
t k
1 i A dt = i
( 1)n nk!(n
k Z n Yn
k)!
0
(t
i=0,i6=k
i)dt
能源学院，航空航天中心
15
低阶 Newton-Cotes 公式
当 n=1 时，Cotes 系数为
C0(1) =
Z1t 00
1 1
dt
=
1 2
C1(1)
对 n 个子区间积分求和：
余项：
nX1
nX1
⇡ I
Ik =
h 2
f (xk) + f (xk+1)
k=0
k=0
⇣
nX1
⌘
=
h 2
( )+ ( )+2 fa fb
f (xk)
k=0
R=
h3 12
nX1
f 00(⇠k)
=
k=0
(b a)3 12n2
f
00(⌘)
能源学院，航空航天中心
h= b a n
⌘ 2 [a, b]
Zb
f
( )d xx
⇡
( b
) a
Xn
(n)
Ck f
(xk
)
a
k=0
(b a)Ck(n) = Ak
称为 Newton-Cotes 求积公式。其中
(n)
Ck
=
1
Z
0
b
Yn
@
x
b a a i=0,i6=k xk
1
xi A d x
xi
称为
Newton-Cotes
系数。令
x
=
+ ，有 a th
t 2 [0, n]
数值积分
能源学院，航空航天中心
1
概述
求某函数、变量给定区间上的积分是经常会遇到的问题。
很多情况下无法利用 Newton-Leibniz 公式获得解析解，用数 值方法计算积分也能获得满足精度要求的结果。
数值积分从形式上看就是对节点处的函数值进行线性组合来 获得积分的近似解。
• 代数精度与插值型求积法 • Newton-Cotes 求积公式 • 复合求积公式，变步长梯形求积公式 • 龙贝格算法 • 高斯型求积公式
xk
= x0 + kh,
k
= 0 1 ··· ,,
,n
(2) 在第 k 个区间 [xk, xk+1] 上使用求积公式求得 Ik 。 nX1
(3) 求 n 个子区间积分之和 I = Ik ，作为整个区间积
分的近似值。
k=0
能源学院，航空航天中心
21
复合梯形求积公式
由梯形求积公式：
Ik =
h 2
f (xk) + f (xk+1)
的积分余项为零。
(n+1)( ) = ( + 1)! f xn
)
Rn(f
,
) x
=
(n+1)( ) f⇠ ( + 1)!
!n+1
() x
=
!n+1
() x
n
Zb
Z n Yn
Rn(f ) =
!n+1
( )d xx
=
n+2
h
( )d t jt
a
0 j=0
令 n = 2k ，t = u + k ，则 u 2 [ k, k]
并对其进行积分，这样获得的积分公式为插值型求积
公式。
Zb
对于积分 ( )d ，在区间 [a, b] 上给定 n+1 个节点 xi , i = 0, · · · , n fx x
a
和对应的函数值
f (xk)
，构造函数
() fx
的
Lagrange
插值多项式：
Xn
( )= fx
(
k=0 x
!n+1(x) xk)!n0 +1(xk
i
=
0 ,
1 ,
·
·
·
,
m
m+1个方程
i 次多项式的近似积分公式
能源学院，航空航天中心
6
代数精度─待定系数法
Xn
Zb
i+1 i+1
i
Ak xk
=
k=0
id = b xx
a
a +1 i
,
i
=
0 ,
1 ,
·
·
·
,
m
一一般而言，当给定代数精度后，待定系数法取 n=m ，并取m+1 个节
点。那么有 m+1 个方程，求解 m+1 个待定系数 Ak 。若 m+1 个节
待定系数法构造求积公式：以代数精度为出发点构造求积公式的方法。
要使得求积公式具有 m 次代数精度，只要使它对于 ( ) = i ( = 0 1 ··· )
f x x i , , ,m
都成立，即要求：
i 次多项式的积分精确公式
Xn
Zb
i+1 i+1
i
Ak xk
=
k=0
id = b xx
a
a +1 i
,
⇡b
6
a
() fa
+
4 f
✓
a
+ 2
b
◆
+
f
! () b
Simpson求积公式
能源学院，航空航天中心
16
Newton-Cotes求积─误差估计
定理7.3 当 n 为偶数时，n 阶 Newton-Cotes公式的代数精度至少是 n+1.
证明：
需要验证当 n 为偶数时，Newton-Cotes公式对于 ( ) = (n+1) fx x
能源学院，航空航天中心
17
Newton-Cotes求积─误差估计
Z k Y 2k
Rn(f ) = hn+2
(u + k j)du
k j=0
Y 2k
Yk
被积函数 H(u) = (u + k j) = (u
项
j=0
j= k
j) 有2k+1
( u j) = (u + j) = (u ( j)) H(u) 是奇函数。
+ 2
b
,
f
(
a
+ 2
b
))
的抛物线来
代替曲线
() fx
，得到
Simpson
积分公式
I
=
Zb
a
( )d fx x
⇡
b
6
a
✓ ()
fa
+
4 f
✓ a
+ 2
◆ b
+
◆ () fb
求积公式的一一般形式：
Zb
Xn
= I
( )d ⇡ fx x
Akf (xk)
a
k=0
能源学院，航空航天中心
5
代数精度
定义 如果某求积公式对于次数不超过 m 的多项式，均能精确成 立，但对于 m+1 次多项式不能精确成立，则称该求积公式有 m 次代数精度。
Xn
Zb
lim
n!1,h!0
Akf (xk) =
k=0
( )d fx x
a
稳定： 对于任意给定的 " > 0 ，若 9 > 0 ，当 |f (xk)
˜ fk
|
时
|In(f )
In(f˜)| =
Xn ⇣ Ak f (xk)
k=0
则称求积公式稳定。
⌘
˜ fk
"
Zb
Xn
定理7.2 若求积公式
( )d ⇡ fx x
a
k=0
Zb
其中 Ak =
a
( x
!n+1
() x
xk)!n0 +1
(xk
)
d x
是对多项式每一一项的积分，
Zb
Rn(f ) =
Rn
( f,
)d xx
是插值余项的积分，也是积分公式的误
a
差。略去误差项，就得到插值型求积公式：
对应的余项估计为
Zb
Xn
f
( )d xx
⇡
Akf (xk)
a
k=0
Zb
|Rn
Rn
( f
,
)d xx
，其中
a
Rn(f,
) x
=
(n+1)( ) f⇠ ( + 1)!
!n+1
() x
n
而对于函数
() fx
=
i( xi
=
0 ,
1 ,
·
··
) ,n
， f (n+1)
=
0
Rn(f ) = 0
插值型求积公式对次数小于 n+1 的多项式精确成立，故插 值型求积公式至少具有 n 次代数精度。
能源学院，航空航天中心
19
复合求积公式
能源学院，航空航天中心
20
复合求积公式
如果积分区间比较大，直接使用上述单步求积公式的话精度 难以保证（高阶插值不稳定）。一一般将积分区间 [a,b] 分成 n 个小区间，对 ( ) 进行分段低阶插值，然后计算积分近似
fx
值。
通常的步骤：
(1) 将求积区间 n 等分，步长为 h = b n a ，节点为
能源学院，航空航天中心
11
插值型求积公式─与待定系数法
充分性证明见课本227页
• 待定系数法用解线性方程组的方法来构造求积公式。其方程组
系数矩阵是病态矩阵，当节点数较多时求解困难。
• 插值型求积公式通过对插值多项式积分获得求积公式，计算稳
定，代数精度高。
能源学院，航空航天中心
12
收敛性与稳定性
收敛：
能源学院，航空航天中心
2
基本思想
计算定积分：
Zb
= ( )d
I
fx x
a
解析法（ Newton-Leibliz 公式）：
Zb
= ( )d = ( )
I
fx x Fb
a
() Fa
在解决实际问题的过程中，可能会遇到以下情况：
（1）函数
f
() x
的原函数
( ) 无法获得。 Fx
（2）函数
f
() x
的表达式未知。
Zb
在区间 [a,b] 上使用梯形公式及复合形式逐次计算积分
( )d fx x
T0 = b
a 2
f (a) + f (b)
a
此时步长 h0 = b a 。将区间 [a, b] 二等分，使用 n=2 时的复合
( f
)|
Mn+1 ( + 1)!
n
|!n+1(x)|dx
a
能源学院，航空航天中心
10
插值型求积公式─代数精度
定理7.1 证明：
n+1 个节点的求积公式为插值型的充要条件是该公式至少
有 n 次代数精度。
Zb
Xn
(必要性) 对于插值型求积公式
( )d ⇡ fx x
Akf (xk)
Zb
a
k=0
它的余项为 Rn(f ) =
⇡
A0f
(0)
+
A1f
(1)
+
A2f
(2)
+
A3f
(3)
0
代数精度和节点给定，只需对应以下公式求出系数 Ak
Xn
Zb
i+1
i+1
i
Ak xk
=
k=0
id = b xx
a
a +1 i
,
i
=
0 ,
1 ,
·
·
·
,m
A0 + A1 + A2 + A3 =3,
A1
+
2A2
+
3A3
=
9 2
,
A1 + 4A2 + 9A3 =9,
)
f
(xk
)
+
Rn(f
,
) x
Yn
其中
!n+1
() x
=
( x
i=0
是多项式插值余项。
xi)
，
Rn(f,
) x
=
(n+1)( ) f⇠ ( + 1)!
!n+1(x),
⇠
2( ) a, b
n
能源学院，航空航天中心
9
插值型求积公式
对以上
() fx
的表达式进行积分，得到：
Zb
Xn
( )d = fx x
Akf (xk) + Rn(f )
=
Z1
0
tdt
=
1 2
Zb
a
f
( )d xx
⇡
b
a 2
( )+ ( ) fa fb
梯形求积公式
当 n=2 时，Cotes 系数为
C0(2)
=
1 4
Z2 (t
0
1)(t
2)dt
=
1 6
C1(2)
=
1 2
Z2 t(t
0
2)dt
=
2 3
C2(2)
=
1 4
Z2 t(t
0
1)dt
=
1 6
Zb
a
( )d fx x