几种定积分的数值计算方法摘要:本文归纳了定积分近似计算中的几种常用方法,并着重分析了各种数值方法的计算思想,结合实例,对其优劣性作了简要说明.关键词:数值方法;矩形法;梯形法;抛物线法;类矩形;类梯形Several Numerical Methods for Solving Definite IntegralsAbstract: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 trapezoid1. 引言在科学研究和实际生产中,经常遇到求积分的计算问题,由积分学知识可知,若函数)(x f 在区间],[b a 连续且原函数为)(x F ,则可用牛顿-莱布尼茨公式⎰-=baa Fb F x f )()()(求得积分.这个公式不论在理论上还是在解决实际问题中都起到了很大的作用. 在科学研究和实际生产中,经常遇到求积分的计算问题,由积分学知识可知,若函数)(x f 在区间],[b a 连续且原函数为)(x F ,则可用牛顿-莱布尼茨公式 ⎰-=baa Fb F x f )()()(求得积分.这个公式不论在理论上还是在解决实际问题中都起到了很大的作用.另外,对于求导数也有一系列的求导公式和求导法则.但是,在实际问题中遇到求积分的计算,经常会有这样的情况:(1)函数)(x f 的原函数无法用初等函数给出.例如积分 dx ex ⎰-102, ⎰10sin dx xx等,从而无法用牛顿-莱布尼茨公式计算出积分。

(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,其中nab h -=表示小区间的长度. 2.1矩形法矩形法就是用小矩形面积近似代替各个小曲边梯形面积,从面积得到S的近似值.若取小区间左端点的函数值为小矩形的高,如图1中所示,则∑=-=niixfnabA1).(图1 分割曲边矩形近似积分2.2 梯形法梯形法则用小直边梯形的面积近似代替小曲边梯形面积,见图2,从而得到S的近似值,即⎥⎦⎤⎢⎣⎡++-=∑-=11)(2)()(niixfbfafnabA.图2 分割曲边梯形近似积分2.3抛物线法抛物线法以抛物线为曲边梯形的曲边,曲边梯形的面积近似代替小曲边梯形的面积,如图3所示.图3 抛物线积分210,,x x x 对应的曲线上的点210,,P P P 可以唯一地确定一条抛物线c bx ax y ++=2,这条抛物线将作将代替从0x 至2x 的曲线段,此时积分可以转化为对抛物线积分,而抛物线的积分可以利用牛顿—莱布尼玆公式.第1、2个小区边梯形的面积:)]()(4)([3)2102120x f x f x f hdx c bx ax A x x ++=++=⎰（上面利用了条件210,,P P P 是抛物线上的点以及等式1022x x x =+.同理可证： )]()(4)([3h4322x f x f x f A ++=……)]()(4)([3122/n n n n x f x f x f hA ++=--所以,})(2)(4)]()({[12/122/11232/21∑∑-==--+++=+++≈n i i n i i na b n x f x f b f a f A A A S3.概率意义上的数值算法概率算法是定积分问题数值求解的一类常用方法,其设计思想简单,易于实现 .尽管算法要耗费较多计算时间,但是往往能得到问题的近似解,并且近似程度能随计算时间的增加而不断提高.概率算法可用于计算定积分的近似值.3.1平均值法考虑定积分⎰=ba dx x f I )(的近似计算,其中)(x f 在[]b a ,内可积,用平均值法计算该积分,首先随机产生n 个独立的随机变量,且服从在[]b a ,上均匀分布,即),2,1(n i i =ξ;其次,计算I 的近似值I ,∑=-=ni i f n a b I 1)(ξ.由中心极限定理知,若{}),2,1(n i i =ξ相互独立、同分布,且数学期望及标准差0>σ存在,则当n 充分大时,随机变量nII Y σ-=渐近服从正态分布)1,0(N ,即对任意的0>αt ,}t I -I P{}t Y P{nσαα<=<这表明,用平均值法计算定积分的收敛速度较慢,在概率意义下的误差阶仅为)1(n O .3.2“类矩形”Monte-Carlo 方法由于平均值法计算定积分的收敛速度较慢,且在概率意义下的误差阶仅为)1(n O ,就有对平均值法的改进，“类矩形” Monte-Carlo 方法,改进过程为：先将积分区间[]b a ,n 等分, 随机产生n 个相互独立且服从[]1,0上均匀分布的随机变量序列),2,1(},{n i i =ξ；然后由这n 个随机点类似于矩形公式构造计算公式,即作变换 n i i nab a i i ,2,1),1(=-+-+=ξη将}{i ξ映射到子区间 []n i b a a b nia ab n i a ,,2,1,,)}(),(1{ =⊂-+--+最后,计算I 的近似值I ~,∑=-=ni i f n a b I 1)(~η． 下面用两个命题证明“类矩阵”方法的可行性. 命题1 设[][][]有记,,,)(max ,,)(0,1b a x x f M b a C x f b a x ∈∀'=∈∈M a b a b x f dx x f ba2)())(()(20-≤--⎰证明:由Lagrange 中值定理得)())(()()(000之间与介于x x x x f x f x f ξξ-'+=上式两边在[]b a ,积分,得⎰⎰-'+-=babadx x x f a b x f dx x f ))(())(()(00ξ由)(x f '得连续性,得.2)()(21)())(())(()(222020000M a b b a x b a x M dxx x M dx x x f a b x f dx x f b ababa-≤⎥⎦⎤⎢⎣⎡++--=-≤-'≤--⎰⎰⎰ξ命题2 设[],,,)(1nab h b a C x f -=∈[][]n i x f M x f M ih a h i a x i b a x ,.2.1,)(max,)(max ,)1(, ='='=+-+∈∈I ~与I 如上,则I ~与I 的误差满足)1(~nO I I =-.证明: ⎰∑=--=-ba ni i f n a b dx x f I I 1)()(~η ∑⎰∑=+-+=-=ni iha hi a ni i f h dx x f 1)1(1)()(η∑⎰=+-+-≤n i iha hi a i hf dx x f 1)1()()(η由命题1得, n i h M hf dx x f i iha hi a i ,,2,1,2)()(2)1( =≤-⎰+-+η 于是∑=-≤≤-ni i a b n M h M I I 122)(22~即)1(~nO I I =-.3.3“类梯形”Monte-Carlo 方法再给出平均值法的另一种改进.首先将[]b a ,n 等分,再在每个子区间上随机产生n 2个相互独立且服从]1,0[上均匀分布的随机变量序列,并两两分组,得),,3,2,1(},,{212n i i i =-ξξ;做变换)12(2)22(2221212i i i i i nab a i nab a ξηξη+--+=+--+=--将12-i ξ,i 2ξ分别映射到子区间ni a b ni a a b n i a a b ni a a b n i a ,,3,2,1)],(),(212[)](212),(1[ =-+--+--+--+然后在每个等分子区间上)](),(1[a b nia ab n i a -+--+利用i i 212,ηη-两点类似于梯形公式构造“类梯形”公式 )]()([212i i i f f nab S ηη+-=- 近类似⎰+-+ih a hi a dx x f )1()(.最后计算I 的近似值I ~~,∑=-+-=n i i i f f n a b I 12122)()(~~ηη． 下面证明“类梯形”方法可行性的两个命题：命题3 设()[]2,f x a b ∈C ，记[](),max ''x a b f x ∈M=，则()1212,x x a x x ∀≤≤，有()()()()312212bab a b af x dx f x f x M ---+≤⎡⎤⎣⎦⎰． 证明： 过()()()()1122,,,x f x x f x 两点的直线方程为()()()211121()()f x f x P x f x x x x x -=+--所以 ()(),1,2.i i P x f x i ==令12()()()()()()R x f x P x k x x x x x =-=-- （1）将x 看成[],a b 上的一个定点,构造辅助函数12()()()()()()t f t P t k x t x t x φ=----由于12()()()0x x x φφφ===,由Rolle 中值定理,'()t φ在(),a b 内至少有两个零点,对'()t φ再用Rolle 中值定理,知''()t φ在(),a b 内至少有一个零点,即存在(),a b ξ∈,使''()''()2()0f k x φξξ=-=,所以''()()2f k x ξ=．将它代入（1）式,并两段同时从a 到b 积分,得()()()12121212121212122''()()()2()()2()()()()()()2ba babax x b a x x b af x dx f x f x f x x x x dx M x x x x dxM x x x x dx x x x x dx x x x x dx ξ--+⎡⎤⎣⎦=--≤--⎡⎤=--+--+--⎢⎥⎢⎥⎣⎦⎰⎰⎰⎰⎰⎰记121212121212(,)()()()()()()x x bax x L x x x x x x dx x x x x dx x x x x dx =--+--+--⎰⎰⎰不妨设12a x x b <<<,则将12(,)L x x 分别对求偏导数,得1222212121()()()02()()()02x x a bL b a x x x a bL b a x x x +=----=+=--+-=解得唯一驻点：121(3)41(3)4x a b x a b ⎧=+⎪⎪⎨⎪=+⎪⎩ 又3333()()(,),(,)44166a b a b b a b a L L a b ++--==故当12a x x b ≤<≤时,()()()()31212(,)2212bab a b a M f x dx f x f x L x x M ---+≤≤⎡⎤⎣⎦⎰结论成立．命题4 [],,)(2b a C x f ∈设[][]n i x f M x f M nab h ih a h i a x i b a x ,2,1,)(max ,)(max ,,)1(,=''=''=-=+-+∈∈I 与I ~~如上,则I 与I ~~ 的误差满足:)1(~~2nO I I =-.证明：∑⎰∑⎰∑⎰∑=+-+-=+-+=-=-+-≤+-=---=-n i iha hi a i i ni ih a h i a ni i i ba n i i i hf f dx x f f f h dx x f f n a b dx x f I I 1)1(2121)1(121212122)()()(2)()()(2)()(~~ηηηηηη由命题3,得n i h Mh f f dx x f i iha hi a i i ,,2,1,122)()()(3)1(212 =≤+-⎰+-+-ηη于是∑=-≤≤-n i i n M a b h M I I 123312)(12~~即)1(~~2nO I I =-.4.例题对于积分dx 14102⎰+x ,该积分精确值为 3.1416.下面分别给出本文所涉及计算方法对它的计算结果:4.1用三种基于几何意义的算法:矩形算法,梯形法,抛物线法作比较,结果如表1:表1 几何意义算法的比较分割数算法 近似值 误差110矩形3.14241294 13.8-e 梯形 3.1399398 37.1-e 抛物线 3.1415569 81.4-e 210矩形3.1415528 13.8-e梯形 3.1416496 57.1-e 抛物线3.14160128110.4-e4.2用平均值法,及其改进“类矩形”Monte-Carlo 方法, “类梯形”Monte-Carlo 方法计算结果如表2:表2 概率意义算法的比较5．结语本文介绍的几种求积公式各有特点：梯形求积公式和抛物线法求积公式是低精度公式,但对于光滑性较差的被积函数有时比用高精度方法能得到更好的效果,尤其是梯形求积公式．当被积函数为周期函数时,效果更为突出.由表1分析,一般情形下,三种基于几何的算法中矩形算法的误差最大,梯形法次之,抛物线法最高.抛物线法的积分精度远远高于另外两种方法,特别是在积分区间分割份数较小的情况下,仍然保持较高的近似程度.“类矩形”Monte-Carlo 方法; “类梯形”Monte-Carlo 方法是平均值法的改进,提高了平均值法的精确度．通过表2可以看出,直接用平均值法计算定积分,410节点的计算已经很可观了,但计算结果只有2位有效数字,而选取同样的节点数,计算量几乎不变,类矩阵法就达到了4位有效数字,类梯形法则达到了8位有效数字,恰好与上述定理中误差阶的估计是一致的,从而也验证了“类矩形”Monte-Carlo 方法和”类梯形”Monte-Carlo 方法的高效性．从表2中也可以看出随着节点数的增大,积分精度会不断提高,当然计算复杂度就会增加.节点数算法 近似值 误差410平均值法3.13849728 31.3--e类矩形法 3.1416903 50.7-e类梯形法3.1416002990.2--e参考文献[1] 费祥历,刘奋,马铭福.高等数学(第2版上册)[M].山东:石油大学出版社,2008: 211-287．[2] 徐萃薇,孙绳武.计算方法引论(第三版).北京:高等教育出版社,2007.[3] 王晓东.计算机算法分析与设计[M].北京:电子工业出版社,2001:197-228.[4] 徐钟济.蒙特卡罗方法[M].上海:上海科学技术出版社,1985.[5] 张平文,李铁军,数值分析[M].北京:北京大学出版社,2007.[6] 阮宗利.计算一元定积分的若干数值算法及其比较[J].中国石油大学学报(科技教育),2010:182-184.[7] 朱长青. 计算方法及其应用.北京:科学出版社,2006.[8] 张威,刘志军,李艳红.数值分析与科学计算.北京:清华大学出版社,200,5[9] 明万元,郑华盛.求解数值积分的两类新的Monte-Carlo方法[J].南昌航空大学学报,2010,40(10):180-186.[10]刘长虹,关永亮等.蒙特卡洛在数值积分上的应用[J].上海工程技术大学学报,2010,24(1):43-46.[11]王岩.Monte-Carlo方法应用研究[J].云南大学学报(自然科学版),2006,28(SI):23-26.10。