数值积分与数值微分

1、a1a2在工程问题和科学实验中，常常需要计算积分。例如：力学和电学中功和功率的计算，电流和电压的平均值和有效值的计算以及一些几何图形的面积、体积和弧长的计算等等。另外微分方程的求解也是以积分计算为基础的。在高等数学中，曾用Newton-Leibniz公式 xfxFaFbFabxFdxxfba（其中 是 的一个原函数）来计算定积分，虽然这一公式比较简单，但常常会遇到下面情况：的原函数不能用初等函数表示。 如 其原函数不能用初等函数表示。 xfxfxf、321dxedxxSinxx10102,的结构复杂，求原函数困难。的表达式不知道，只给出了一张由试验提供的函数表。a3利用积分中值定理：即以底长b

2、-a，高为 的矩形的面积恰等于所求曲边梯形的面积I. fbafabdxxfba,.xfy fxba0这样，只要只要对平均高度 提供一种算法，便获得便获得一种数值积分方法。如果用两端点的“高度” 与 取算术平均作为平均高度 的近似值 ，这样导出积分近似公式： bfaff 12bafTf af b对于这些情况，要计算积分的准确值都是十分困难的，这就要求建立积分的近似计算方法。此外，积分的近似算法又为其它一些数值方法，例如微分方程数值解、积分方程数值解等，提供了必要的基础。问题是 的具体位置一般是不知道的，因而难以准确算出 的值。称 为区间a,b上的平均高度平均高度。 f fa4kNkkNbaxfd

3、xxf)(lim)(1一般地，由定积分定义： 积 分为和式的极限)2()2()(bafabR的“高度”f(c)近似取代高度 ， 则导出中矩形公式中矩形公式： 这是熟知的梯形公式梯形公式（如下图）而若改用区间中点 2bac)(fabxy)(xfa5 下面通过最简单的情形做进一步分析。如果把区间a,bN等份。节点：(0,1,1)kkkkxakhbahxx kNN取定积分的基本分析步骤是四步：分割、近似、求和、取极限。分割就是把总体（整块梯形面积）分成若干分量（小曲边梯形面积），近似就是在每个分量中用容易计算的量去代表小曲边梯形的面积（这是用矩形面积近似曲边梯形的面积）。求和就是把分量加起来得到总近

4、似值，最后取极限就得到积分的准确值。以上分析看出，前三步比较容易，最后一步的计算比较困难，但现在只是要求近似值，因而可以省掉求极限这一步，只要经过前三步就可求得积分 近似值，这就是建立数值积分方法的基本步骤。a6xy)(xfy ax 01x2x3xbx 4hkaxxkkk)21(21若在同一分割下取,11() ( )2kkkbax xf akhf xxN就是小区间中用近似 即用小矩形面积近似曲边梯形面积。如下图示。 )3()()(10NkbakhafNabdxxf所有相加即得kkkxxfNabxf)()(在每个小区间上这称为左矩形公式左矩形公式。a7把每个小区间的近似值相加就得到新的积分公式：

5、 称为中矩形求积公式中矩形求积公式1( )() )(4)2bab af x dxf akhN11101( )()1( ) ( ( )( )( ) (5)22NNbkkkakkf xf xb ab af x dxf af bf xNN 再把每个小区间的分量相加就得梯形求积公式梯形求积公式： 若在每个小区间 中用梯形面积近似曲边梯形面积 ，,1kkx x2)()(1kkxfxfNab即用kxxf)(近似a8 由上述知，计算积分的问题归结为计算被积函数 f(x) 在 a,b 上某些点的值 ，这就使积分的计算变得很容易了。现在的问题是如何 提高精确度。从求近似积分的三步看来，前两步分割与近似是从求近似

6、积分的三步看来，前两步分割与近似是提高求积公式精确度的途径。第三步求和无须专门讨论提高求积公式精确度的途径。第三步求和无须专门讨论。)(kxf 先讨论不分割区间a,b 时，哪一种近似得到的求积公式精确度更高。事实上，对于不分割区间 而言，用零次多项式零次多项式 近似f(x), 然后取积分而得到的矩形公式1, 0 xx)(0 xf 直观上看，它比左矩形公式（3）精确些，但总的说来，上述矩形法和梯形法精度都较低，往往不能满足实际计算的要求，因此需要建立精确度更高的求积公式。 分割主要是细分，从理论上讲，分得越细，精确度越高，但必然增加计算量，增大舍入误差。因此，怎样在一种合适的分割下，使用容易计算

7、的分量去近似曲边梯形面积就成为提高求积公式精确度的主要问题。a9)()()(01001010 xxxfdxxfdxxfxxxx同样，用一次多项式近似一次多项式近似f (x), 取积分而得到梯形公式：)(2)()()()()()()(01100010101010 xxxfxfdxxxxxxfxfxfdxxfxxxxkKkxAx为 求 积 节 点 ，为 求 积 系 数 ， 亦 称 之 为 伴 随 节 点的 权 。 而不依赖于被积函数f(x) 的具体形式。我们的目的就是根据一定原则, ,kkAx 和系数选择求积节点使得求积一般公式（6）具有较高的精确度同时又计算简单。)6()()(0kNkkbaxf

8、Adxxf,)()(babadxxfdxxp的近似值的值作为积分可以设想，用更精确的插值多项式更精确的插值多项式 p(x)近似近似 f (x), 然后以或许会提高近似积分的精确度，这就是我们建立求积 公式的基本思路。这样得到的公式一般可以写成:,的选取有关仅仅与节点权kkxAa10 数值求积方法是 近似方法，为要保证精确度，自然希望求积公式能够对“尽可能多”的被积函数f(x)都准确成立，在计算方法中，常用代数精度这个概念来描述它。)()(0knkKbaxfAdxxf定义定义1：若求积公式却不一定能准确成立，则称该求积公式的代数精度则称该求积公式的代数精度m。1mx而对于对于任意不高于m 次的代

9、数多项式都准确成立，上式右端=即f(x)=1时上面求积公式准确成立。 )()(2)(bfafabdxxfbaabdxba1abab 11 2)(2122abxdxba)(21)(2122ababab例如：梯形求积公式的代数精度m =1，这是因为：当 f(x)=1时，上式左端= 上式右端=当f(x) = x 时，上式左端=a11)()(0knkkbaxfAdxxf02201101()(7)21()1nkknkkknmmmkkkAbaA xbaA xbam 一般的，欲使求积公式 具有 m 次代数精度。都能准确成立mxxxxf, 1)(2只要令它对于即要求右端=上式左端= 综上即得此求积公式对f(x

10、)=1和f(x) = x 的任一线性组合，即不高于 1次的代数多项式都准确成立。即此求积公式代数精度至少为m=1。)(31332abdxxba)(222baab即 f(x) = x 时，上面求积公式也准确成立。 时但是当2)(xxf 当 ab 时，左端右端 因此，由定义知梯形求积公式代数精度为m=1。 同理可证：矩形求积公式也具有一次代数精度。a12如果事先选定求积节点，如，以区间a,b的等距节点依次为节点，这时取m=n，求解上述线性方程组（）即可确定系数从而使求积公式至少有m=n次代数精度。具体示例在下面一节中介绍。kA由插值余项定理，对于插值型的求积公式其余项为：)8()()()()(11

11、dxxwxxxwdxxlAbanknbakk若记(1)1( )( ) ( )( )( )(9)(1)!nbbnnnaafR fIIf xL x dxx dxn 有关与其中nxxxx,10)()()(101nnxxxxxxx已知), 1 , 0)(,10nkxfbxxxakn)()()(0 xlxfxLknkkn), 1 , 0)(nkxlk其中0( )( )()( )nbbbnkknaaakIf x dxL x dxf xlx dx I设给定一组节点作 f(x)的n次Lagrange 插值公式 为n 次插值基函数，代替被积函数得：用)(xLna13定理定理1： 的求积公式至少有n次代数精度的)

