数值积分

其中a, b为小区间的端点，x2 , , xn1, A1, An 为2n-2个参数，
代数精度可达到2n-3
注意：实际计算中一般采用自适应方法确定步长
矩形 公式
Sum(x)
用MATLAB 作数值积分
n 1
Ln h
fk
k 0
n
Rn h
fk
k 1
输入数组x(即fk),输出x的和(数)
cumsum(x) 输入数组x,输出x的依次累加和(数组)
梯形 公式
n 1
Tn h
h fk 2 ( f0 fn )
k 1
trapz(x) 输入数组x,输出按梯形公式x的积分(单位步长)
trapz(x,y)
输入同长度数组 x,y,输出按梯形公式 y对x的积分(步长不一定相等)
用MATLAB 作数值积分
(
f0

fn)
(3)
k 1
数值积分 2.辛普森(Simpson)公式
（抛物线公式）
梯形公式相当于用分段线性插值函数代替 f (x)
提高精度
分段二次插值函数
抛物线 公式
y
y=f(x)
每段要用相邻两小区间
端点的三个函数值
f2k
f2k+
f2k+2
1
区间数必须为偶数 n 2m
o
a x2k x2k+1 x2k+2 b x (x2k , f2k ), (x2k1, f2k1), (x2k2, f2k2 )
高斯公式的思路
取消对节点的限制，按照代数精度最大
的原则，同时确定节点xk和系数Ak
1
对于 I f ( x)dx 1
构造求积公式
G2 A1 f (x1) A2 f (x2) 使G2的代数精度为3
1
f (x) 1, x, x2, x3 f (x)dx A1 f (x1) A2 f (x2 ) 1
f (x) T(x)
f
(k
2
)
(
x

xk
)(x

xk
1 ),
x,k (xk , xk1)
因为：(x-xk)(x-xk+1)在(xk,xk+1)不变号，所以：
xk1 [ f (x) T (x)]dx
f
(k
2
)
xk 1
(x

xk
)(x

xk
1)dx


