数值计算实习报告——牛顿—柯特斯积分方法及其应用姓名：杨银月学号：139084154班级：数131牛顿—柯特斯积分算法及其应用一、引言●数值积分的必要性现实生活当中往往会遇到这样的问题：拉着一块物体在一粗糙平面上沿着一直线从一点移动到另一点所发生的功是多少？一个不规则的平面图形的面积是多少？已知边际成本—产量函数求在一产量下的总成本……通过物理和几何学以及经济学的知识容易知道这类问题是需要计算积分的，根据微积分定理对于积分⎰=badx x f I )(，只要找到)(x f 的原函数)(x F ，便有牛顿—莱布尼茨公式：)()()(a F b F dx x f ba-=⎰但是，现实生活中往往只能得到一些离散的点，无法得到连续的函数)(x f ，即便是给定了)(x f ，也不一定就是容易找到原函数的（原函数往往非初等），比方说)0(sin ≠x xx，2x e-等，故不能用N-L 公式进行积分运算。

即使能求得原函数的积分有时候计算也是非常困难的。

例如对于被积函数611)(xx f +=，其原函数Cx x x x x x x x F ++-+++-+=1313ln 341)1arctan(61arctan 31)(22计算)(),(b F a F 仍然很困难，因此有必要研究积分的数值计算问题。

●数值积分的基本思想由积分中值定理（如图1）知，在积分区间],[b a 内存在一点ζ，成立)()()(ζf a b dx x f ba-=⎰图1就是说，底为a b -而高为)(ζf 的矩形面积恰等于所求曲边梯形的面积I 。

问题在于点ζ的具体位置一般是不知道的，因而难以准确算出)(ζf 的值，称)(ζf 为区间],[b a 上的平均高度。

这样只要对平均高度)(ζf 提供一种算法，相应地便获得一种数值求积方法。

如果我们用两端点“高度”)(a f 与)(b f 的算术平均值作为平均高度)(ζf 的近似值，这样导出的求积公式)]()([2)(b f a f ab dx x f ba+-≈⎰)11(-便是众所周知的梯形公式。

而如果改用区间的中点2ba c +=的“高度”)(c f 近似地取代平均高度)(ζf ，则又可导出所谓的“中矩形公式”。

更一般地，可在区间],[b a 上适当选取某些节点k x ，然后用)(k x f 的加权平均得到平均高度)(ζf 的近似值，这样构造出的求积公式具有下列形式：∑⎰=≈nk k k bax f A dx x f 0)()(，)21(-式中k x 称为求积节点；k A 称为求积系数，亦称伴随节点k x 的权。

权k A 仅仅与节点k x 的选取有关，而不依赖于被积函数)(x f 的具体形式。

这类数值积分通常称为机械求积，其特点是将积分求值问题归结为被积函数值的计算，这就避开了N-L 公式需要寻求原函数的困难。

梯形公式即机械型求积公式1=n的情况，而由公式)21(-，给定不同的n 可以得到不同的求积公式。

下面讨论不同n 的情况下，求积公式的形式与性质。

二、N-C 公式及其简单形式（2,1==n n ）●N-C 公式推导对于大部分的函数)(x f 都是不容易求积的，但是可以用插值法的思想将)(x f 用多项式插值表示，因此可以定义插值型的求积公式：∑⎰=≈nk k k bax f A dx x f 0)()(其中⎰=bak k dx x l A )(。

等距节点的插值型求积公式称为牛顿—柯特斯公式，推导如下：取等距节点：ih a x i +=，nab h -=,n i ,...,2,1=，令th a x +=，得：⎰∏⎰∏⎰∏⎰≠-≠≠----=--=--==nji in ni j ba i j j i jba i i dt j t i n ni ab hdt ji j t dx x x x x dx x l A 00)()!(!)1)(()(，称⎰∏≠----=-n ji i n i dt j t i n ni a b A 0)()!(!)1(为柯特斯系数，记做)(n i C ，所以有牛顿—柯特斯公式：∑⎰=-≈ni i n i bax f C a b dx x f 0)()()()(由于是多项式积分，柯特斯系数的计算不会遇到困难。

当1=n 时，21)1(1)1(0==C C ,这时的求积公式就是梯形公式；当2=n ，这时的柯特斯系数为:61)2)(1(4120)2(0=--=⎰dt t t C ；64)2(2120)2(1=--=⎰dt t t C ；61)1(4120)2(2=-=⎰dt t t C .相应的求积公式是如下辛普森公式)].(2(4)([6b f b a f a f a b S +++-=)12(-●稳定性分析下表是柯特斯公式的系数表（n=1,2，…，8）（程序2—1）n)(n i C 10.50000.500020.16670.66670.166730.12500.37500.37500.125040.07780.35560.13330.35560.077850.06600.26040.17360.17360.26040.066060.04880.25710.03210.32380.03210.25710.048870.04350.20700.07660.17300.17300.07660.20700.043580.03490.2077-0.03270.3702-0.16010.3702-0.03270.20770.0349从表中可以看出当8≥n 时，柯特斯系数)(n i C 出现负值，于是有10)(0)(=≥∑∑==ni n i ni n iC C,如果0))((~)(>-i i n i f x f C，且δ=-~)(i i f x f ，则有∑∑==-=-=-ni i i n ini i i n in n f x f C f x f Cf I f I 0~)(0~)(~])([])([)()(δδ>=-=∑∑==ni n i ni i i n iC f x f C)(0~)()(.这表明初始数据误差将会引起计算结果误差增大，即计算不稳定，故8≥n 时的牛顿—柯特斯公式是不用的。

牛顿—柯特斯公式通常只用4,2,1=n 时的公式，下面只讨论2,1=n 时的误差。

●梯形公式与辛普森公式的误差分析代数精度：如果某个求积公式对于次数不超过m 的多项式均能准确地成立，但对于1+m 次多项式就不准确成立，则称该求积公式具有m 次代数精度。

定理：若求积公式的代数精度为m ，则求积公式的余项可表示为：)()()(][)1(0η+==-=∑⎰m ni i i baKf x f A dx x f f R ,)22(-其中K 为不依赖于)(x f 的待定参数，),(b a ∈η。

不难得到梯形公式的代数精度为1，辛普森公式的代数精度为3，所以梯形公式的误差余项可表示为：),(),(][b a f K f R ∈''=ηη，其中：332233)(121])(61[21)](2)(31[21a b a b b a a b a b K --=--=+---=,所以梯形公式的余项为),(),(12)(][3b a f a b f R ∈''--=ηη，)32(-同理可得，辛普森公式的误差余项为()),(,)2(180][)4(4b a f a b a b f R ∈---=ηη，)42(-●复合求积公式（梯形）（程序2—2）由上述分析知当8≥n 时N-C 公式不具有稳定性，因此不能通过提高求积公式的阶数来提高求积精度，要换一种思路来解决提高精度的问题。

在前面的分析中，求积运算都是在一整个给定的区间上进行的，既然能够在一个大区间上面进行求积计算，那么把区间分成若干个小的区间时也能够在每一个区间上进行求积计算的。

不妨将区间],[b a 分成n 等分，若我们令分点为ih a x i +=，步长为nab h -=，n i ,...,2,1,0=，在每一个子区间)1,...,1,0](,[1-=+n i x x i i 上采用梯形公式，则得)()]()([2)()(1111f R x f x f h dx x f dx x f I n n i i i n i x x bai i++===∑∑⎰⎰-=+-=+,记11101[()()][()2()()],22n n n i i i i i h hT f x f x f a f x f b --+===+=++∑∑)52(-称之为复合梯形公式，其余项为：∑-=+∈''-=-=113),()],(12[][n i i i i i n n x x f h T I f R ηη。

由于],[)(2b a C x f ∈，且)(max )(1)(min 11010i n i n i i i n i f f n f ηηη''≤''≤''∑-=-≤≤-≤≤，所以),(b a ∈∃η使得∑-=''=''1)(1)(n i i f n f ηη，于是复合梯形公式的余项为：)(12][2ηf h a b f R n ''--=，)62(-可以看出误差是2h 阶，且当],[)(2b a C x f ∈，有0][lim =∞→f R n n ，所以有⎰=∞→ban n dx x f T )(lim ，因此复合梯形公式是收敛的。

事实上，只要设],[)(b a C x f ∈，则可得到收敛性，因为只要把n T 改写为])()([21110∑∑=-=-+-=ni i n i i n x f n a b x f n a b T当∞→n 时，上式的右端括号内两个和式均收敛到积分dx x f ba⎰)(，所以复合梯形公式收敛。

此外，n T 的求积系数为正，故复合梯形公式是稳定的。

●复合求积公式（辛普森）（程序2—3）将区间],[b a 分为n 等分，在每个子区间],[1+i i x x 上采用辛普森公式，若记,22/1hx x i i +=+则得：][)]()(4)([6)()(112/111f R x f x f x f h dx x f dx x f I n n i i i i n i x x bai i+++===∑∑⎰⎰-=++-=+,记],)()(2)(4)([6)]()(4)([611102/11012/1∑∑∑-=-=+-=+++++=++=n i i n i i n i i i i n b f x f x f a f hx f x f x f h S )72(-称为复合辛普森求积公式，其余项为：),(,)()2(180][110)4(4+-=∈-=-=∑i i i n i i n n x x f h h S I f R ηη，于是当],[)(4b a C x f ∈时，与复合梯形公式相似有),(),(2(180][)4(4b a f h a b S I f R n n ∈--=-=ηη，)82(-可以看出，误差阶为4h ，收敛性是显然的，实际上，只要],[)(b a C x f ∈，则可以得到收敛性，即.)(lim ⎰=∞→ban n dx x f S 此外，由于n S 中的求积系数均为正数，故知复合辛普森公式计算稳定。