12、(0nkkknxfAI形如充分必要条件充分必要条件是：它是插值型的。因此，对于插值型求积公式，由上面余项公式（9）可见，对于次数n 的多项式f(x)，其余项R(f)=0，因而这时求积公式至少具有n 次代数精度。至少具有n次代数精度，则它必是插值型的，这是因为： 若它具有n次代数精度，则对于n 次插值基函数 应准确地有下式成立：)(0knkknxfAI)(xlk)()(0jknjjbakxlAdxxl反之，若求积公式kAjkjkxlkjjk, 0, 1)(再注意到，即知上式右端实际上等于因而（8）成立。故为插值型求积公式。综上所述有： a14则令为确定khaxthaxcknk,)()!( !)

13、1()()2)(1(1)2)(1()()()()()()2)(1()()()(1110,1111knkhnkkkkhxxxxxxxxxxxwntttthxxxxaxxwnknnnkkkkkkknnnnk)(nkcnabh步长khaxk选取等距节点) 1 ()()(0)(nkknknxfcabI设将积分区间a,b 划分n 等份，构造出的插值型求积公式：称作称作Newton-Cotes 公式公式，其中 称作称作Cotes 系数系数。a15( ),11( )()()bnkkakkxcAdxbabaxxx则，,1100( )()()( 1)(1)()() !()!( 1)(1)(1)(1)()!()!

14、bkakkn knnnn knxAdxxxxht ttnhdthtk k nkht ttktktn dtk nk而，( )000,11( 1)(1)(1)(1)()!()!( 1)()(2)!()!nkkkn knn knnjj kcAAbanht ttktktn dtnk nktj dtnk nk则，a16211)1(1)2(0CCn时，当),) 1(10)1(110)1(0tdtCdtkC即为梯形公式）：因此由（)3()(21)(21)()(1101xfxfabIdxxfba)4()(61)(64)(61)()(61) 1(4164)2(2161)2)(1(412210220)2(220)

15、2(120)2(0 xfxfxfabIdxxfdtttCdtttCdtttCnba时，当由于是多项式积分，Cotes 系数计算不会遇到实质性困难。a1740)4(240)4(3)4(140)4(4)4(09012)4)(3)(1(1619032)4)(3)(2(! 41907)3)(2)(1(! 4414dtttttCdtttttCCdtttttCCn时，当)(907)(9032)(9012)(9032)(907)()(1432104xfxfxfxfxfabIdxxfba）有：由（这就是抛物线公式，又称辛浦生辛浦生 Simpson 公式公式 。几何意义就是用抛物线下的面积近似曲线f(x)下的面

16、积。a18这就是柯特斯公式柯特斯公式（Cotes ) 当n 较大时，例如 n=8 时，系数 中出现负数，而且有正有负会使舍入误差增大，数值稳定性较差，因此实际计算并不用 n较大的 公式，而是将区间a,b 分割成若干个小区间，对每个或几个小区间应用n 较小的公式去计算。Cotes系数表详见教材.)(nkC先看Simpson公式，它是二阶Newton-Cotes公式, 因此至少具有二次代精度。进一步用 进行检验。 按Simpson公式计算得：3)(xxf)2(46333bbaaabS作为插值型求积公式，n阶Newton-Cotes公式至少具有n次的代数精度（定理1），那么实际的代数精度可否进一步提

17、高呢？a19另一方面，直接求积分得：5)2(46554444abdxxIbbaaabSba4443abdxxIba4)(xxf而对 易验证S=I, 即Simpson 公式对次数不超过三次的多项式均能准确成立。 易验证此时SI （ b a 时） 因此Simpson公式实际上具有3次代数精度。一般的有： 定理定理2：当阶n 为偶数时， Newton-Cotes公式至少有n+1 次代数精度。a20dxxxdxxdxxnfdxxLxfdxxLdxxffRnxfjbanjbannbannbababann)()()()!1()()()()()()()!1()(011)1()1(而由于引进变换x=a+ht

18、, 并注意到有：jhaxjdtjthfRnnjn)()(0021( )nf xx的余项为零。证明：只要证明当n为偶数时， Newton-Cotes公式对a21上面已说明 Newton-Cotes公式的余项为 ),(,)()!1()()()()(1)1(badxxnfdxxLdxxffRnbannbaba为整数则2n，进一步有：再令2nutdujnuhfRnnnjn)2()(2202若 n 为偶数，而积分区间关于原点对称，因此积分为0，即 证毕。为奇函数由于被积函数)()2()(220jujnuuHnnnj0)(fR（n+1) 个乘积,(u-j)(u+j)为偶函数， j=0时，u-j 为奇函数。