同理可得：
|
R(
f
,
Sn )
|
h4 180
M 4 (b

a)
(6)
其中 M4 max f (4) (x) , x (a,b)
辛普森公式Sn的误差是h4阶的。
梯形公式和辛普森公式的收敛性
若对I某个数值积分In有
lim
n
I In hp
c（非零常数）
则称 In是 p 阶收敛的。
梯形公式 2 阶收敛，辛普森公式 4 阶收敛。
k 0,1,L ,m 1
2.辛普森(Simpson)公式（抛物线公式） 用(x2k , f2k ), (x2k1, f2k1), (x2k2 , f2k2 )构造 二次插值函数sk(x)
x2k 2
h
sk (x)dx
x2k
3 ( f2k
4 f2k 1
f2k 2)
对k求和(共m段)，得（复合）辛普森公式：
c=0: 至少p阶收敛（超p阶收敛）
积分步长的自动选取 选定数值积分公式后，如何确定步长h以满足给定的误差
梯形公式
I
Tn
h2 12
( f ' (b)
f ' (a))
h h 2 (n 2n)
1 I T2n 4 (I Tn )
1 I T2n 3 (T2n Tn )
代数精度为3 节点加密时，原计算信息无法利用
改进的高斯公式
思路：将积分区间分小，在小区间上用n不太 大 的 Gn 。而在节点加密一倍时能够利用原节点的函 数值，可以把区间的端点作为固定节点。
Gauss-Lobatto求积公式
n 1
Gn A1 f (a) Ak f (xk ) An f (b) k 2
数学建模讲座(三)
数值积分
胡支军 贵州大学数学系 贵州赛区组委会
数值 积分
为什么要作数值积分
• 积分是重要的数学工具，是微分方程、概率 论等的基础；在实际问题中有直接应用。
• 许多函数“积不出来”,只能用数值方法，如
b x2 e 2 dx , a
bsin x dx
ax
• 对于用离散数据或者图形表示的函数， 计算积分只有求助于数值方法。
m1
m1
Sm

h 3
(
f0

f2m

4
f2k 1 2
f2k ),
h ba 2m
(4)
k 0
k 1
梯形公式
的误差估计
n1
h
Tn

h
k 1
fk

2 ( f0

fn )

b
f (x)dx
a
b
R( f ,Tn ) I Tn f (x)dx Tn
a
梯形公式在每小段上是用线性插值函数T(x)代替 f(x)
确定 x1, x2, A1, A2
将f(x)代入计算得
A1 A2 2 A1x1 A2 x2 0 A1x12 A2 x22 2 / 3 A1x13 A2 x23 0
x1 1/ 3, x2 1/ 3, A1 A2 1
G2 f (1/ 3) f (1/ 3) 用n个节点，Gn的代数精度可达2n-1, 但是需解
y
y=f(x) a x0 x1 xk xn b,
h ba, n
fk f (xk )
n 1
fk
Ln h
fk
(1)
o
a
xk-1 xk xk+1 b x
k 0 n
Rn h
fk
(2)
Ln, Rn 平均，得到
k 1
n1
梯形公式
Tn h
fk

h 2
辛普森公式
m1
m1
Sn

h 3
( f0

f2m
4
f2k 1 2
f2k ),
h ba 2m
k 0
k 1
用自适应辛普森公式计算
quad(@fun,a,b,tol,trace) [I,fn]=quad(…)
tol为绝对误差，缺省时为10-6
n1
Gauss-Lobatto公式 Gn A1 f (a) Ak f (xk ) An f (b) k 2
长方体上计算三重积分的命令： triplequad(@fun,xmin,xmax,ymin,ymax, zmin,zmax,tol) 注：fun是被积函数，本身可以有自己的参数
用MATLAB 作数值积分

例. 计算 4
1
dx
0 1 sin x
1）矩形公式和梯形公式 将(0, /4)100等分
x a ~ 长半轴, b ~ 短半轴,由s1, s2, r决定
b
y a
s1=439km, s2= 2384km, r=6371km
2a 2r s2 s1
a r s2 s1 7782.5 2
oc
x 焦距c a r s1 c s2 s1
2
s2
r
s1
b a2 c2 7721.5
h3 12
n 1 k 0
f (k )
估计
| R( f
M2 max f (x) , x (a,b)
,Tn )
|
h2 12
M 2 (b

a)
(5)
因为n b a h
梯形公式Tn的
h2
误差是h2阶的
I Tn 12 ( f' (b) f' (a)) (5' )
辛普森公式的误差估计
h3 12
f (k )
xk
xk
I Tn
b

a
f ( x)dx Tn
h3 n1 12 k 0
f (k )
I Tn
h2
1 b 12 a
f '' ( x)dx 1 ( f ' (b) 12
f '(a))
梯形公式
的误差
|
R(
f
,Tn ) |
数值积分实例 人造卫星轨道长度

L 4 2 a2 sin2 t b2 cos2 t dt 0
用梯形公式和辛普森公式计算
Jifen2.m
轨道长度 L=4.8707104千米
只将区间5等分，梯形公式就给出很好的结果
数值积分的基本思路
回忆定积分的定义
b
I f (x)dx lim In,
a
n
n
In
f
(k
)
b
n
a
k 1
n充分大时In就是I的数值积分
各种数值积分方法研究的是
k 如何取值，区间 (a,b)如何划分，
使得既能保证一定精度，计算量又小。 （计算功效：算得准，算得快）
数值积分 1.从矩形公式到梯形公式
复杂的非线性方程组，实用价值不大。
常用的高斯公式
将(a,b)分小，把小区间变换为(-1，1), 再用G2
b
m
f (x)dx h 2
[ f (zk(1) ) f (zk(2) )]
a
k 1
zk(1)

xk
xk 1 2

h 23
,
zk(2)

xk
xk 1 2

h 23
h (b a) / m, xk a kh, k 0,1, m
用二分法只要 T2n Tn
I T2n
且T2n可在Tn
基础上计算
n 1
T2n

Tn 2

h 2
fk 12
k 0
其中fk+1/2是原分点xk,xk+1的中点（记xk+1/2）的函数值
高斯(Gauss)求积公式
矩形公式(1)、(2) 梯形公式(3)
辛普森公式(4)
n
In
2）辛普森公式和Gauss-Lobatto公式 精确、方便
无法计算用数值给出的函数的积分
精确值为 2
Jifen1a.m Jifen1b.m
数值积分的应用
y
s2
s1
实例
人造卫星轨道长度
r o
轨道长度
x 近地点s1=439km,远地点s2= 2384km 地球半径r=6371km a ~ 长半轴
x a cost, y b sint b ~ 短半轴
Ak f (xk ) (7) Newton
k 1
Ak是与f无关的常数
-Cotes 方法
代数 精度
b
设 f (x) xk , 用(7)计算 I f ( x)dx, a
若对于 k 0,1, m 都有 In I,
而当 k m 1, In I, 则称In的代数精度为m.
梯形公式的代数精度（考察T1）
k=1
T1

h[ 2
f
(a)
f
(b)]
(b a)(a b) 2
f(x)=x
I
b
xdx
b2
a2
a
2wk.baidu.com
T1 I
k=2 f(x)=x2
T1

(b a)(a2 2
b2)
I
b
x2dx
b3

a3
a
3
T1 I
梯形公式的代数精度为1 辛普森公式的代数精度为3
quadl(@fun,a,b,tol,trace) 用自适应Gauss-Lobatto公式计算
[I,fn]=quadl(…)
tol为绝对误差，缺省时为10-6
注意：fun.m中应以自变量为矩阵的形式输入(点运算)
广义积分、二重和三重积分
广义积分： 通过分析和控制误差，转换成普通积分
向量值积分： quadv(@fun,a,b,tol,trace) 矩形域上计算二重积分的命令： dblquad(@fun,xmin,xmax,ymin,ymax,tol)
(0 t 2 )
由s1, s2 , r决定


L 4 2 x2 y2 dt 4 2 a2 sin2 t b2 cos2 tdt
0
0
需要作数值积分
数值积分实例 人造卫星轨道长度
y
s2
r o

s1
L 4 2 a2 sin2 t b2 cos2 t dt 0