若干数值积分的计算方法
黄海琼
（广西民族大学数计学院04数本1班 南宁 530006）
摘 要: 本文讨论了若干数值积分的计算方法。在一维情形下，介绍了Newton-Cotes 公式，Gauss 型等求积法则； 在二维情形下, 主要介绍了二元Newton-Cotes 积分方法。最后，对几类数值积分方法及其数值实验进行比较评述。
关键词: 牛顿-柯特斯公式；Gauss 型求积法则；二元数值积分；数值实验
Some Computational Methods of numerical integration
Huang Haiqiong
(College of Mathematics and Computer Science,Guangxi University for Nationalities, Nanning 530006)
Abstract: In this paper, some computational methods about numerical integration are discussed. under the univariate situation, the quadrature rule of Newton-Cotes formula, Gauss formula and so on is introduced. Under the two-dimensional situation, it mainly introduced the dual Newton-Cotes integral method. Finally, the numerical integration methods and numerical experiment were discussed.
Key word: Newton-Cotes formula; Gauss integration principle; dual numerical integration; numerical experiment.
1 引 言
数值积分是积分计算的重要方法，是数值逼近的重要内容，是函数插值的最直接应用，也是工程技
术计算中常常遇到的一个问题[1]。在应用上，人们常要求算出具体数值，因此数值积分就成了数值分析的一个重要内容。在更为复杂的计算问题中，数值积分也常常是一个基本组成部分。
在微积分理论中，我们知道了牛顿-莱布尼茨(Newton-Leibniz )公式
()()()b
a
f x dx F b F a =-?
其中()F x 是被积函数()f x 的某个原函数。但是随着学习的深入，我们发现一个问题: 对很多实际问题，上述公式却无能为力。这主要是因为:它们或是被积函数没有解析形式的原函数,或是只知道被积函数在一些点上的值,而不知道函数的形式,对此,牛顿—莱布尼茨(Newton-Leibniz)公式就无能为力了。此外,即使被积函数存在原函数,但因找原函数很复杂,人们也不愿花费太多的时间在求原函数上,这些都促使人们寻找定积分近似计算方法的研究,特别是有了计算机后,人们希望这种定积分近似计算方法能在计算机上实现,并保证计算结果的精度,具有这种特性的定积分近似计算方法称为数值积分。
2 一维数值积分的计算方法
2.1 一维数值积分方法的基本思想[2]
对于一维数值积分法的思想来源于定积分的定义，即
1
()lim ()n
b
i i a
i f x dx f x x λ→==?∑?
其中
max{}i i
x λ=?，一般的提法是: 对于给定的权函数()0,[,]x x a b ρ≥∈，用()f x 在点
01n a x x x b =<
00110
()()()()()n
n n n i i i I f A f x A f x A f x A f x ==++??????+=∑
作为积分()()()b
a
I f x f x dx ρ=
?
的近似值，即
()()()()n
b
i i a
i I f x f x dx A f x ρ==≈∑?
(2.1)
并称此为求积公式，称
()()()n R f I f I f =-
为求积公式(2.1)的余项或误差，i x 及(0,1,2,,)i A i n =???分别称为求积公式(2.1)的求积节点及求积系数，这里求积系数(0,1,2,,)i A i n =???只与权函数()x ρ及积分区间[,]a b 有关，而与()f x 无关。
为保证机械求积公式的精度, 自然希望它对尽可能多的简单函数是准确的，即要求它对一切m 次多项式是准确的, 而对1m +次多项式不一定准确。则得到关于系数i A 的1n +阶线性方程组：
0122
001111
0011
21
n n n n n n n n n n
A A A b a b a A x A x A x b a A x A x
A x
n +++++=-??-?+
++
=
????-+++=
?+?
L L
M M
M M
L
由于系数行列式为Vandermonde 行列式,不为零,则解i A 是唯一存在的。
定义2.1.1 如果当()(0,1,,)k
f x x k m ==L 时，求积公式(2.1)精确成立，而当1
()m f x x +=时，
求积公式(2.1)不精确成立，那么称求积公式(2.1)具有m 次代数精度。
定理2.1.1 任意给定1n +个节点01,,,n x x x L ，如果()f x 是次数不超过n 的多项式，那么一定
存在常数(0,1,,)i A i n =L ，使求积公式(2.1)精确成立，即
()()()n
b
i i a
i x f x dx A f x ρ==∑?
.
证明 设()n P x 是()f x 关于节点01,,,n x x x L ，的n 次Lagrange 插值多项式, 即
()()()()()()n
n n i i n i f x P x r x l x f x r x ==+=+∑
其中()
(),(0,1,,)()()i i i x l x i n x x x ωω=
='-L 是Lagrange 基函数，0
()()n
j j x x x ω==-∏，()n r x 是Lagrange 插值余项。于是
()()()()()()()n
b
b
b
i i n a
a
a
i x f x dx x l x f x dx x r x dx ρρρ==+∑?
??
不妨令
()(),(0,1,,)b
i i a
A x l x dx i n ρ=
=?
L (2.2)
则有
()()()()()n
b
b
i i n a
a
i x f x dx A f x x r x dx ρρ==+∑?
?
因为()f x 是次数不超过n 的多项式，所以()()n f x P x ≡，这意味着0n r ≡，于是()()0b
n a
x r x dx ρ=?
，
从而
()()()n
b
i i a
i x f x dx A f x ρ==∑?
.
这个定理告诉我们，具有一定代数精度的求积公式是存在的。
定理2.1.2 形如(2.1)的求积公式的代数精度d n ≥的充要条件是:它是插值型求积公式。 证明 充分性上面已证. 现在来证必要性. 设求积公式
()()()n
b
i i a
i x f x dx A f x ρ=≈∑?
(2.3)
的代数精度d n ≥. 因为Lagrange 基函数(),(0,1,,)j n l x P j n ∈=L , 所以
()()(),(0,1,,)n
b
j i j i j a
i x l x dx Al x A j n ρ====∑?
L
故求积公式(2.3)是插值型求积公式。
定义2.1.2 如果01,,,n x x x L 属于区间[,]a b ，那么称
()()()n
b
i i a
i x f x dx A f x ρ=≈∑?
为内插型求积公式，其中求积系数(0,1,,)i A i n =L 由(2.2)决定。
2.2 几种常用数值积分方法
2.2.1 插值型求积公式[3]
在[,]a b 上,用以i x 为节点的n 次Lagrange 插值多项式()n L x 作为()f x 的逼近函数,即可得到插值型求积公式:
()()()()n
b
b
b
n i i a
a
a
i f x dx L x dx f x l x dx =≈=∑?
??
即
()()n
b
i i a
i f x dx A f x =≈∑?
其中01110111()()()()()
()()()()()()
b
b
i i n i i a
a
i i i i i i i n x x x x x x x x x x A l x dx dx x x x x x x x x x x -+-+-----=
=-----?
?
L L L L .
插值型积分公式具有n 次代数精度,且0i A >时公式是稳定的。 2.2.2 一元Newton- Cotes 公式
[4]
将区间[,]a b 等分成n 份，步长()
b a h n
-=, 求积节点为,(0,1,2,,)i x a ih i n =+=L 令x a th =+,则Lagrange 插值基函数为
0,0,,(0,1,2,,)n
n
i j i i j i i j j i x x t i
l j n x x j i
=≠=≠--===???--∏∏
求积系数j A 可表示为
00,(1)()(),(0,1,2,,)!()!n j n
b
n j j a i i j
h A l x dx t i dt j n j n j -=≠-==-=???-∏?? 令
00,(1)()!()!n j
n n j
j i i j
A C t i dt b a n j n j -=≠-==--?-∏?
称j C 为Cotes 系数，则求积公式可化为
(1)00
()()()()(),(,)(1)!n n
n
b
b j j i a
a j i f f x dx
b a C f x x x dx a b n ηη+===-+-∈+∑∏?
? 若令()1f x ≡,可得出
1n
j
j C
=≡∑. 记
(1)(1)0
()()(),()[,](1)!n n b n n i a i f R f x x dx f x C a b n η++==-∈+∏? (2.4) 称(2.4)为Newton-Cotes 公式的截断误差。Cotes 系数j C 与被积函数和积分区间都无关,只要给出区间等分数n 即可求出。
由于n 点Newton-Cotes 公式是特殊的插值型求积公式,故其代数精度至少是1n -,但由于其节点的等距性,我们有
定理2.2.1 当求积节点个数n 为奇数时,对应的Newton-Cotes 求积公式的代数精度至少为n -1. 证明 考虑求积余项. 由于是插值型求积公式,故有
()0,0,1,,1,k R x k n ==-L
对n
x 有
()
1
()()()!
n n n
b
x n
k
a
k x R x x x dx n ξ
===-∏?
11
1
(1)n
n n k h
t k dx -+==-+∏? (x a th =+换元)
21
1
1
(1)l l n l k h
s l k dx ++-==+-+∏? (令,21s t l n l =-=+)
记211
()(1)(),11l l
k m l
s s l k s m m k ?+==-=
+-+=+=-+∏∏, 则
21
()()(1)
()()l
l
l m l
m l
s s m s m s ??+=-=--=-+=-+=-∏∏
故()s ?是奇函数,得()0n
R x =,定理得证.
定理2.2.2 2n 阶Newton-Cotes 公式具有21n +次代数精度.
证明 设21
210
()n i
n i i P x a x
++==
∑为21n +次多项式，由(2.4)有
(21)2212210
()
()()(21)!n n
b
n n n i a
i P R P x x dx n η+++==-+∏?
2210
()n
b n i
a
i a x x dx +==-∏?
在2n 阶Newton-Cotes 公式中，[,]a b 是分成2n 等分的，所以()2a b c +=必是分点.令()
2b a h n
-=， 于是
221211
1
()()()()n
n
b
n n n a
i i R P a x c x c ih x c ih dx ++===----+∏∏?
令2
a b
x c x u +-=-
=,则有 2222
221211
2
()()b a n
b a
n n n i R P a
u u i h du --++-==-∏?
222211
(){()}n
b
n a
i a x c x c i h dx +==---∏?
此积分的被积函数为奇函数，积分区间关于原点对称，故其值为零，即221()0n n R P +=. 定理得证.
下面我们给出一些常见的Newton-Cotes 公式及其余项: 1、令1n =, 即得梯形公式
1
1010011
(1),22C t dt C tdt =--===?? ,
()[()()]2
b a
b a f x dx f a f b -≈+? 当2
()[,]f x C a b ∈时,
3
()(),(,)12
f R f h a b ηη''=-
∈. 2、令2n =, 即得Simpson 公式
()[()4()()]62
b
a
b a a b
f x dx f a f f b -+≈
++?
， 当4
()[,]f x C a b ∈时,
3(4)
()(),(,)90
h R f f a b ηη=-∈.
3、令4n =, 即得Cotes 公式
01234()[7()32()12()32()7()]90
b
a
b a
f x dx f x f x f x f x f x -≈
++++?
当6
()[,]f x C a b ∈ 时，
7(6)
8()(),(,)945
R f h f a b ηη=-
∈.
2.2.3 复化公式
由定积分知识,定积分只与被积函数和积分区间有关,而在对被积函数做插值逼近时,多项式的次数越高,对被积函数的光滑程度要求也越高,且会出现Runge 现象。如7n >时,Newton-Cotes 公式就是不稳定的。因而,人们把目标转向积分区间,类似分段插值,把积分区间分割成若干小区间,在每个小区间上使用次数较低的Newton-Cotes 公式,然后把每个小区间上的结果加起来作为函数在整个区间上积分的近似,这就是复化的基本思想。常用的复化公式有:
1、复化梯形公式
1
1
()[()2()()]2n b
n i a
i h
f x T f a f x f b -=≈=++∑?
,
2
()()(),(,)12
T b a R f h f a b ηη-''=-
∈. 2、复化Simpson 公式
1
21211
()[()4()2()()]3m m b
m i i a
i i h
f x dx S f a f x f x f b --==≈=+++∑∑?
,
5(4)4(4)
(2)()()(),(,)2880180
S h b a R f mf h f a b ηηη-=-=-∈.
2.2.4 逐次分半技术与Romberg 公式
如何确定适当的n 使得近似值与真值之差在允许范围内,一般来说是比较困难的。而逐次分半技术是在求积过程中根据精度的要求,自动确定n 的选择是否满足精度要求,以二分后前后两次之差来估计误差,这样既缩小了步长,又能保留原有的计算结果,减少计算量。
对复化公式,,n n n T S C , 可得到
221
()3n n n I T T T ≈+-,
221
()15n n n I S S S ≈+-,
221
()63
n n n I C C C ≈+-,
且
24133n n n S T T =-,
2161
1515n n n C S S =-,
2641
6363
n n n R C C =-.
称之为Romberg 公式. Romberg 公式的加速效果是极其显著的,在相同精度要求下,计算量较小. 2.2.5 Gauss 型求积公式
为进一步提高求积公式的代数精度,可通过适当选择插值节点和求积系数,使得代数精度最高达到
21n +.把求积节点和求积系数视为同等参数求解，既可利用方程组得到,也可借助正交多项式的零点来确定。一般可设n 个节点的求积公式为
1
()()()()n
b
i i n a
i x f x dx A f x I f ρ=≈=∑?
. （2.5）
定理2.2.3 插值型求积公式（2.5）具有21n -次代数精度，必须且只须插值节点01,,,n x x x L 是
[,]a b 上以()x ρ为权的n 次正交多项式的零点.
证明 设插值型求积公式（2.5）具有21n -次代数精度，记n 次多项式1()()()n x x x x x ω=--L ,对任意给定的多项式1()n q x P -∈. 因为（2.5）具有21n -次代数精度，()()x q x ω是21n -次多项式，所以
1
()()()()()0n
b
i i i a
i x x q x dx A x q x ρωω=≈=∑?
即()x ω与()q x 关于权()x ρ正交.
反之,设()x ω与任意次数1n ≤-的多项式关于权()x ρ正交,任意给定21()n f x P -∈恒可表示为
()()()()f x x q x r x ω=+,
其中()q x 与()r x 均属于1n P -,且
()()()()()()()b
b b
a
a
a
x f x dx x q x x dx x r x dx ρρωρ=+?
??
1
()()()n
b
i i a
i x r x dx A r x ρ===∑?
1
1
[()()()]()n n
i i i i i i i i A q x x r x A f x ω===+=∑∑.
于是定理得证.
定义2.2.1 具有最高代数精度为21n -的插值型求积公式称为Gauss 型求积公式，求积节点
01,,,n x x x L 称为Gauss 点.
定理2.2.4 求积公式(2.1)是Gauss 型的，当且仅当Gauss 点01n a x x x b <<<< 证明 充分性定理2.2.3已证.下面证必要性. 令 101()()()()n n x x x x x x x ω+=---L 任取次数n ≤的多项式()q x 有 110 ()()()()()0n b n i n i i a i x x q x dx A x q x ρωω++===∑? 即1(,)0n q ω+=对一切次数不高于n 的代数多项式()q x 成立，从而1()n x ω+是[,]a b 上关于权()x ρ的 1n +次正交多项式. Gauss 点,0,1,,i x i n =L 是1n +次正交多项式1()n x ω+的根. 证毕. 定理2.2.5 若（2.5）为Gauss 型内插求积公式，则其求积系数i A 皆为正，且 2 2 2 1 ()[()][()]() b i a i i x x A dx x x x ρωω= '-? (2.6) 其中1()()()n x x x x x ω=--L . 证明 显然,多项式2() ()( )i x f x x x ω=-的次数不超过22n -,且 2 0, ()[()],k i k i f x x k i ω≠?=?'=? 又由假设知（2.5）是Gauss 型求积公式，所以 21 ()()()[()]n b k k i i a k x f x dx A f x A x ρω='==∑? 注意到()0f x >,所以 ()()0b a x f x dx ρ>? ,于是得 2 1 ()()0,(1,,)[()]b i a i A x f x dx i n x ρω= >='? L 再将()f x 的表达式代入上式，便可得到公式(2.6). 证毕. 定理2.2.6 如果()[,]f x C a b ∈,则是Gauss 型求积公式收敛,即 lim ()()()b n a n I f x f x dx ρ→∞ =?. 定理2.2.7 Gauss 型求积公式是数值稳定的。且对有限区间上的连续函数，Gauss 求积的数值随节点数目的增加而收敛到准确积分值. 定理2.2.8 若(22) [,]()n a b f x C +∈，则Gauss 型求积公式(2.5)的余项为 (22)2 101()[]()[()()()],(,)(22)! n b n n a f R f x x x x x x x dx a b n ξρξ++=?---∈+?L (2.7) 证明 取()f x 的Hermite 插值多项式为21()n H x +，满足插值条件 212 1()(),()(),0,1,,n k k n k k H x f x H x f x k n ++''===L 由于 210 ()()()()n n b n k k k k a k k x H x dx A H x A f x ρ+====∑∑? 可得 121[]()[()()]b n n a R f x f x H x dx ρ++=-? (22)201() [()()()](22)! n b n a f x x x x x x dx n ξ+=---+? L 利用2 01[()()()]0n x x x x x x ---≥L 及积分中值定理即得式(2.7). 证毕. 利用具有不同权函数的正交多项式, 就能得到不同类型的Gauss 型求积公式: 1、Gauss-Legendre : 21 21112(1) ()();[()]n i i i i n i xi f x dx A f x A nL x -=--≈=∑?. 2、Gauss-Laguerre : 2 2 1 (!)()();[()]n x i i i i i n i n e f x dx A f x A x L x ∞ -=≈='∑? . 3、Gauss-Hermite : 2 12 1 2()();[()]n n x i i i i n i n e f x dx A f x A H x +∞ --∞ =≈='∑? . 4、Gauss-Chebyshe ： 1 1 21();cos ;2(1)(1) n i i i i i i A f x x A n n π π-=+≈==++∑? . 利用上述几种公式,一维数值积分的近似计算问题已成功获得解决。 3 二元Newton-Cotes 积分方法 对于二元数值积分 21()() ((,))b f y a f y I g x y dx dy =? ? 的计算问题, 我们在本节做些讨论。 采用和一元数值积分类似的方法, 可构造出二元数值积分的Newton-Cotes 公式。 首先, 对积分区域进行简化. 由于12()()f y f y ，为非规则边界, 不好找到原函数, 于是把积分区域等分成若干个规则的矩形(如图1) 图1 积分区域 则每个矩形区域内的积分变成 ((,))b d a c g x y dx dy ?? ，故 21()() 00 ((,))((,))j b i j i N M b f y b d a f y a c j i I g x y dx dy g x y dx dy === =∑∑?? ? ? （3.1） 其中,M 表示等分纵坐标范围[,]a b 所得到的曲边矩形的条数,j N 表示在第j 个曲边矩形中等分横坐标范围12[()()]f x f y ,所得到的小矩形个数.可见,二元数值积分转换成一些在j j a y b ≤≤,i i c x d ≤≤矩形区域内的二元数值积分的二重累加和,但是,依然是二元数值积分,还需要进行改进。对单个的矩形区域内的二元数值积分 (,))b d a c g x y dx dy ? ?（，要尽量做到 * (,))()()b d b d a c a c g x y dx dy I y dy I x dx =? ???（ 从而实现一元化。 解决积分（3.1）的数值方法有很多, 常用插值多项式构造数值积分公式的方法. 具体做法如下: 给定区域,a y b c x d ≤≤≤≤上一组节点0101m n a y y y b c x x x d ≤<<<=≤<<<=L L ，以及节点处的(,)0,1,,,0,1,,i j g x y i m j n ==L L ，, 作(,)g x y 的插值多项式 00 (,)()()(,)(,)m n i j i j i j g x y I x I y g x y r x y ===+∑∑, 其中 0;0;(),()m n k k i j k k i k k j i k j k x x y y I x I y x x y y =≠=≠--==--∏∏ 代入得 (,))b d a c g x y dx dy ? ?（ 00 [()()(,)(,)]m n b d i j i j a c i j I x I y g x y r x y dxdy ===+∑∑? ? 00()()(,)(,)]m n b d b d i j i j a c a c i j I x I y g x y dxdy r x y dxdy ===+∑∑? ? ? ? 00 (,)()()(,)]m n b d b d i j j i a c a c i j g x y I y dy I x dx r x y dxdy ===+∑∑??? ? 记 ()()()(()(),(),,[](,)]d b b d m n mn m n i i j j ij i j c a a c A I x dx A I y dy A A A R g r x y dxdy ====??? ? ） 则 () () 00 00 ((,))(,)[](,)m n m n b d mn mn ij i j ij i j a c i j i j g x y dx dy A g x y R g A g x y =====+≈∑∑∑∑?? 在等距节点i x x c i h =+?, 其中步长x d c h m -= , 如果作变量替换x x c t h =+?,那么 0 (1)(1)(1)() ()(1)!()! m i x m i t t t i t i t m A m h dt i m i ---+---=--? L L 0(1)()(1)(1)(1)()!()! m i m i d c A m t t t i t i t m dt m i m i ---=--+----?L L 记 ()0(1)(1)(1)(1)(),0,1,,!()! m i m m i C t t t i t i t m dt i m mi m i --=--+---=-?L L L L 则 ()()()m m i i A d c C =-, 其中，() m i C 称为柯特斯系数, 部分柯特斯系数列在表1中。 同理 ()()()n n j j A b a C =-， 故 () ()()()()()()()()mn m n m n ij i j i j A d c C b a C d c b a C C =--=--， 于是 ()()00(,))()()(,)[]m n b d m n i j i j a c i j g x y dx dy d c b a C C g x y R g ===--+∑∑? ?（ ()()00 ()()(,)m n m n i j i j i j d c b a C C g x y ==≈--∑∑ （3.2） (3.2)称为二元数值积分的Newton-Cotes 公式[7] 。 在1m n ==时 1 (,))()()[(,)(,)(,)(,)]4 b d a c g x y dx dy d c b a g c a g c b g d a g d b ≈--+++??（. 在2m n ==时 1(,))()()[(,)4(,)(,)4(,)3622 b d a c a b c d g x y dx dy d c b a g c a g c g c b g a ++≈--+++??（ 16(,)4(,)(,)4(,)(,)]2222 c d a b c d a b g g b g d a g d g d b +++++++++. 表1 柯特斯系数 4 几类数值积分方法的比较评述 由于实际问题计算时要考虑到精度要求和计算量, 因而不同求积公式各有特点. 梯形求积公式和Simpson 求积公式虽然计算简单、使用方便, 但是精度较差, 但对于光滑性较差的被积函数有时比高精度方法更为有效。尤其梯形公式对被积函数是周期函数的效果更为突出。7n > 时, Newton-Cotes 公式是不稳定的, 而复化梯形公式和复化Simpson 公式在保留了低阶公式的优点, 又能获得较高的精度, 因此在实际计算中应用的最为广泛。 利用二分技术得到的Romberg 方法的算法简单, 易于编程实现。当节点加密提高积分近似程度时, 前面计算的结果可以为后面所用, 对减少计算量很有好处, 并有比较简单的误差估计, 能得到若干积分序列, 如果在做收敛性控制时, 同时检查各行、各列, 对于不同性态的函数可以用其中最快的收敛序列来逼近积分。 高斯型求积公式的精度较高, 特别对计算无穷积分和瑕积分方面是其它方法所不能比拟的。但因为 其节点不规则, 当节点增加时, 前面计算的函数值不能为后面所用, 上机计算时还要预先存入不同N 值的节点值和系数表, 比较麻烦。 5 数值积分方法的比较数值实验 本节通过Matlab上机实验的数值例子, 比较各种数值积分方法的优缺点。 (1)积分节点数相同, 不同算法的精度比较。 在理论上可以确定各种方法的精度高低, 但是不利于理解。运用目前流行的数学软件的强大功能进行上机操作, 能够非常准确、直观地认识到不同数值积分方法的准确程度。 例如: 对于积分 1 2 1 1 I dx x = + ?, 精度要求为6 10-.分别用梯形公式、Simpson公式、Gauss公式与 精确解进行比较。 利用Matlab编写程序运行后结果显示: [ Precision Trapezoidal Simpson Gauss ] [0.7853981635 0.7849814972 0.7853981535 0.7853981635] 显然, 梯形公式精度较低, Simpson公式较好, 而Gauss求积公式的精度最高, 在给定精度下与精确解(Precision)一致。 (2)同一算法,积分区间有不同等分数的比较。 例如:求积分 6 2 1 1 1 I dx x = + ?,采用梯形公式考虑同一算法的不同步长的比较或积分区间有不同等 分数时的计算精度、计算时间的比较,用Matlab编程,运行结果显示: 区间数目: [10 20 40 80 ] 计算结果: [0.6305 0.6228 0.6209 0.6204] (精确解为: 0.6203) 计算时间: [0 0 0 0.0100] 图形演示: 参考文献: [1] 苏成,韩大建. 正交各向异性弹性力学平面问题的样条虚边界元法[J].应用数学和力学,2002,23(4): 402-403. [2] 刘鹏飞,徐乃楠. 数值积分方法的比较教学研究与实验[J]. 长春师范学院学报(自然科学版),2007, 26(6):23. [3] 洪柏娣. 三次Hermite插值的扩展及其应用[J]. 福建电脑,2008(2):97-98. [4] 安梅. Newton-Cotes求积公式新证[J]. 河北大学学报(自然科学版).1994,14(4):67-68. [5] 包玉兰. 复化Newton-Cotes及其误差估计[J]. 内蒙古民族师院学报(自然科学版).1997,12(2)：153-137. [6] 曹丽华. 一类广义Guass型求积公式[J]. 数学物理学报.2007,27(3):524-525. [7] 陈付龙. 二元数值积分的计算方法[J]. 计算工程与应用.2007,43(19):32-33. 几种定积分的数值计算方法 摘要:本文归纳了定积分近似计算中的几种常用方法,并着重分析了各种数值方法的计 算思想,结合实例,对其优劣性作了简要说明. 关键词:数值方法;矩形法;梯形法;抛物线法;类矩形;类梯形 Several Numerical Methods for Solving Definite Integrals Abstract:Several common methods for solving definite integrals are summarized in this paper. Meantime, the idea for each method is emphatically analyzed. Afterwards, a numerical example is illustrated to show that the advantages and disadvantages of these methods. Keywords:Numerical methods, Rectangle method, Trapezoidal method, Parabolic method, Class rectangle, Class trapezoid 1. 引言 在科学研究和实际生产中,经常遇到求积分的计算问题,由积分学知识可知,若函数 )(x f 在区间],[b a 连续且原函数为)(x F ,则可用牛顿-莱布尼茨公式 ?-=b a a F b F x f ) ()()( 求得积分.这个公式不论在理论上还是在解决实际问题中都起到了很大的作用. 在科学研究和实际生产中,经常遇到求积分的计算问题,由积分学知识可知,若函数)(x f 在区间],[b a 连续且原函数为)(x F ,则可用牛顿-莱布尼茨公式 ?-=b a a F b F x f ) ()()( 求得积分.这个公式不论在理论上还是在解决实际问题中都起到了很大的作用.另外,对于求导数也有一系列的求导公式和求导法则.但是,在实际问题中遇到求积分的计算,经常会有这样的情况: (1)函数)(x f 的原函数无法用初等函数给出.例如积分 dx e x ?-1 02 , ? 1 sin dx x x 等,从而无法用牛顿-莱布尼茨公式计算出积分。 (2)函数)(x f 使用表格形式或图形给出,因而无法直接用积分公式或导数公式。 (3)函数)(x f 的原函数或导数值虽然能够求出,但形式过于复杂,不便使用. 由此可见,利用原函数求积分或利用求导法则求导数有它的局限性,所以就有了求解数值积分的很多方法,目前有牛顿—柯特斯公式法,矩形法,梯形法,抛物线法,随机投点法,平均值法,高斯型求积法,龙贝格积分法,李查逊外推算法等等,本文对其中部分方法作一个比较． 2.几何意义上的数值算法 s 在几何上表示以],[b a 为底,以曲线)(x f y =为曲边的曲边梯形的面积A ,因此,计 算s 的近似值也就是A 的近似值,如图1所示.沿着积分区间],[b a ,可以把大的曲边梯形分割成许多小的曲边梯形面积之和.常采用均匀分割,假设],[b a 上等分n 的小区间 ,x 1-i h x i +=b x a x n ==,0,其中n a b h -= 表示小区间的长度. 2.1矩形法 实验二数值方法计算积分 学号：姓名：指导教师：实验目的 1、了解并掌握matlab软件的基本编程、操作方法； 2、初步了解matlab中的部分函数，熟悉循环语句的使用； 3、通过上机进一步领悟用复合梯形、复合辛普森公式，以及用龙贝格求积 方法计算积分的原理。 一、用不同数值方法计算积分 10x ln xdx=-94. (1)取不同的步长h.分别用复合梯形及辛普森求积计算积分，给出误差中关 于h的函数，并与积分精确值比较两个公式的精度，是否存在一个最小 的h，使得精度不能再被改善？ (2)用龙贝格求积计算完成问题（1）。 二、实现实验 1、流程图： 下图是龙贝格算法框图： 2、 算法： (1) 复合梯形公式：Tn=++)()([2b f a f h 2∑-=1 1 )](n k xk f ； (2) 复合辛普森公式：Sn=6h [f(a)+f(b)+2∑-=11)](n k xk f +4∑-=+1 )2/1(n k x f ]； 以上两种算法都是将a-b 之间分成多个小区间（n ）,则h=(b-a)/n,x k =a+kh, x k+1/2=a+(k+1/2)h,利用梯形求积根据两公式便可。 (3) 龙贝格算法：在指定区间内将步长依次二分的过程中运用如下公式 1、Sn= 34T2n-31 Tn 2、 Cn=1516S2n-151 Sn 3、 Rn=6364C2n-631 Cn 从而实现算法。 3、 程序设计 （1）、复合梯形法： function t=natrapz(fname,a,b,n) h=(b-a)/n; fa=feval(fname,a);fb=feval(fname,b);f=feval(fname,a+h:h:b-h+0. 001*h); t=h*(0.5*(fa+fb)+sum(f)); （2）、复合辛普森法： function t=natrapz(fname,a,b,n) h=(b-a)/n; fa=feval(fname,a);fb=feval(fname,b);f1=feval(fname,a+h:h:b-h+0 .001*h); f2=feval(fname,a+h/2:h:b-h+0.001*h); t=h/6*(fa+fb+2*sum(f1)+4*sum(f2)); （3）龙贝格法： function [I,step]=Roberg(f,a,b,eps) if(nargin==3) eps=1.0e-4; end; M=1; tol=10; k=0; T=zeros(1,1); h=b-a; T(1,1)=(h/2)*(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)), 实验二数值积分实验 一. 实验目的 （1）熟悉数值积分与数值微分方法的基本思想，加深对数值积分与数值微分方法的理解。 （2）熟悉Matlab编程环境，利用Matlab实现具体的数值积分与数值微分方法。 二. 实验要求 用Matlab软件实现复化梯形方法、复化辛甫生方法、龙贝格方法和高斯公式的相应算法，并用实例在计算机上计算。 三.实验内容 1. 实验题目 已知x e x f x4 sin 1 ) (- + =的数据表 分别编写用复化梯形法、复化辛甫生公法、龙贝格法、三点高斯法求积分?=1 ) (dx x f I 近似值的计算机程序。 A.复化梯形法: a．编写文件Trapezoid.m,代码如下所示： b．编写文件f2.m: c.运行： B．复化辛甫生公法 a．编写文件FSimpson.m,代码如下所示： b．编写文件f2.m: function f=f2(x) f=1+exp(-x).*sin(4*x); c.运行： C．龙贝格法 a．编写文件Romberg.m,代码如下所示： b.运行： D.三点高斯法 a.编写文件TGauss.m文件，如下所示： b．运行： 2. 设计思想 要求针对上述题目，详细分析每种算法的设计思想。 总体的思想是化复杂为简单的重复 A．复化梯形法使用直接法，通过递归，缩减规模； B．复化辛甫生也是使用直接法，根据公式直接进行编程，通过递归缩减规模； C．龙贝格算法应该在做了的几个中最体现了“化复杂为简单的重复”的思想，多个循环通过变量的适当递增，和一个for循环语句来实现，循环主体只有一句话，但确是整个程序中的亮点和难点； D．三点高斯法直接通过一条简单的公式来编写程序，难度不大； 四.实验体会 对实验过程进行分析总结，对比不同方法的精度，指出每种算法的设计要点及应注意的事项，以及自己通过实验所获得的对数值积分方法的理解。 毕业论文文献综述 信息与计算科学 导数的数值计算方法 一、 前言部分 导数概念的产生有着直觉的起源，与曲线的切线和运动质点的速度有密切的关系．导数用于描述函数变化率，刻画函数的因变量随自变量变化的快慢程度．比如说，物理上考虑功随时间的变化率（称为功率），化学上考虑反应物的量对时间的变化率（称为反应速度），经济学上考虑生产某种产品的成本随产量的变化率（称为边际成本）等等，这些变化率在数学上都可用导数表示． 导数由于其应用的广泛性，为我们解决所学过的有关函数问题提供了一般性的方法，导数是研究函数的切线、单调性、极值与最值等问题的有力工具；运用它可以简捷地解决一些实际问题，导数的概念是用来研究函数在一点及其附近的局部性质的精确工具，而对于函数在某点附近的性质还可以应用另一种方法来研究，就是通过最为简单的线性函数来逼近，这就是微分的方法．微分学是数学分析的重要组成部分，微分中值定理作为微分学的核心，是沟通导数和函数值之间的桥梁， Rolle 中值定理， Lagrange 中值定理， Cauchy 中值定理， Taylor 公式是微分学的基本定理， 统称为微分学的中值定理，这四个定理作为微分学的基本定理，是研究函数形态的有力工具 ] 1[．在微分学中，函数的导数是通过极限定义的，但 当函数用表格给出时，就不可用定义来求其导数，只能用近似方法求数值导数] 2[．最简单 的数值微分公式是用差商近似地代替微商，常见的有 [3] ． ()()() 'f x h f x f x h +-≈ ， ()()() 'f x f x h f x h --≈， ()()() '2f x h f x h f x h +--≈ ． 需要注意的是微分是非常敏感的问题，数据的微小扰动会使结果产生很大的变化] 4[． 目录 第一章数值积分计算的重述 (1) 1.1引言 (1) 1.2问题重述 (2) 第二章复化梯形公式 (3) 2.1 复化梯形公式的算法描述 (3) 2.2 复化梯形公式在C语言中的实现 (3) 2.3 测试结果 (4) 第三章复化simpson公式 (6) 3.1 复化simpson公式的算法描述 (6) 3.2 复化simpson公式在C语言中的实现 (6) 3.3 测试结果 (7) 第四章复化cotes公式 (8) 4.1 复化cotes公式的算法描述 (8) 4.2 复化cotes公式在C语言中的实现 (9) 4.3 测试结果 (10) 第五章Romberg积分法 (11) 5.1 Romberg积分法的算法描述 (11) 5.2 Romberg积分法在C中的实现 (12) 5.3 测试结果 (13) 第六章结果对比分析和体会 (144) 参考文献 (16) 附录 (16) 数值积分?-10 2 dx e x （一） 第一章 数值积分计算的重述 1.1引言 数值积分是积分计算的重要方法，是数值逼近的重要内容，是函数插值的最直接应用，也是工程技术计算中常常遇到的一个问题。在应用上，人们常要求算出具体数值，因此数值积分就成了数值分析的一个重要内容。在更为复杂的计算问题中，数值积分也常常是一个基本组成部分。 在微积分理论中，我们知道了牛顿-莱布尼茨(Newton-Leibniz)公式 ()() () b a f x d x F b F a =-? 其中()F x 是被积函数()f x 的某个原函数。但是随着学习的深入，我们发现一个问题: 对很多实际问题，上述公式却无能为力。这主要是因为:它们或是被积函数没有解析形式的原函数,或是只知道被积函数在一些点上的值,而不知道函数的形式,对此,牛顿—莱布尼茨(Newton-Leibniz)公式就无能为力了。此外,即使被积函数存在原函数,但因找原函数很复杂,人们也不愿花费太多的时间在求原函数上,这些都促使人们寻找定积分近似计算方法的研究,特别是有了计算机后,人们希望这种定积分近似计算方法能在计算机上实现,并保证计算结果的精度,具有这种特性的定积分近似计算方法称为数值积分。由定积分知识,定积分只与被积函数和积分区间有关,而在对被积函数做插值逼近时,多项式的次数越高,对被积函数的光滑程度要求也越高,且会出现Runge 现象。如7n >时,Newton-Cotes 公式就是不稳定的。因而,人们把目标转向积分区间,类似分段插值,把积分区间分割成若干小区间,在每个小区间上使用次数较低的Newton-Cotes 公式,然后把每个小区间上的结果加起来作为函数在整个区间上积分的近似,这就是复化的基本思想。本文主要 基础实验二 定积分数值计算 一、实验目的 学习定积分的数值计算方法，理解定积分的定义，掌握牛顿-莱布尼兹公式。 二、实验材料 2.1定积分的数值计算 计算定积分?b a dx x f )(的近似值，可将积分区间n 等分而得矩形公式 n a b n a b i a f dx x f n i b a ---+≈∑?=]) 1([)(1 或 n a b n a b i a f dx x f n i b a --+≈∑?=][)(1 也可用梯形公式近似计算 n a b b f a f n a b i a f dx x f n i b a -++-+≈∑?-=]2)()()([)(11 如果要准确些，可用辛普森公式 n a b b f a f a b i a f n a b i a f dx x f n i n i b a 6)]()()2)21((4)(2[)(111-++--++-+≈∑∑?=-= 对于?1 0sin xdx ，矩形公式、梯形公式、辛普森公式的Mathematica 程序为 a=0;b=1;k=10; f[x_]:=Sin[x]; d=N[Integrate[f[x],{x,a,b}],k];（计算精确值） s1[m_]:=N[Sum[f[a+i*(b-a)/m]*(b-a)/m,{i,0,m-1}],k];（取小区间左端点的矩形公式） s2[m_]:=N[Sum[f[a+(i+1/2)*(b-a)/m]*(b-a)/m,{i,0,m-1}],k]; （取小区间中点的矩形公式） s3[m_]:=N[Sum[f[a+i*(b-a)/m]*(b-a)/m,{i,1,m}],k]; （取小区间右端点的矩形公式） s4[m_]:=N[Sum[(f[a+i*(b-a)/m]+f[a+(i+1)*(b-a)/m])/2*(b-a)/m,{i,0,m-1}],k]; （梯形公式） s5[m_]:=N[(b-a)/m/6*((f[a]+f[b])+2*Sum[f[a+i*(b-a)/m],{i,1,m-1}] 计算方法实验报告(5) 学生姓名杨贤邦学号指导教师吴明芬实验时间2014.4.16地点综合实验大楼203 实验题目数值积分方法 实验目的●利用复化梯形、辛普森公式和龙贝格数值积分公式计算定积分的 近似植。 实验内容●梯形、辛普森、柯特斯法及其Matlab实现； ●变步长的梯形、辛普森、柯特斯法及其Matlab实现。 ●题目由同学从学习材料中任意选两题 算法分析梯形：function y=jifeng_tixing(a,b,n,fun) fa=feval(fun,a); fb=feval(fun,b); s=0; h=(b-a)/n; for k=1:n-1 xk=a+k*h; s=feval(fun,xk)+s; end y=(h/2)*(fa+fb+2*s); 辛普生：function y=jifeng_xingpu(a,b,n,fun) fa=feval(fun,a); fb=feval(fun,b); h=(b-a)/n; s=0; s2=feval(fun,a+0.5*h); for k=1:n-1 xk=a+k*h; s=feval(fun,xk)+s; s2=feval(fun,xk+(h/2))+s2; end 与源程序y=(h/6)*(fa+fb+2*s+4*s2); 龙贝格：function r2=jifeng_long(fun,a,b,e) h=b-a; t1=(h/2)*(feval(fun,a)+feval(fun,b)); k=1; r1=10; r2=0; c2=0; while abs(r2-r1)>e; s=0; x=a+h/2; while x=3 r1=r2; c2=s2+(1/15)*(s2-s1); r2=c2+(1/63)*(c2-c1); k=k+1;h=h/2; t1=t2;s1=s2; c1=c2; end end 第七章 数值积分 如果函数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 )，然后通过求?b a n dx x P )(得?b a dx x f )(的近似值。 7.1 插值型求积公式 设?=b a dx x f I )(*，插值型求积公式就是构造插值多项式P n (x )，使 ?=≈b a n dx x P I I )(* 。 构造以a ，b 为结点的线性插值多项式)()()(1b f a b a x 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 b a b a +-=?? ? ???--+--==??称为梯形公式。 6 数值积分 如果函数)(x f 在区间],[b a 上连续，且原函数为)(x F ，则可用牛 顿―莱布尼兹公式：)()()(a F b F dx x f b a -=?计算定积分。然而很多函数 无法用牛顿―莱布尼兹公式求定积分。 一个简单被积函数，例如错误!未找到引用源。dx cx bx a ?++2，其不定积分可能很复杂，见下面的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 所以有必要研究简单、高效的计算定积分的方法（即数值积分方法）。数值积分的基本思想是构造一个简单函数)(x P n 来近似代替被积分函数)(x f ，然后通过求?b a n dx x P )(得?b a dx x f )(的近似值。 6.1 插值型求积公式 设?=b a dx x f I )(* ，插值型求积公式就是构造插值多项式)(x P n ，使 ?=≈b a n dx x P I I )(*。 构造以a ，b 为结点的线性插值多项式)()()(1b f a b a x 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 b a b a +-=?? ? ???--+--==??称为梯形公式。 第七章 微积分的数值计算方法 7.1 微积分计算存在的问题/数值积分的基本概念 1. 微分计算问题 求函数的导数(微分)，原则上没有问题。当然，这是指所求函数为连续形式且导数存在的情形。但如果函数一表格形式给出，要求函数在某点的导数值；或者是希望某点的导数值只用其附近离散点上的函数值近似地表示，这就是新问题了，它称为微分的数值计算，或称为数值微分。 2．定积分计算问题 计算函数f 在],[b a 上的定积分 dx x f I b a ?= )( 当被积函数f 的原函数能用有限形式)(x F 给出时，可用积分基本公式来计算： )()()(a F b F dx x f I b a -==? 然而，问题在于：① f 的原函数或者很难找到，或者根本不存在；②f 可能给出一个函数表；③仅仅知道f 是某个无穷级数的和或某个微分方程的解等等。这就迫使人们不得不寻求定积分的近似计算，也称数值积分。 3．数值积分的基本形式 数值积分的基本做法是构造形式如下的近似公式 ∑?=≈n k k k b a x f A dx x f 0 )()( (7.1.1) 或记成 ∑?=+=n k n k k b a f R x f A dx x f 0 ][)()( (7.1.2) ∑==n k k k x f A I 0 * )( 和 ][f R n 分别成为],[b a 上的f 的数值求积公式及其 余项(截断误差)，k x 和k A ),,1,0(n k =分别称为求积节点和求积系数(求积系数与被积函数无关)。 这种求积公式的特点是把求积过(极限过程)程转化为乘法与加法的代数运算。构造这种求积公式需要做的工作是：确定节点k x 及系数 k A ),,1,0(n k =，估计余项][f R n 以及讨论* I 的算法设计及其数值稳定 性。 4．插值型求积公式 如何构造求积公式呢？基本的技术是用被积函数f 的Lagrange 插值多项式 )(x L n 近似代替f ，也即对],[b a 上指定的1+n 个节点 第四章 数值积分 一．问题提出： （1）针对定积分()b a I f x dx =? ，若()5 f x x =,a=0,b=1，即有1 61 500166 x I x dx == =?，但当()sin x f x x = ，()2sin f x x =，……，时，很难找到其原函数。 （2）被积函数并没有具体的解析形式，即()f x 仅为一数表。 二．定积分的几何意义 定积分()b a I f x dx =?的几何意义为，在平面坐标系中I 的值即为四条曲线所围图形的面 积，这四条曲线分别是()y f x =，y=0，x=a ，x=b 。 x y 三．机械求积公式 1.中矩形公式 ()()2b a a b I f x dx b a f +?? =≈- ??? ?； 几何意义：用以下矩形面积替代曲边梯形面积。 x y 2 2.梯形公式 ()()()2b a b a I f x dx f a f b -=≈ +??? ?? 梯形公式的几何意义：用以下梯形面积替代曲边梯形的面积： x y 3.辛普生公式 ()()()462b a b a a b I f x dx f a f f b -? +??? =≈ ++ ? ?????? ? 辛普生公式的几何意义：阴影部分的面积为抛物线曲边梯形，该抛物线由 ()(),(),,,,()22a b a b a f a f b f b ?++? ?? ? ????? 三点构成。 x y a+b 2 4.求积公式的一般形式 ()()0 n b k k a k f x dx A f x =≈∑?，其中k x 称为节点，k A 称为求积系数，或权。 5.求积公式的代数精度（衡量求积公式准确度的一种方法） 含义：衡量一个积分公式的好坏，要用具体的函数来衡量，寻找怎样的函数来衡量呢？简单的多项式函数是一个理想的标准。 定义：若某积分公式对于()0,1,,k x k m = 均能准确成立，但对于1m x +不能准确成立。则称该公式具有m 次代数精度。 解释：代数精度只是衡量积分公式好坏的1种标准。 例1．研究中矩形公式()()2b a a b f x dx b a f +?? ≈- ??? ?的代数精度及几何意义。 解：当()01f x x ==时，公式左边()1b b a a f x dx dx b a ===-??，公式右边b a =-，左=右； 当()1 f x x =时，公式左边()2222 2 b b b a a a x b a f x dx x dx -=== =?? ， 公式右边()22 22a b b a b a +-??=-= ??? ，左=右； 当()2f x x =时，公式左边()3332 33 b b b a a a x b a f x dx x dx -====?? ， 公式右边()2 2a b b a +?? =- ??? ，左≠右； 标题：数值计算及其应用随着计算机的迅速发展和广泛应用，在众多领域内，人们越来使越认识到科学计算是科学研究的第三种方法，数值计算是研究数学问题的数值解及其理论的一个数学分支，它涉及面很广，如：代数、微积分、微分方程、无穷级数、概率论等多方面数学基础知识。自计算机成为数值计算的主要工具来，人们主要研究适合于在计算机上用的数值计算方法及与此相关的理论，包括方法的敛散性、稳定性及误差分析，还要根据计算机的特点研究计算时间最短、需要内存最少的计算方法。它除了具有数学的抽象性与严格性外，还具有应用的广泛性与实际实验的技术性。 数值计算有很多重要的应用，下面举例说明： 1. 在科学技术工程和实验中，经常需要从实验数据中寻找拟合直线，如：天文学家通过对天体运行的观测数据进行分析和处理得到天体的运动轨迹，这就需要用到“多项式逼近”理论和“曲线拟合”的相关知识。 2. 现实生活中经常遇到最优化问题，如：商家寻求最大收益、投资者寻求最小风险等。这就需要用到“数值优化”的知识。 3. 很多数学物理问题都涉及到偏（常）微分方程、科学工程领域建立的许多数学模型也经常用到微分方程，但通常我们无法计算其解析解（事实上也没有必要计算解析解），那么此时数值近似解就具有重要的意义，要求得其数值解就要用到“微分方程求解”的相关理论。 4. 在很多关键领域：如航天领域要研究系统的稳定性，实际上就是研究“收敛”和“发散”，对与这些问题就要用到“方程根的求解”的相关知识。 5. 现实中还有很大一类问题需要求解线性方程组，这就需要“线性方程组求解”及“特征值与特征向量”理论。 综上所述：数值计算在现实生活中发挥着重要的作用，在高科技领域占中有举足轻重的地位！ 若干数值积分的计算方法 黄海琼 （广西民族大学数计学院04数本1班 南宁 530006） 摘 要: 本文讨论了若干数值积分的计算方法。在一维情形下，介绍了Newton-Cotes 公式，Gauss 型等求积法则； 在二维情形下, 主要介绍了二元Newton-Cotes 积分方法。最后，对几类数值积分方法及其数值实验进行比较评述。 关键词: 牛顿-柯特斯公式；Gauss 型求积法则；二元数值积分；数值实验 Some Computational Methods of numerical integration Huang Haiqiong (College of Mathematics and Computer Science,Guangxi University for Nationalities, Nanning 530006) Abstract: In this paper, some computational methods about numerical integration are discussed. under the univariate situation, the quadrature rule of Newton-Cotes formula, Gauss formula and so on is introduced. Under the two-dimensional situation, it mainly introduced the dual Newton-Cotes integral method. Finally, the numerical integration methods and numerical experiment were discussed. Key word: Newton-Cotes formula; Gauss integration principle; dual numerical integration; numerical experiment. 1 引 言 数值积分是积分计算的重要方法，是数值逼近的重要内容，是函数插值的最直接应用，也是工程技 术计算中常常遇到的一个问题[1]。在应用上，人们常要求算出具体数值，因此数值积分就成了数值分析的一个重要内容。在更为复杂的计算问题中，数值积分也常常是一个基本组成部分。 在微积分理论中，我们知道了牛顿-莱布尼茨(Newton-Leibniz )公式 ()()()b a f x dx F b F a =-? 其中()F x 是被积函数()f x 的某个原函数。但是随着学习的深入，我们发现一个问题: 对很多实际问题，上述公式却无能为力。这主要是因为:它们或是被积函数没有解析形式的原函数,或是只知道被积函数在一些点上的值,而不知道函数的形式,对此,牛顿—莱布尼茨(Newton-Leibniz)公式就无能为力了。此外,即使被积函数存在原函数,但因找原函数很复杂,人们也不愿花费太多的时间在求原函数上,这些都促使人们寻找定积分近似计算方法的研究,特别是有了计算机后,人们希望这种定积分近似计算方法能在计算机上实现,并保证计算结果的精度,具有这种特性的定积分近似计算方法称为数值积分。 山西大学计算机与信息技术学院 实验报告 姓 名 学 号 专业班级 课程名称 计算方法 实验日期 成 绩 指导老师 批改日期 实验名称 实验六 数值积分 一． 实验目的： 利用复化梯形公式、复化辛普生公式和龙贝格数值积分公式计算?b a dx x f )(的近似 值。 二． 实验方法： （1） 将[a ，b]区间n 等分，记分点为),1,0,(,n i n a b h ih a x i =-= +=，并在每个小区间[1,+i i x x ]上应用梯形公式 )]()(2)([21 1 b f x f a f h T n i i n ++=∑-= （2） 在每个小区间[1,+i i x x ]上，用辛普生公式 )]()(2)(4)([61 1102 1b f x f x f a f h S n i i n i i n +++=∑∑-=-=+ 式中2 1+ i x 为[1,+i i x x ]的中点，即h x x i i 2 12 1 +=+ （3） 先用梯形公式计算)]()([2/)(1b f a f a b T +?-=，然后，将求积区间(a,b)逐次折半 的方法，令区间长度).,,1,0(,2)(n i a b h i =?-= 计算∑=-?++?=n k n n k h a f h T T 1 2))21 ((2/2/1，式中i n 2=。 于是，得到辛普生公式)/3T -(n 22n n n T T S +=。 柯斯特求积公式15/)(22n n n n S S S C -+=。 最后，得龙贝格求积公式63/)(22n n n n C C C R -+=。 利用上述各公式计算，直到相邻两次的积分结果之差满足精度要求。 三． 实验内容 利用复化梯形公式、复化辛普生公式和龙贝格数值积分公式计算 ? ?+==1 2 21 2 14 dx x dx xe e x π和 的近似值，要求误差为7102 1 -?=ε，将计算结果与精确值比较，并对计算结果进行 分析（计算量、误差） 四． 实验程序： 复合梯形公式： #include 微积分的数值计算方法 7.1 微积分计算存在的问题/数值积分的基本概念 微分计算问题 求函数的导数(微分)，原则上没有问题。当然，这是指所求函数为连续形式且导数存在的情形。但如果函数一表格形式给出，要求函数在某点的导数值；或者是希望某点的导数值只用其附近离散点上的函数值近似地表示，这就是新问题了，它称为微分的数值计算，或称为数值微分。 2．定积分计算问题 计算函数在上的定积分 当被积函数的原函数能用有限形式给出时，可用积分基本公式来计算： 然而，问题在于：① 的原函数或者很难找到，或者根本不存在；②可能给出一个函数表；③仅仅知道是某个无穷级数的和或某个微分方程的解等等。这就迫使人们不得不寻求定积分的近似计算，也称数值积分。 3．数值积分的基本形式 数值积分的基本做法是构造形式如下的近似公式 (7.1.1) 或记成 (7.1.2) f ],[b a dx x f I b a ?=)(f )(x F )()()(a F b F dx x f I b a -==?f f f ∑?=≈n k k k b a x f A dx x f 0 )()(∑?=+=n k n k k b a f R x f A dx x f 0 ][)()( 和 分别成为上的的数值求积 公式及其余项(截断误差)，和分别称为求积节点和求积系数(求积系数与被积函数无关)。 这种求积公式的特点是把求积过(极限过程)程转化为乘法与加法的代数运算。构造这种求积公式需要做的工作是：确定节点及系数，估计余项以及讨论 的算法设计及其数值稳定性。 4．插值型求积公式 如何构造求积公式呢？基本的技术是用被积函数的Lagrange 插值多项式近似代替，也即对上指定的个节点 及相应的函数值 ，作 代入(7.1.2)式等号左边有 或写成形如(7.1.2)式的一般形式： ∑==n k k k x f A I 0 * )(][f R n ],[b a f k x k A ),,1,0(n k Λ=k x k A ),,1,0(n k Λ=][f R n *I f )(x L n f ],[b a 1+n b x n ≤ 数值积分 第1章 理论依据 逼近论——构造一个简单函数p(x)近似表示f(x),然后对 p(x)求积分得到 f(x)的积分的近似值。基于插值原理，推导出数值积分的基本公式。 §1插值求积公式 为了用数值方法求 b a I(f)=f(x)dx ? ，对被积函数f(x)在给定的n+1个节点 上作Lagrange 插值，用插值函数Pn(x)代替f(x)，就可用I （Pn(x)）构造求积公式，近似地计算定积分I(f(x))。 §2Newton —Cotes 公式 §2.1Newton —Cotes 公式的推导 当§1.1插值求积公式的插值节点为等距节点时，就得到Newton —Cotes 公式。 将区间[a,b]n 等分， b a h n -= ，n+1个节点为 x k =a+kh (k=0,1,…,n) 在节点上对f(x)的Lagrange 插值多项式是： 0()()() n n j n k k j k j j k x x p x f x x x ==≠-=-∑∏ 用P n (x)代替f(x)构造求积公式： 0()()()n n b b j n n k a a k j k j j k x x I p x dx f x dx x x ==≠-==-∑∏?? 记，(k=0,1,…,n) 作代换x=a+th 带入上式，变为： () 00()n n n n k k j j k b a t j A dt b a C n k j =≠? --==--∏? 其中： (k=0,1,…,n) （1-1） 这个积分是有理多项式积分，它与被积函数f(x)和区间[a,b]无关。只要确定n 就能计算出系数。 于是得到称为Newton —Cotes 公式的求积公式： ()0 ()n n n k k k I b a C y ==-∑ （1-2） 其中称为Newton —Cotes 系数。如表1所示。 §2.2Newton —Cotes 公式误差和稳定性 在积分公式中用插值多项式Pn(x)代替f(x)的插值误差是 (1)0 ()()()()()(1)!n n n n k k f R x f x p x x x n ξ+==-=-+∏ 因此，Newton —Cotes 公式的截断误差是 (1)0 ()()()(1)!n n b k a k f R f x x dx n ξ+==-+∏? （1-3） 讨论舍入误差对计算结果产生的影响，设（1-2）式近似计算()b a f x dx ? 其中计算函数值f(xn)有误差值（k=0,1,2, …,n ）。在（1-2）式中令 ? 设计算无误差，舍入误差也忽略，则，由（1-2）式 计算时引式的误差为 () ()()() 0000()[()(())()(...) n n n n n n n k k k k n n n k k e b a C f x C f x b a C C εεε===--+=--++∑∑ 如果皆为正，并设，则 ，故 有 学科分类号110.3420 本科毕业论文 题目几种常用数值积分方法的比较 姓名潘晓祥学号1006020540200 院（系）数学与计算机科学学院 专业数学与应用数学年级2010 级 指导教师雍进军职称讲师 二〇一四年五月 贵州师范学院本科毕业论文（设计）诚信声明 本人郑重声明：所呈交的本科毕业论文（设计），是本人在指导老师的指导下，独立进行研究工作所取得的成果，成果不存在知识产权争议，除文中已经注明引用的内容外，本论文不含任何其他个人或集体已经发表或撰写过的作品成果。对本文的研究做出重要贡献的个人和集体均已在文中以明确方式标明。本人完全意识到本声明的法律结果由本人承担。 本科毕业论文作者签名： 年月日 贵州师范学院本科毕业论文(设计)任务书 毕业设计题目几种常用数值积分方法的比较 作者姓名潘晓祥学号1006020540200 年级2010级 所属学院数学与计算机科 学学院 专业数学与应用数学班级四班 指导教师签名雍进军讲师职称讲师开题日 期 2013年7月10日 主要目标 1．了解什么数值积分基本思想和一些常用的数值积分方法； 2．对各种数值积分方法的误差以及代数精度进行分析； 3．对各积分方法进行比较总结出优缺点。 主要要求 通过对几种常用的数值积分方法进行了的分析，并用这几种方法对被积函数是普通函数做了数值积分，并在计算机上进行实验。数值积分是计算方法或数值分析理论中非常重要的内容，数值积分方法也是解决实际计算问题的重要方法，对几种常用数值积分方法的分析很必要。 主要内容 本文通过对复化求积公式, Newton—Cotes求积公式, Romberg求积公式,高斯型求积公式进行分析讨论并在计算机上积分实验,从代数精度,求积公式误差等角度对这些方法进行分析比较，并总结出每种求积分法的优缺点以及实用性。 摘要 本文应用插值积分法和逼近论的思想，简单重述了推导Newton-Cotes公式和Gauss-Legendre求积公式的过程，以及这两个公式的系数、精度等问题。并以这两种数值积分的求解方法为基础，应用quad、guass函数编写具体Matlab程序，通过计算机软件计算出所给题目的近似数值积分。对二者所得的结果进行比较，从而研究了用Newton-Cotes和Gauss-Legendre公式求积分的方法和二者的精确度问题。得知，这两种求积公式所得的结果在精度上的确存在差异，结合理论部分更加充分地说明了，n相同时Gauss-Legendre公式比Newton-Cotes公式具有更高的代数精度，但当代数精度相同时，二者计算的结果仍存在细微的差异。 关键字：插值积分、Newton-Cotes公式、Gauss-Legendre公式 数值积分 1 理论依据 逼近论——构造一个简单函数p(x)近似表示f(x),然后对 p(x)求积分得到 f(x)的积分的近似值。基于插值原理，推导出数值积分的基本公式。 1.1插值求积公式 为了用数值方法求 b a I(f)=f(x)dx ，对被积函数f(x)在给定的n+1个节点上 作Lagrange插值，用插值函数Pn(x)代替f(x)，就可用I（Pn(x)）构造求积公式，近似地计算定积分I(f(x))。 2 Newton —Cotes 公式 2.1 Newton —Cotes 公式的推导 当§1.1插值求积公式的插值节点为等距节点时，就得到Newton —Cotes 公式。 将区间[a,b]n 等分， b a h n -= ，n+1个节点为 x k =a+kh (k=0,1,…,n) 在节点上对f(x)的Lagrange 插值多项式是： 0()()() n n j n k k j k j j k x x p x f x x x ==≠-=-∑∏ 用P n (x)代替f(x)构造求积公式： 0()()()n n b b j n n k a a k j k j j k x x I p x dx f x dx x x ==≠-==-∑∏?? 记错误！未找到引用源。，错误！未找到引用源。(k=0,1,…,n) 作代换x=a+th 带入上式，变为： () 00()n n n n k k j j k b a t j A dt b a C n k j =≠? --= =--∏? 其中：错误！未找到引用源。 (k=0,1,…,n) （1-1） 这个积分是有理多项式积分，它与被积函数f(x)和区间[a,b]无关。只要确定n 就能计算出系数错误！未找到引用源。。 于是得到称为Newton —Cotes 公式的求积公式： ()0()n n n k k k I b a C y ==-∑ （1-2） 其中错误！未找到引用源。称为Newton —Cotes 系数。如表1所示。 表1 Newton —Cotes 系数 n 1 1/ 2 1/2 2 1/6 4/6 1/6几种定积分的数值计算方法
(精选)实验二 数值方法计算积分
计算方法-数值积分实验
导数的数值计算方法[文献综述]
数值积分 (论文)
基础实验二 定积分数值计算
matlab计算方法实验报告5(数值积分)
计算方法讲义：七 数值积分
工程中的计算方法课件6 数值积分
微积分的数值计算方法
数值计算方法教案_数值积分
数值计算方法课程总结
若干数值积分的计算方法
计算方法实验六 数值积分
微积分的数值计算方法
数值积分-计算方法
几种常用数值积分方法的比较
数值积分的计算方法论文
相关主题
文本预览