19、a22dxbxaxffRba)()(! 21)(, ,)(, ,f由于0)()(bxaxx当当n=1 即梯形公式的余项为即梯形公式的余项为：是依赖于x的函数，在 a,b上连续，而且使：),(ba)()(61)()()()(, ,3, , ,fabdxbxaxfdxbxaxfbaba),(),(12)()(, ,3bafabfR故可应用中值理，从而梯形公式的截断误差为：当当n=2 即即Simpson公式的余项公式的余项。为此，先构造一个次数的插值多项式H(x) ，使满足：3a23)()(),()(),()(),()(,cfcHcfcHbfbHafaH2bac其中)()2(4)(6)()(4)(6

20、)(bfbafafabbHcHaHabdxxHba 由于已证明Simpson公式具有三次代数精度，因此， Simpson公式对于这样构造的三次插值函数H(x) 是准确的：即（由插值条件）此即为Simpson公式求得的f(x)的积分。因此Simpson公式余项为：dxxHxfSIfRba )()()( 又可以证明：满足上面插值条件的多项式H(x)的插值余项为：a24 )()(! 4)(2) 4(bxcxaxfxHxf,故积分余项为：dxbxcxaxffRba)()(! 4)(2)4(当n=4，即Cotes公式的积分余项为： (9)(4945)(2)6(6fababCIfR )(218044fab

21、ab由于 是依赖于x 的函数，在a,b上连续，且)()4(f 0)()(2bxcxaxxdxbxcxaxffRba)()(! 4)(2)4(故可利用中值定理： （8)a25 由上面截断误差易见，当积分区间a,b较大时，直接使用N-C公式的截断误差增大，积分近似值的精度难以保证。因此，在实际应用中，为了既提高结果的精度，又使算法简便且易在电子计算机上实现，往往采取复合求积的方法复合求积的方法，即：先将积分区间分成几个小区间，并在每个小区间上利用低阶Newton-Cotes公式计算积分近似值。然后对这些近似值求和，从而得所求积分的近似值。由此得到一些具有更大实用价值的数值求积公式，统称为复合求积复

22、合求积公式公式。 然后在每个小区间 上应用梯形公式即：kxakhnabh11( )()()2kkxkkxhfx dxfxfx例如：先将区间a,bn等分，记分点为：（k=0n)，其中步长为（k=1,2n)1,kkxxa26可导出复合梯形公式：注意 (10) )()(0afxf)()(bfxfnnnkkknkkbankxxTbfxfafhxfxfhdxxfdxxfkk)()(2)(2)()(2)()(111111若在 上应用Simpson公式：,1kkxx111( )()4()62kkxkkkkxxxhf x dxf xff x从而可得复化Simpson公式：a27 同理可得复化Cotes公式)(

23、)(2)(4)(6)(111021bfxfxfafhSdxxfnkknkkban11104310211041)(7)(14)(32)(12)(32)(790)(nkknkknkknkkbanbfxfxfxfxfafhCdxxfhxxkk4141hxxkk2121hxxkk4343注：（12）（11）a28定理定理：若f(x)在积分区间a,b上分别具有二阶、四阶和六阶连续导数，则复化梯形公式、复化Simpson公式和复化Cotes公式的余项分别为：其中a,b )()4(945)(2)()()()2(180)()()( 12)()()6(6)4(42fhabCdxxffRfhabSdxxffRfh

24、abTdxxffRbancbansbanT（13）（14）（15）)()()4(9452)()( )( )2(1801)()()(121)()5()5(642afbfhfRafbfhfRafbfhfRcsT（16）（17）（18）当 h 充分小时又有：a29证明：只证明复合梯形公式的余项，其余由学生类似可做。 由于f (x)在a,b连续，从而在每个小区间 上做积分 使用梯形公式时的截断误差为：（前面已证）,1kkxxkkxxdxxf1)(而 ，)( )( )( 112)(212nTfffnhabfR)( max)( )( )( 1)( min,21,xffffnxfbaxnbax故 （*）)(

25、 )( )( 121)()(213nbanTfffhTdxxffR即),(1kkkxx)( 12)(31kkkfxx1kkxxh（ ），而a30)( )( )( 1)( 21nfffnf故 # )( 12)(2fhabfRT)( )( 121)( 121)( 12lim)(lim1020afbfdxxfhfhTdxxfnkbakhbanh)( )( 121)(2afbfhfRTf(x)在a,b上的连续性及介值定理知：有(a,b),使又由上面（*）式及定积分定义有：可见，当h充分小时有： 其余两类公式证明与此完全类似，略。 a31由上定理中余项及近似余项公式可见，只要所涉及的各阶导数在积分区间上

26、连续，则当n （即h 0)时， 、 和 都收敛于积分真值 （所谓积分真值是指每一位数字都是有效数字的积分值），而且收敛速度一个比一个快。nTnSnCbadxxf)(则称则称 是是P阶收敛的阶收敛的。定义定义：对于复合求积公式 ，若当h 0时有dxxfIban)(0)(cchIdxxfpbannIa32 对于数值求积公式来说，收敛阶越高，近似值 收敛到其值 的速度就越快，在相近的计算工作量（顺便提一下，数值求积计算工作量的大小，主要取决于计算函数值次数的多少）下，有可能获得较准确的近似值。nIbadxxf)(例：利用复合Newton-Cotes公式计算 的近似值。10214dxx定理：复合求积公

27、式：梯形、Simpson、Cotes分别具有二阶、 四阶和六阶收敛性证：复合梯形公式的收敛性由上面证明易见。 其余两个也可类似证明。 #a33解：用两种方法求解。先将积分区间0，1八等分（分点及分 点处的函数值见下表），用复合梯形公式得：（h=1/8）138988. 3) 1 (87868584838281 2) 0 (1618fffffffffTxx04.00000005/8 2.876404491/83.938461546/8 2.560000001/43.764705887/8 2.265486733/83.5068493212.000000001/23.20000000)1/(4)(2

28、xxf)1/(4)(2xxf141593. 3) 1 (434241 287858381 4)0(6414fffffffffS再将区间0，1四等分，利用复合Simpson公式得：a34两种方法都用到表中九个点上的函数值，计算工作量基本相同。但所得计算结果与积分真值=3.14159265相比较，复合Simpson公式所得近似值 远比复合梯形公式所得近似值 要精确。因此在实际计算中，较多的应用复合Simpson公式。为了便于上机计算，常将复合Simpson公式改写成： 4S8TbankkknxfxfbfafhSdxxf121)()(22)()(6)(19)a35 虽然,我们上面已得到Newton-

