第五章数值积分和微分

1、第五章第五章 数值积分和数值微分数值积分和数值微分n 定积分问题的数值解法定积分问题的数值解法主要内容主要内容n导数或微分数值计算导数或微分数值计算 传统方法的困境传统方法的困境 数值积分的基本思想数值积分的基本思想 数值积分的一般形式数值积分的一般形式 代数精度问题代数精度问题求函数求函数 f f( (x x) ) 在区间在区间 a a, ,b b 上的定积分上的定积分 ( )baIf x dx是微积分学中的基本问题。是微积分学中的基本问题。 5.1 数值积分概述数值积分概述badxxffI)()(对于积分对于积分公式有公式有则由则由的原函数的原函数如果知道如果知道LeibnizNewton

2、xFxf ),()( badxxf)()()()(aFbFxFba 但是在工程技术和科学研究中但是在工程技术和科学研究中, ,常会见到以下现象常会见到以下现象: :的的一一些些数数值值只只给给出出了了的的解解析析式式根根本本不不存存在在)(,)()1(xfxf传统方法的困境传统方法的困境2120sind ,d ,sindbbxaaxexxxxx（2）原函数）原函数F（x）表达不出来表达不出来， F（x）不是初等函数不是初等函数以上这些现象以上这些现象,Newton-Leibniz,Newton-Leibniz很难发挥作用很难发挥作用! !只能建立积分的近似计算方法只能建立积分的近似计算方法-数

3、值积分数值积分正是为解决这样的困难而提出来的，正是为解决这样的困难而提出来的，不仅如此，数值积分也是微分方程数值解法的不仅如此，数值积分也是微分方程数值解法的工具之一。工具之一。求求原原函函数数较较困困难难的的表表达达式式结结构构复复杂杂 ,)()3(xf数值积分的基本思想数值积分的基本思想 数值积分数值积分-指计算定积分近似值的各种计算方法。指计算定积分近似值的各种计算方法。常用一个简单函数代替原来的复杂函数求积分。常用一个简单函数代替原来的复杂函数求积分。 从几何上看，就是计算曲从几何上看，就是计算曲边梯形面积的近似值。边梯形面积的近似值。 最简单的办法，是用最简单的办法，是用直线直线、抛

4、物线抛物线等代替曲边，使得等代替曲边，使得面积容易计算。面积容易计算。( )baIf x dxf(x)abf(a)f(x)abf(x)abf(b)f(a)(a+b)/2左矩公式左矩公式 abafdxxfba )(中矩公式中矩公式 abbafdxxfba 2)(梯形公式梯形公式 bfafabdxxfba 2)(用用直线直线代替曲边代替曲边抛物线公式抛物线公式( ) ( )4 ()( )62babaabf x dxf aff b 用用抛物线抛物线代替曲边代替曲边又称辛卜森公式又称辛卜森公式数值积分的一般形式数值积分的一般形式 正是由于权系数的构造方法不同，从而决定了数值积正是由于权系数的构造方法不

5、同，从而决定了数值积分的不同方法。分的不同方法。 上述的近似求积公式都是取上述的近似求积公式都是取a,ba,b上若干点上若干点 处的高度通处的高度通过加权过加权 后再进行求和得到积分的近似值后再进行求和得到积分的近似值，写成一般形式：，写成一般形式：kAix 0nbkkakfx dxA fx 0nbkkakfx dxA fxR f 或写为：或写为：其中，其中， - 称为称为求积节点求积节点 A Ak k-称为节点称为节点 x xk k 上的上的权系数权系数。 -是函数是函数f f( (x x) )在节点在节点 x xk k 上的函数值，上的函数值， -称为求积公式的称为求积公式的截断误差或余项

6、截断误差或余项。kx ixf fR利用插值多项式来构造数值求积公式利用插值多项式来构造数值求积公式, ,具体步骤如下具体步骤如下: :上上取取一一组组节节点点在在积积分分区区间间,babxxxan 10次次插插值值多多项项式式的的作作nxf)( nkkkxlxfxP0)()()(为插值基函数为插值基函数其中：其中：), 1 , 0)(nkxlk 插值型求积公式插值型求积公式思思路路利用利用插值多项式插值多项式 则积分易算。则积分易算。)()(xfxPn 以拉格朗日插值多项式为例以拉格朗日插值多项式为例有有的近似的近似作为被积函数作为被积函数用用,)()(xfxP badxxf)( badxxP

7、)( bankkkdxxlxf0)()( nkbakkdxxlxf0)()(Ak由由 决定，决定，与与 无关。无关。节点节点 f (x) bankjjxxxxkdxAjkj0)()(称为称为求积系数。求积系数。 knkkbaxfAdxxf 0)(求积公式求积公式定义定义其系数其系数 bakkdxxlA)(， 为拉各朗日插值基函数为拉各朗日插值基函数 这种求积公式称为这种求积公式称为插值型积分公式插值型积分公式)(xlk bankkxnbababankkkdxxxnfdxxRdxxxfxfAdxxffR0)1(0)()!1()()()(P)()()( 插值型的插值型的求积公式余项求积公式余项 为

8、了保证数值求积公式的精度为了保证数值求积公式的精度, ,我们自然希望我们自然希望求积公式能够对尽可能多的函数求积公式能够对尽可能多的函数f(x)f(x)都准确成立都准确成立, ,这在数学上常用这在数学上常用代数精度代数精度这一概念来说明。这一概念来说明。解：解：逐次检查公式是否精确成立逐次检查公式是否精确成立代入代入 f(x) = 1： baabdx1112 ab=代入代入 f(x) = x ：=代入代入 f(x) = x2 ： 222abdxxba 2baab 3332abdxxba 222baab 例：例：对于对于a, b上上1次插值，有次插值，有)()()(bfabaxafbabxxP

9、)()(2)(221bfafabdxxfabAAba 考察考察 时其求积误差。时其求积误差。mxxxxf, 1)(2梯形公式梯形公式因此梯形公式只对一次多项式精确成立。因此梯形公式只对一次多项式精确成立。代数精度代数精度定义定义 如果某个求积公式对于次数不超过如果某个求积公式对于次数不超过m m的一切多的一切多项式都准确成立项式都准确成立, ,而对而对 某个某个m+1m+1次多项式并不准确成次多项式并不准确成立，则称该求积公式的代数精度为立，则称该求积公式的代数精度为m m。显然，梯形公式与中矩形公式均具有一次代数精度。显然，梯形公式与中矩形公式均具有一次代数精度。 一般来说，代数精度越高，求

10、积公式越精确。一般来说，代数精度越高，求积公式越精确。定理定理 对于对于n+1n+1节点的插值型求积公式节点的插值型求积公式至少具有至少具有n n 次代数精度。（插值余项证明）次代数精度。（插值余项证明）例例试确定下面积分公式中的参数使其代数精确度尽量高试确定下面积分公式中的参数使其代数精确度尽量高. .)()()0()()0(2)()(120fIhffahhffhdxxffIh hdxxI00解解: :221hI 202231hahhI 0)(xxf 对对于于hI 1hhdxxI011)(xxf 对对于于22h hdxxI022)(xxf 对对于于33h 3)221(ha 1II 令令121

11、 a3022241hahhI hdxxI033)(xxf 对对于于44h 44h4023251hahhI hdxxI044)(xxf 对对于于55h 65h 3 , 2 , 1 , 0)()(1 jxIxIjj)()(414xIxI 因此因此所以该积分公式具有所以该积分公式具有3 3次代数精确度次代数精确度 5.2 牛顿柯特斯求积公式牛顿柯特斯求积公式Newton-CotesNewton-Cotes公式是指公式是指等距节点等距节点下使用下使用LagrangeLagrange插值插值多项式建立的数值求积公式多项式建立的数值求积公式,等等份份分分割割为为将将积积分分区区间间nbankkhaxk,

12、1 , 0, ,nabh 步步长长为为各节点为各节点为一、公式推导：一、公式推导： kknkkbaAxfAdxxf重点是求求积系数重点是求求积系数根据插值型求积公式根据插值型求积公式 0)(，以此分点为节点以此分点为节点,构造出的插值型求积公式。构造出的插值型求积公式。 bakkdxxlA)(dxxxxxbakjnjjkj 0:的计算的计算kA注意是等距节点注意是等距节点thax 令令,bax 由由, 0nt 可可知知kAdxxxxxbakjnjjkj 0dthhjkhjtnkjnj 00)()(dtjtknkhnkjnjkn 00)()!(!)1(00nj nj kathajhhdtakha

13、jh 000()()nj nj nj kj khtjdtkj 00( 1)()()! ()!n knkj nj kAbatj dtn knk )()(nkkCabA nkkknxfAfI0)()( nkknkxfCab0)()()(所以所以Newton-CotesNewton-Cotes公式为公式为Cotes系数系数()nkC注：注：Cotes 系数仅取决于系数仅取决于 n 和和 k，可查表得到。，可查表得到。与与 f (x) 及区间及区间a, b均无关。均无关。二、低阶二、低阶Newton-CotesNewton-Cotes公式及其余项公式及其余项在在Newton-CotesNewton-C

14、otes公式中公式中,n=1,2,4,n=1,2,4时的公式是最常用也时的公式是最常用也最重要三个公式最重要三个公式, ,称为低阶公式称为低阶公式(1). 梯形公式及其余项梯形公式及其余项abhbxaxn , 110则则取取101(1)1 0! 1!tdt )1(0CCotesCotes系数为系数为21 dtt 10)1(1C21 00( 1)()! ()!n knnkj nj kctj dtn knk )(1fI 10)1()()(kkkxfCab)()(210 xfxfab )()(2bfafab )(1fI即即上式即为上式即为梯形求积公式梯形求积公式, ,也称也称两点公式两点公式，记为，

15、记为)()(2)(bfafab )(1fIT 梯形公式的余项为梯形公式的余项为)()(1IRTR badxxR)(1求积公式为求积公式为dxbxaxfTRba )(2)()( dxbxaxfba )(2)( ,ba广义积分广义积分中值定理中值定理6)(2)(3abf )()!1()()(1)1(xnfxRnnn)(12)(3 fab 2312)(| )(|MabTR | )(|max,2xfMbax 故故(2).(2). 辛卜生公式及其余项辛卜生公式及其余项2,2, 2210abhbxabxaxn 则则取取CotesCotes系数为系数为dtttC 20)2(0)2)(1(4161 dtttC

16、 20)2(1)2(2164 dtttC 20)2(2)1(4161 求积公式为求积公式为)(2fI 20)2()()(kkkxfCab)(61)(64)(61)(210 xfxfxfab )()2(4)(6bfbafafab )(2fI-0.500.511.500.511.522.533.544.5上式称为上式称为辛卜生求积公式辛卜生求积公式，也称，也称三点公式或抛物线公式三点公式或抛物线公式记为记为)(2fIS SimpsonSimpson公式的余项为公式的余项为)()(2IRSR badxxR)(2)()2(180)4(4 fabab (3). 柯特斯柯特斯公式公式柯特斯公式柯特斯公式，

17、为四阶牛顿，为四阶牛顿取取4 n)(4fI)(7)(32)(12)(32)(79043210 xfxfxfxfxfab上式特别称为上式特别称为柯特斯求积公式柯特斯求积公式，也称，也称五点公式五点公式 柯特斯系数可由柯特斯系数表得到（柯特斯系数可由柯特斯系数表得到（P197P197）。由此可）。由此可以得到任意阶数的牛顿柯特斯求积公式。但实际计算时以得到任意阶数的牛顿柯特斯求积公式。但实际计算时一般不用高阶的公式，因为高次插值有一般不用高阶的公式，因为高次插值有Runge Runge 现象。现象。 由此可自然会得出以下结论：由此可自然会得出以下结论：1.1. 梯形规则简单，有梯形规则简单，有1

18、1阶代数精度；阶代数精度；2.2. 再增加一个节点，就是具有再增加一个节点，就是具有3 3阶代数精度的阶代数精度的辛卜生规则；辛卜生规则；三、牛顿科特斯公式的代数精度三、牛顿科特斯公式的代数精度牛顿科特斯公式实际上是插值求积公式，因此牛顿科特斯公式实际上是插值求积公式，因此n n阶牛顿科特斯公式至少有阶牛顿科特斯公式至少有n n次代数精度。次代数精度。由于牛顿科特斯公式是等距插值，因此，有定理：由于牛顿科特斯公式是等距插值，因此，有定理：当当n n为偶数时，牛顿科特斯公式有为偶数时，牛顿科特斯公式有n n1 1次代数精度次代数精度。四、复化求积法四、复化求积法固固定定时时而而节节点点个个数数的

19、的长长度度较较大大当当积积分分区区间间1, nba直接使用直接使用Newton-CotesNewton-Cotes公式的余项将会较大公式的余项将会较大。增增加加时时即即而而如如果果增增加加节节点点个个数数1, n公式的舍入误差又很难得到控制公式的舍入误差又很难得到控制为了提高公式的精度为了提高公式的精度, ,往往使用往往使用复化求积法复化求积法。分分成成若若干干个个子子区区间间即即将将积积分分区区间间,ba然后在每个小区间上使用低阶然后在每个小区间上使用低阶Newton-CotesNewton-Cotes公式公式最后将每个小区间上的积分的近似值相加最后将每个小区间上的积分的近似值相加 复化梯形

20、公式复化梯形公式：),., 0(,nkhkaxnabhk 在每个在每个 上用梯形公式：上用梯形公式：,1kkxx nkxfxfxxdxxfkkxxkkkk,., 1,)()(2)(111 11)()(2)(2nkkbfxfafh bankkkxfxfhdxxf11)()(2)(= Tn),(),()(12)()(12)(1221213bafabhnfabhfhfRnkknkk /*中值定理中值定理*/ 复化复化 Simpson 公式公式：),., 0(,nkhkaxnabhk )()(4)(6)(1211 kkkxxxfxfxfhdxxfkkkx21 kx1 kx44444 )()(2)(4)

21、(6)(1010121 nknkkkbabfxfxfafhdxxf= Sn)(2180)4(4 fhabfR 求积公式的余项比较求积公式的余项比较我们知道我们知道, ,两个求积公式的余项分别为两个求积公式的余项分别为)(TR)(12)(3 fab )(SR)()2(180)4(4 fabab 单纯的求积公式单纯的求积公式复化求积公式复化求积公式2()12bahf 4(4)()1802bahf 复化求积公式精度提高。复化求积公式精度提高。2 Composite QuadratureLab 13. Composite Trapezoidal RuleUse the Composite Trapez

22、oidal rule with a given n 0 to approximate a given integral .You are supposed to write a functiondouble CTR ( int n, double a, double b, double (*f)( ) )to approximate the integral from a to b of the function f using the trapezoidal rule on n equal-length subintervals.InputThere is no input file. In

23、stead, you must hand in your function in a *.h file. The rule of naming the *.h file is the same as that of naming the *.c or *.cpp files.Output For each test case, you are supposed to return the approximation of the integral. badxxf)(Sample Judge ProgramSample Judge Program#include #include #includ

24、e 98115001_13.h double f1 ( double x ) return (1.0/(1.0+sin(x)*sin(x); double f2 ( double x ) return (x*log(x); void main( ) FILE *outfile = fopen(out.txt, w); int n; double a, b;a = 0.0; b = 1.0; n = 10;fprintf(outfile, %lfn, CTR(n, a, b, f1);a = 1.0; b = 2.0; n = 4;fprintf(outfile, %lfn, CTR(n, a,

25、 b, f2); fclose(outfile);Sample Output Sample Output 0.8090930.639900nTI 2()12kbahf nSI 4(4)( )1802bahf 比较两种复合公式的的余项比较两种复合公式的的余项nnSTI即 比趋于定积分 的速度更快2 4h分别是 的 ，阶无穷小量为此介绍收敛阶的概念!)(2ho)(4ho定义定义1. 00,nnIpcII对于复合求积公式 若存在及使其余项满足0limnphIIchnIp则称复合求积公式 是 阶收敛的)(pnhII阶收敛的概念也等价于显然 p,不难知道,复合梯形、Simpson公式的收敛阶分别为2阶、

26、阶、4阶阶 复化求积法通过将积分区间分成复化求积法通过将积分区间分成n n等份，来减小截断等份，来减小截断误差，因此误差，因此n n越大积分精度越高。但越大积分精度越高。但n n太大，运算量也增大，太大，运算量也增大，舍入误差也增大；舍入误差也增大；n n太小，精度可能达不到。如何确定适太小，精度可能达不到。如何确定适当的当的, ,使得计算结果达到预选给定的精度要求呢？使得计算结果达到预选给定的精度要求呢？ 在实际计算中在实际计算中, ,常采用积分步长的自动选择。具体地讲常采用积分步长的自动选择。具体地讲, ,就就是在求积过程中是在求积过程中, ,将步长逐次折半将步长逐次折半, ,反复利用复合

27、求积公式反复利用复合求积公式, ,直直到相邻两次的计算结果之差的绝对值小于允许误差为止。这实到相邻两次的计算结果之差的绝对值小于允许误差为止。这实际上是一种事后估计误差的方法际上是一种事后估计误差的方法变步长求积算法变步长求积算法。5.3 变步长求积和龙贝格算法变步长求积和龙贝格算法 问题问题 综合前几节的内容综合前几节的内容, ,我们知道我们知道梯形公式梯形公式,Simpson,Simpson公式公式,Cotes,Cotes公式的代数精度分别为公式的代数精度分别为1 1次次,3,3次和次和5 5次次复化梯形、复化复化梯形、复化SimpsonSimpson、复化、复化CotesCotes公式的

28、收敛阶分别为公式的收敛阶分别为2 2阶、阶、4 4阶和阶和6 6阶阶5.3.1 变步长梯形求积法变步长梯形求积法 对于复合梯形公式，若将积分区间对于复合梯形公式，若将积分区间a,bna,bn等分等分, ,积分近似积分近似值记为值记为T Tn n , ,积分精确值记为积分精确值记为I I , ,则有则有: :把每个子区间分半把每个子区间分半, ,也就是将积分区间也就是将积分区间a,b 2na,b 2n等分等分, ,则有则有 bafnababTInnn ,)(12 2 bafnababTInnn 22 22,)2(12 则有则有当当 在连续在连续, ,且函数值变化不大时且函数值变化不大时, ,即有

29、即有)( xf)()(2 nnff 24nnITIT 给定求积精度给定求积精度 ，如何取，如何取n n ? ?244nnITIT2233nnnITTT223nnnTTIT 整整理理得得可用来判断可用来判断迭代迭代是否停止。是否停止。变步长梯形法计算过程变步长梯形法计算过程 |2nnTT迭迭代代停停止止条条件件： 22111afbfabTTn时，计算时，计算取取 112222122222 xfabTbafafbfabTTn 算算时，把区间二等分，计时，把区间二等分，计取取12abx 3 , 2 , 1,4, 421 22444 31232144 iabhihaxxfxfabTxfxfxfafbf

30、abTTni式中式中等分，计算等分，计算时，把区间时，把区间取取 可以看到可以看到, ,每次都是在前一次的基础上将子区间再对分。原分每次都是在前一次的基础上将子区间再对分。原分点上的函数值不需要重复计算点上的函数值不需要重复计算, ,只需计算新分点上的函数值即可只需计算新分点上的函数值即可, ,一般地计算公式为：一般地计算公式为：1201(21)222nnnibabaTTf ainn 由上节由上节变步长变步长梯形公式梯形公式得到的积分近似值的误得到的积分近似值的误差大致是差大致是 ,因此人们期望因此人们期望, 如果用这个误如果用这个误差作为对差作为对 的一种补偿的一种补偿,则得到的求积公式的代

31、数精则得到的求积公式的代数精度会有所提高。度会有所提高。)(312nnTT nT2nnnnnbaTTTTTdxxf3134)(31)(222 （1）5.3.2 龙贝格公式龙贝格公式 龙贝格算法是在复化梯形公式误差估计的基础上，龙贝格算法是在复化梯形公式误差估计的基础上，应用线性外推的方法构造出的加速算法。应用线性外推的方法构造出的加速算法。通过直接验证可知通过直接验证可知 也就是说也就是说 , 用梯形公式二分前后的两个积分值用梯形公式二分前后的两个积分值 与与 按照公式按照公式 (1) 做线形组合做线形组合,其结果正好是用抛物线其结果正好是用抛物线公式得到的积分值公式得到的积分值 。nTnT2

32、nS )()(2)(4)(6313410101221nknkkknnbfxfxfafhTTnSnnnTTS31342 (2)即即同理可知同理可知,用用抛物线公式抛物线公式得到的积分近似值得到的积分近似值 的误差的误差大致是大致是 ,因此对抛物线公式进行修正因此对抛物线公式进行修正, 得到得到nS215)(2nnSS nnnnnbaSSSSSdxxf1511516)(151)(222 (3)也就是说也就是说 , 用抛物线公式二分前后的积分值用抛物线公式二分前后的积分值 与与 按照公式按照公式(3) 作线形组合作线形组合, 其结果正好是用柯特斯公式其结果正好是用柯特斯公式得到的积分值得到的积分值

33、。nSnS2nC通过直接验证可知通过直接验证可知nnnSSC15115162 (4)同理可知，同理可知，用柯特斯公式得到的积分近似值用柯特斯公式得到的积分近似值 的误差的误差大致是大致是 ,因此因此, 对柯特斯公式进行修改对柯特斯公式进行修改,得得到求积公式到求积公式nC263)(2nnCC nnnnnbaSCCCCdxxf6316364)(631)(222 (5)为此为此,构造求积公式构造求积公式nnnCCR63163642 (6)称称(6)式为式为龙贝格龙贝格(Romberg)公式。公式。龙贝格公式是一种计算积分的方法。在变步长的求龙贝格公式是一种计算积分的方法。在变步长的求积过程中积过程

34、中, , 运用运用(2),(4),(6)(2),(4),(6)式可以将精度低的梯式可以将精度低的梯形值形值 逐步加工成精度较高的抛物线逐步加工成精度较高的抛物线 , ,柯特斯柯特斯值值 与龙贝格值与龙贝格值 。nTnSnCnR总之有：总之有：nnnSTT 1442nnnCSS 144222nnnRCC 144323Romberg 序列序列(1)计算计算f(a) , f(b) , 算出算出 。 1T(2) 把把 a , b 2等分等分, 计算计算 , 算出算出 与与 。 2baf2T1S(3) 把把 a , b 4等分等分,计算计算 算出算出 与与 。,43,4 abafabaf24,ST1C龙贝格求积的计算步骤如下：龙贝格求积的计算步骤如下：(4) 把把 a , b 8 等分等分, 计算计算 算出算出 与与 与与 。7, 5, 3, 1,8 iabiaf48,ST1R2C(5) 把把 a , b 16等分等分, 计算计算15, 5,