29、Cotes低阶公式的近似值的误差估计，也可根据精度要求用这些公式确定积分区间的等份数，即确定步长 h ，但由于余项公式中含有被积函数 f(x) 的高阶导数，在具体计算时往往会遇到困难，因此，在实际应用时，常常利用误差的事后估计法来估计近似值的误差，或确定步长 h。此方法的大致做法是：将积分区间逐步将积分区间逐步分半，每分一次就用同一复合求积公式算出相应的积分近似分半，每分一次就用同一复合求积公式算出相应的积分近似值，并用前后两次计算来判断误差的大小。值，并用前后两次计算来判断误差的大小。a36 其原理和具体做法如下：nT2 对复化梯形公式（10），由其余项公式（13）或（16）可见，当f (x

30、)在积分区间上变化不大或积分区间a,b的等分数n较大（即步长h较小）时，若将a,b的等分数改为2n（即将步长缩小到原来步长的一半），则新近似值 的余项约为原近似值余项的1/4，即：（令 ）badxxfI)( 此式表明：若用 作为积分真值I的近似值，则其误差约为412nnTITI)(3122nnnTTTI （20）即有nT2)(312nnTTa37先算出 和 ，若 （为计算结果的允许误差），则停止计算，并取 为积分近似值。否则，将区间再次分半后算出新的近似值 ，并检验不等式 是否成立，直到得到满足精度要求的结果为止。 nT2nT232nnTT324nnTTnTnT4 故将区间逐次分半进行计算的过

31、程中，可以用前后两次计算结果 和 来估计误差和确定步长。具体做法是：nT2nT对于复合Simpson公式（11），复合Cotes公式（12），由它们的积分余项（14）、（15）或（17）、（18）可见，当所涉及的高阶导数在积分区间上变化不大或积分区间的等份数n较大时，a38有： 和 1612nnSISI6412nnCICI则 （21） 及 （22）)(15122nnnSSSI)(63122nnnCCCI 因此，也可以像使用复合梯形公式求积分近似值那样，在将积分区间逐次分半进行计算的过程中，估计新近似值 和 的误差，并判断计算过程是否需要继续进行下去。上述过程很容易在计算机上实现。nS2nC2a

32、39 先看复合梯形公式（10），计算 时，需计算n+1个点（它们是积分区间a,b的n等分的分点）上的函数值，当 不满足精度要求时，根据上面提供的方案，就应将各小区间分半，计算出新近似值 。若仍用（10）计算 ，就需求出2n+1个点（a,b的2n等分点）上的函数值。nTnTnT2nT2 上面介绍的步长变化的计算方案，虽然提供了估计误差与选取步长的简便方法，但是还没有考虑到避免在同一节点上重复计算函数值的问题，故有进一步改进的余地。a40 而实际上，在这2n+1个分点中，包含有n+1个n分点，对应的函数值在计算 时已算出，现重新来计算这些点上的函数值，显然是极不合理的。nT 为避免这种重复计算，我

33、们来分析新近似值 与原有近似值 之间的关系。由复合梯形公式（10）知：nT2nT1212)()2(2)(4nknbfnabkafafnabT注意到在2n分点 （k=1,2,2n-1)中，当k取偶数时， 即为n分点，k为奇数时， 才是新增加的分点，将新增加的分点处的函数值从求和记号中分离出来，就有: nabkaxk2kxkxa41 由递推复化梯形公式（23）可见，在已计算出 基础上再计算 时，只要计算n个新分点上的函数值就行了，这与直接利用复合梯形公式（10）相比，计算工作量几乎节省一半。nTnT2例例2 2：由（23）重新计算 的近似值，误差不超过 。10214dxx6101211111( )

34、 222(21)( )422 ( ) 2()( )(21)422nnnkknnkkb ab ab aTf af akf akf bnnnb ab ab ab af af a kf bf aknnnn (23)11111 ( ) 2()( )(21)2 2221(21) )222nnkknnkhb ab af af a khf bf aknnhhb aTf akhn 其中a42 再将各个小区间二等分，出现两个新分点x=1/4, x=3/4,由 （23）得：13117647. 3)43()41(412124ffTT解：在积分区间逐次分半的过程中顺序计算积分近似值 ， 并用是否满足 （= ）来判断计

35、 算过程是否需要计算下去。1T2T84,TT32nnTT6101 . 3)21(212112fTT3)24(21)1 ()0(211ffT 然后将区间二等分，出现新分点是x=1/2,由递推公 （23） 得： 先对整个区间0，1使用梯形公式得：a43为了便于上机计算，我们将积分区间a,b的等分点依次取为： 如上例，并将递推公式改写为：,2 ,2 ,22101212212) 12(221)()(21kikkabiafabTTbfafabTkkk=1,2,3, （24） 这样不断将各小区间二分下去可利用递推公式（22）依次算出计算结果如右表，因为：故 =3.14159202是满足要求的近似值。321

36、68,TTT6256512103 TT512Tnn13323.1414298923.1643.1415519643.131176471283.1415824883.138988492563140941615123.14159202nTnTa44 Romberg算法是在积分区间逐次分半的过程中，对用复合梯形法产生的近似值进行加权平均，以获得准确度较高的近似值的一种方法，具有公式简练，使用方便，结果较可靠的优点。 上面介绍的递推公式（23）或（24），虽然具有结构简单，易在计算机上实现等优点，但是由它产生的梯形序列 ，其收敛速度却非常缓慢。如上例中用（23）计算 的近似值

37、时，要一直算到 才获得误差不超过 的近似值。因此，用这种方法计算更复杂的高精度要求的近似值，显然是费时、费力甚至是不可能的 。10214dxx512T610kT2a45 如何提高收敛速度，节约计算工作量，自然是我们所最 关心的问题。 nnnTTT31342（1） 有可能比 更好地接近于积分 的真值 InT2 dxxfba 如上例中，13117647. 34T 和17898849. 38T 为两个精度很差的近似值，但如果将它们按（1）作线性 组合， 所得到的新近似值：14159250. 33134484TTT由（20）易见，nT2作为积分真值I的近似值，其误差约为)(312nnTT 。 如果它作

38、为nT2 的一种补偿，则可期望新)(3122nnnnTTTT即近似值。a46的、却具有7位有效数字，其准确度比 还高。但计算512T4T8T仅涉及到9个点上的函数值，其计算工作量仅为计算512T571。与这表明在收敛缓慢的梯形数列 的基础上，若将 kT2nT2nT按（1）作线性组合就可产生收敛速度较快的Simpson序列：kS21S：、2S、3S关于（2）的证明：由（23）可知：这种新近似值 实质上又是什么呢？可以验证：nTnnST 即：144313422nnnnnTTTTS （2）a4731211223412nknnnnabkafnabTTT 故（2）式成立 nknkkhkafhbfxfaf

39、h1112122231 nnknkkkSbfxfxfafh111021426(由上节（10）（由上节（11）a4814463163643232nnnnnCCCCR（4） 由近似式（22） 类似推导可得：)(63122nnnCCCInnnSSSI22151144151151622222nnnnnSSSSC同理，由上节近似式（21）类似推导可得：（3）即将Simpson序列 按（3）作线性组合就可产生收敛kS2速度更快的新序列,:4222CCCCka49可以证明可以证明：当f（x）满足一定条件时，Romberg序列 kR2比Cotes序列 能更快地收敛到积分 的真值I。kC2 dxxfba即在Co

40、tes序列 的基础上，产生了一个称为RombergkC2 :2kR1R、2R、4R.序列的新序列序列的新序列越精确的近似值 也就是将收敛速度缓慢综上可知：在积分区间逐次分半过程中利用公式可以将粗糙的近似值 逐步地“加工”成越来的梯形序列 逐步地“加工”成收敛速度越来越快的新 kT2新序列 4,3,2,nnnRCS,222kkkRCSa50 例例3 3: 利用公式(2).(3).(4)“加工”例2中得到的 近似值dxx10214.,8421TTTT解:依公式(2),(3),(4)按上图中计算顺序可得结果如下表:其中k为二次分数.这种加速的方法称为 Romberg算法。其“加工”过程如下图，其中圆

41、圈中号码表示计算顺序。 1T2T4T8T16T1S2S4S8S1C2C4C1R2Ra5112kS22kC32kR注意 教材中介绍的Richardson外推法外推法，为便于上机计算，引用记号 来表示各近似值，其中k仍代表积分区间的二次分数，而下标m则指出了近似值 所在序列的性质。如m=0在梯形序列中，m=1在Simpson序列中，m=2在Cotes序列中， 引入上面记号后，Romberg算法所用到的各个计算公式可统一化为： kmT kmT可见，“加工”效果十分显著，而“加工”计算量因只需 做少量四则运算没有涉及到求函数值，故可忽略不计。0 3k1 3.1 3.1333332 3.131176 3

42、.141596 3.14121183 3.138988 3.141593 3.141594 3.141588kT2a52由此可逐行构造 出一个三角形数表-称为T数表K0123. 5 bfafabT200 3 , 2 , 1,212221121100labiafabTTliilll 3 , 2 , 1, 2 , 1 , 0,144111mkTTTmkmkmmkm kT011kT22kT33kT 00T 10T 20T 30T 01T 11T 21T 02T 12T 03Ta53 ITkmklim ITmk0lim(m固定），可以证明：若 f（x）充分光滑，则T数表每一列的元素及对角线均收敛到所求

43、的积分值I。即：这是因为： 时，梯形积分公式余项可写成： baCf,2462123( )(1,2,3,)kkkT hIhhhhk其中系数均与步长h无关.若令: hThThThThTmmmmmm1101412144,经过m（m=1，2，)次加速后，余项取下面形式：2(1)2(2)12( )mmmT hIhh 此即为Richardson外推法外推法。其中 由（5）即知上面两个极限成立。可见加速后余项数量级下降很快。kabh2a541、在上面“加工”过程中的系数 和 ，当 m 4 时， 而另一个 系数则接近于1，也就是144mm141m2551141m新公式与原公式差别不大， 但工作量却大增。 km

44、kmTT1因此，在实际计算中常规定 m 3，即计算到出现Romberg序列为止。 2、 可用二维数组来存放并参加运算，也可用一维数组。 kmTa553、 对于积分限为无穷的积分 ，可利用变量代换化成有限区 间的积分然后再进行计算。例如：4、若被积函数有奇异点（间断点）存在于积分区间内，则 可得积分 区间分成小部分，使间断点在子区间的端点处。 也可用变量代换法处理。 1321xxdxWFM81.51tx1令则dttdx2代入得：1030132213211111tttdttttdxxxdx a56有可能使求积公式具有使求积公式具有2n+1次代数精度次代数精度，这类求积公式称为 Gauss公式公式、

45、相应的节点 为Gauss点点。), 1 , 0(nkxk机械求积公式含有2n+2个待定参数 ，适当选择这些参数 kbankkxfAdxxf0), 1 , 0(,nkAxkk 1 前面我们限定把积分区间的等分点作为求积节点，从而构造出一类特殊的插值求积公式，即NewtonCotes公式。这种做法虽然简化了算法，但却降低了所得公式的代数精度。例如：在构造形如 的两点积分公式时，如果限定求积节点 ，那么 所 得插 值型求积公式 ，其代数精度仅为1。 110011xfAxfAdxxf(*) 1111ffdxxf011,1xx a57 但是，如果我们对公式（*）中系数 和节点 都不加 限制，那么就可以适

46、当选取 使所得公式的代数精 度 m1。 事实上，若要求（*）对f（x）=1，x， 都准确成立，1, 0AA10,xx32,xx1010,AAxx和210AA01100 xAxA32211200 xAxA0311300 xAxA满足方程组：1010,AAxx和只要易验证，这是代数精度为m=3的插值型求积公式。解之得：33,33, 11010 xxAA代入（*）得：1133( )()()33f x dxffa58 证明： 设p(x)是任意次数不超过n的多项式。则 p(x)w(x) 的次 数不超过2n+1。因此，如果 是Gauss点，则求 积分公式（1）对 p(x)w(x) 能准确成立。即有：nxx

47、x,10 0dxxwxpba 2与任意次数不超过n的多项式P（x）均正交，即： nkkxxxw0 定理：对于插值型求积公式（1），其节点 是 Gauss节点的充分必要条件充分必要条件是：以这些点为零点的多项式), 1 , 0(nkxk 因此，可以想象对机械求积公式（1），只要适当选择2n+2个待定参数 ，是可以使其代表精度达 到 2n+1的。可以证明，Gauss型求积公式是插值型的。), 1 , 0(,nkAxkka59 对于任意给定的次数不超过2n+1的多项式f（x），用w（x）除 f（x），记商为p(x)，余式为Q（x），P（x）与Q（x）都是次数不超过n的多项式。 f(x)=p(x)w(

48、x)+Q(x),由 而求积公式（1）是插值型的，代数精度至少为m=n，因此 对Q（x）能准确成立： kknkkbaxwxpAdxxwxp0nkxwk, 1 , 0, 0但故（2）成立。kkkxQxfxw 0由（2）得： dxxQdxxfbaba knkkknkkbaxfAxQAdxxQ00a60（2）、再利用Gauss点确定求积系数), 1 , 0(nkAk 利用解方程组求Gauss点 和权 以确定Gauss公式，需 解2n+2个未知数的方程组，工作量太大。比较简便的做法是：kxkA也就是求积公式（1）对一切次数不超过2n+1的多项式均能 准确成立，因此为Gauss点（k=0，1，n)。 kn

49、kkbaxfAxf0即（有关Legendre多项式见教材！ ）（1）、先利用 上的n+1次正交多项式确定Gauss点ba, , ,0,1,2,kxa bkna61 令它对f（x）=1准确成立。即可得出 。这样构造出的一点Gauss-Legendre公式是中矩形公式。若取 的零点 作节点构造求积公式 xxp100 x 0011fAdxxf20AGaussLegendre公式余项为： 122!2!2!222nnnnffRnn1 , 1SG28 .51371P由于Legendre多项式是 的正交，因此可取Legendre多 多项式 的零点作为求积分公式（*）的高斯点，形如 （*）的Gauss积分公式

50、特别地称作GaussLegendre公式 。 knkkxfAdxxf011 *1 , 1 xpn 1先考虑前面的例子：再取 的两个零点 构造求积公式：31221( )(31)2P xxa62 31311011fAfAdxxf令它对f（x）=1，x均准确成立，即210AA0313110AA110AA从而得两点GaussLegendre公式： 313111ffdxxf类似地，三点GaussLegendre公式形式为： 515950985159511fffdxxfGaussLegendre公式节点和系数详见教材。a63对于一般区间对于一般区间 上的积分上的积分，也可以利用教材中表格的数据写出Gaus

51、s型求积公式。其原理与方法是其原理与方法是：先作变量替换，令 tabbax22则将 上积分化为 上的积分：ba,1 , 1 dttabbafabdxxfba11222 A tabbaftg22记 ，则上式化为： dttgabdxxfba112 B利用表中数据，对于给定的n=1，2，3，4，可以写出Gauss型公式： knkktgAdttg011 Cba,a64knkktabbafAdttabbaf2222011即 代入（A）得： knkkbatabbafAabdxxf2220 D其中系数 和节点 可查表得出，由变量替换公式kAkttabbax22 易见，由于求积公式（C）对变量t不高于2n+1

52、 的多项式准确成立，从而求积公式（D）对自变量x的不高于2n+1 的多项式也准确成立，即（D）是Gauss型求积公式。 例例：利用四点Gauss求积公式计算 的近似值。dxx10214解：由教材中表示Gauss型求积公式（D）得： kkkbatabbafAabdxxf22230a65其中a=0，b=1， 214xxf03, 12tttt0t= -0.86113631，1t= -0.33998104，65214515. 0,34785485. 01203AAAA将上述各数据代入上公式中有：141624.314102dxx优点优点：在此例计算过程中，只涉及到 四个点上的函数值，可见Gauss型求积

53、公式具有计算工作量小，所得近似值精确度高的优点，是一种高精度的求积公式。 Gauss型求积公式的缺点缺点是：当n改变大小时，系数和节点几乎都在改变。虽然可以通过其他资料查到较大n的系数和节点，但应用时却十分不便。同时，余项却涉及到高阶导数（被积函数的），要利用它们来控制精度也十分困难。a66 为克服这些缺点，在实际计算中较多地采用复合求积的方法。例如，先把积分区间 分成m个等长的小 区间 然后，在每个小区间上使用同一低阶（如二点的，三点的， ）高斯型求积公式算出积分近似值，再相加即将积分 的近似值： ba,mixxii, 2 , 1,1 dxxfba kminkkmhthiafAhG21210

54、 EkktAmabh, 其中， 由查表可得。同时在实际计算中，还常用相邻两次计算结果 和 的关系式 mG1mG111mmmGGG1,mmGG11mG相当于相对误差）即算出 后，观察 是否成立，以判定是否终止计算。请同学们据此编程计算。 来控制运算（当 时，a67证明：以 为节点构造次数 的多项式H(x)使满足条件：由第二章Hermite插值知，H(x)为Hermite插值多项式。由于Gauss公式具有2n+1次代数精度，则它对H(x)能准确成立，即nxxx10,12n), 1 , 0(),( )( ),()(nixfxHxfxHiiii)()()(00knkkknkkbaxfAxHAdxxH0

55、(22)2( )( )()( )( )( )( )( )( )(22)!nbnnaknbbbbaaaaR xf x dxA f xff x dxH x dxf xH x dxx dxn bannnknbadxxnfxfAdxxffR)()22()()()(2)22(0 )()(0nxxxxxxx其中定理：对于Gauss公式（1），其余项为：a68再注意到函数 在 上保号，应用积分中值定理即可得结论。)(2xba, 对比Newton-cotes公式，Gauss不但具有高精度，而且是数值稳 定的。Gauss公式之所以能够保证数值稳定性，是由于其求积函数 具有稳定性。证明:由于nkjjjkjkxxx

56、xxl0)(是n次多项式，因而 是2n次)(2xlk多项式，故Gauss公式对 能准确成立。即)(2xlkniikibakxlAdxxl022)()(nkkkbaxfAdxxf0)()(的求积函数kA全是正的。即), 1 , 0( , 0nkAk定理定理：Gauss求积公式), 1 , 0(nk注意到ikikxlkiik01)(从而上式右端实际上等于kA从而有： ), 1 ,0.(0)(2nkdxxlAbakk故定理得证。a69nkkknxfAI0)(，由于实际问题中，通常不一定提供准确的数据 ，而只是给出含有误差的数据 （如：由于计算机字长的限制而产生的舍入误差）。因而实际求得的积分值为：*

57、kfnokkknfAI*那么，原始数据的误差对积分值的影响能否加以控制呢？ 上面已经证明Gauss公式的求积系数 则 ), 1 ,0.(0nkAkkknknkknkkkknnffAffAII*000*max又由于Gauss公式对f(x)=1准确成立,即 nkkbabaAdxdxxfab01)()(kkxff下面讨论求积过程的数值稳定性问题。利用公式 计算积分时a70则 kknknnffabII*0*max)(由此可以断定Gauss公式是稳定的。，其中 称为权函数。 考察积分badxxfxI)()(0)(x当1)(x时即为普通积分 。 可以仿照处理普通积分的方法讨论带权的积分。如，求积分公式nk

58、nnbaxfAdxxfx0)()()(如果对于任意次数不超过2n+1的多项式均能准确地成立，则称之为Gauss型的。上述Gauss公式求积节点 仍称为Gauss点。同样地， 是Gauss点的充要条件为：kx), 1 , 0(nkxknkkxxx0)()(是ba,区间上关于权函数 的正交多项式)(x若 1, 1ba.且权函数211)(xx则所建立的Gauss求积公 20( )()*1nbkkakf xdxA f xx称之为（高斯-切比雪夫公式）式为：a71 确定了Gauss点后，再利用Gauss公式2n+1次代数精确度关系定 值得指出的是值得指出的是：运用正交多项式的零点构造Gauss求积公式，

59、这种方法只针对某些特殊的权函数才有效。 构造Gauss公式的一般方法是利用代数精度的概念用待定系数法求解。 现举一例加以说明： 设要构造下列形式的Gauss公式：kA)()()(110010 xfAxfAdxxfx令它对32, 1)(xxxxf准确成立，得：9272,52,32131030121020110010AxAxAxAxAxAxAA211x的正交多项式是切比雪 夫多项式，因此求积公式（*）的Gauss点是n+1次切比雪夫多项式的 零点，即为： ), 1 ,0(),2212cos(nkkkxk 由于区间-1，1上关于权函数a72从而所构造的Gauss求积公式为：).289949. 0(2

60、77556. 0)821162. 0(389111. 0)(10ffdxxfx 例1：选取常数a，求使积分公式)()0()()0(2)(20hffahhffhdxxfh的代数精度尽量高，并问其代数精度为几次？解：取f(x)=1.则对上求积公式，左端右端。112hh2( ),( )2f xxfxx解此方程组得：01010.8261162,0.2899490.389111,0.277556xxAAf(x)=x. 则 ，左端= ( )1fx 212h 右 端a733322322202,31ahhhahhhh右端左端令左端=右端121 a233)(.)(xxfxxf又取左端=4223441)30(12