第7章 数值积分与数值微分

第七章数值积分与数值微分积分问题最早来自于几何形体的面积、体积计算，也是经典力学中的重要问题（例如计算物体的重心位置）. 在现实应用中，很多积分的结果并不能写成解析表达式，因此需要通过数值方法来计算. 数值微分是利用一些离散点上的函数值近似计算某一点处的函数导数，它针对表达式未知的函数. 本章介绍一元函数积分（一重积分）和微分的各种数值算法，它们也是数值求解积分方程、微分方程的基础.7.1数值积分概论7.1.1基本思想考虑如下定积分的计算：I(f)≡∫f(x)dxba,(7.1) 其中函数f: ℝ→ℝ，首先应想到的是微积分中学习过的牛顿-莱布尼兹（Newton-Leibniz）公式：∫f(x)dxba=F(b)−F(a) ,其中F′(x)=f(x),即F(x)为f(x)的原函数. 但是，诸如e x2,sinxx,sinx2等表达式很简单的函数却找不到用初等函数表示的原函数，因此必须研究数值方法来近似计算积分. 另一方面，某些函数的原函数虽然可以解析表示，但其推导、计算非常复杂，此时也需要使用数值积分方法.一般考虑连续的、或在区间[a,b]上可积①的函数f(x)，则根据积分的定义有：lim n→∞,ℎ→0∑(x i+1−x i)f(ξi)ni=0=I(f) , (7.2)其中a=x0<x1<⋯<x n+1=b,ξi∈[x i,x i+1],i=0,⋯,n, ℎ=max0≤i≤n(x i+1−x i).公式(7.2)实际上也反映了近似计算积分的思路，就是取充分大的n, 用函数值的加权和逼近准确的积分值. 研究数值积分方法主要是探讨如何用相对较少的计算成本得到准确度较高的结果，这里的成本用计算被积函数值的次数来衡量.上述讨论表明，近似计算积分I(f)的数值积分方法（numerical quadrature）一般具有如下形式：I n(f)≡∑A k f(x k)nk=0, (7.3)其中a≤x0<x1<⋯<x n≤b, A k,(0≤k≤n)为一组系数. 形如公式(7.3)的求积公式称为机械求积公式，其中系数A k,(0≤k≤n)称为积分系数，自变量取值x k,(0≤k≤n)称为积分节点. 根据积分节点和积分系数的不同设置，可得到各种具体的求积公式. 针对实际的数值积分问题，往往还需使用多种求积公式构造算法，以达到满意的效果.推导求积公式的一种方法是用多项式函数p(x)来近似f(x)，则可期望有以下的近似关系：①连续函数在闭区间内一定有界、可积，而可积函数则可能不连续.∫f (x )dx b a ≈∫p (x )dx ba ,其中多项式函数的积分很容易通过牛顿-莱布尼兹公式求出. 假设使用拉格朗日插值法构造p (x )，区间[a, b]内的插值节点为x 0,x 1,…,x n ，则p (x )=L n (x )=∑f (x k )l k (x )nk=0,l k (x)为拉格朗日插值基函数. 由此得到求积公式为：I n (f )=∫∑f (x k )l k (x )n k=0dx b a =∑f (x k )∫l k (x )dx b an k=0 .(7.4) 由于l k (x )为拉格朗日插值基函数，一旦插值节点确定，积分∫l k (x )dx b a 可方便地计算出来. 这种用多项式插值近似被积函数得到的求积公式(7.4)被称为插值型求积公式（interpolatory quadrature ），易知它也是一种机械求积公式，其积分节点就是插值节点，而积分系数 A k =∫l k (x )dx b a ,(k =0,1,⋯,n ).(7.5)下面的例子推导n=0, n=1两种情况下的插值型求积公式，它们分别称为中矩形公式（midpoint rule ）和梯形公式（trapezoid rule ）.例7.1 (中矩形公式与梯形公式)：根据n=0, n=1两种情况对应的拉格朗日插值推导相应的求积公式，假设插值节点分别为区间[a, b]的中点和两个端点.[解] 当n=0时，按题意设x 0=(a +b)/2, 由于0次拉格朗日插值多项式为常数，则L 0(x )=f (x 0)因此，I 0(f )=∫f (x 0)dx b a =(b −a )f (a +b 2) . (7.6)当n=1时，按题意设x 0=a , x 1=b , 利用线性拉格朗日插值基函数和公式(7.5)，求出：A 0=∫l 0(x )dx b a =∫x −b dx b a=b −a , A 1=∫l 1(x )dx b a =∫x −a b −a dx b a =b −a 2. 因此，I 1(f )=∑A k f (x k )1i=0=b −a 2[f (a )+f(b)] . (7.7)中矩形公式(7.6)和梯形公式(7.7)具有很直观的几何意义，即分别用矩形面积和梯形面积来近似函数曲线和横轴围成区域的面积（如图7-1）.7.1.2求积公式的积分余项与代数精度定义7.1：对于计算积分I(f)的求积公式I n (f)，称I(f)−I n (f)为该公式的积分余项，常记为R [f ].积分余项反映了求积公式的截断误差，是衡量求积公式准确度的重要依据. 假设I n (f)为某个插值函数p(x)的积分，则R [f ]=∫[f (x )−p (x )]dx ba,即积分余项等于插值余项的积分. 对于插值型求积公式(7.4)，有R [f ]=∫[f (x )−L n (x )]dx b a =∫f (n+1)(ξ)(n +1)!ωn+1(x )dx b a, (7.8) 其中ξ依赖于x .下面介绍代数精度的概念，它是衡量求积公式准确度的另一个重要标准.定义7.2：如果某求积公式对于次数不超过m 的多项式均准确成立，但对于m+1次多项式可能不准确，则称该求积公式具有m 次代数精度（degree of exactness ）.上述定义表明，一个求积公式具有较高次的代数精度，就意味着它能准确计算次数较高的多项式的积分②. 应注意，在某些情况下代数精度并不是越高越好.要判断一个机械求积公式的代数精度，最直接的方法是考察当被积函数分别为1,x,…,x m 时求积公式的准确性. 下面给出一个定理，其证明留给感兴趣的读者思考. 定理7.1：机械求积公式I n (f )=∑A k f (x k )n k=0至少有m 次代数精度的充要条件是当f (x )分别为1,x,…,x m 时，I (f )=I n (f ) .我们讨论的所有求积公式都至少具有0次代数精度，因此根据定理7.1它们应对f (x )=1的积分准确，则推出：∑A k n k=0=∫1dx ba =b −a ,(7.9)这说明积分系数之和等于区间长度.根据插值型求积公式的含义，容易得出如下定理，其证明留给读者思考.定理7.2：机械求积公式I n (f )=∑A k f (x k )n k=0是插值型求积公式(7.4)的充要条件是它至少有n 次代数精度.考察例7.1中的中矩阵公式和梯形公式，可得出它们都具有1次代数精度. 一般地，在部分积分节点、积分系数已知的情况下，利用定理7.1可建立方程求解剩余的积分系数或节点，使其达到一定的代数精度. 而且定理7.1中求积公式的形式还可以更一般，只要是函数值或其导数的线性组合即可，下面的例子说明了这种情况.例7.2(求积公式的代数精度)：用形如H 2(f )=A 0f (0)+A 1f (1)+B 0f ′(0)的求积公式近似积分I (f )=∫f(x)dx 10，试确定系数A 0, A 1, B 0，使公式具有尽可能高的代数精确度.[解] 根据题意可令f (x )=1，x ，x 2分别代入求积公式使H 2(f )=I(f)精确成立.② 当然，这没有什么实际意义，因为很容易得到多项式函数的原函数.(a) 中矩形公式(b) 梯形公式图7-1 中矩形公式(a)和梯形公式(b)的示意图.当f (x )=1时，得A 0+A 1=∫1∙dx 1=1当f (x )=x 时，得A 1+B 0=∫xdx 1=12 当f (x )=x 2时，得A 1=∫x 2dx =1310联立上述三个方程，解得A 1=13,A 0=23,B 0=16 .当f (x )=x 3时，容易验证上述求积公式不准确，因此H 2(f )最多具有2次代数精度. 7.1.3求积公式的收敛性与稳定性实际使用的求积公式都是机械求积公式(7.3)，下面针对它给出求积公式的收敛性和稳定性的概念.定义7.3：对于n 的值可以任意的一系列机械求积公式I n (f)=Σk=0n A k f (x k )，a ≤x 0<x 1<⋯<x n ≤b ，若lim n→∞,ℎ→0I n (f)=∫f (x )dx ba ,其中ℎ=max 1≤k≤n (x k −x k−1)，则称这一系列求积公式具有收敛性.收敛性说明求积公式在积分节点逐渐增多、且节点间距逐渐变小时，其结果收敛到准确的积分值. 这个概念不同于公式(7.2)，后者反映的是被积函数具有可积性. 在实际应用中，求积公式具有收敛性非常重要，后面还将针对具体的公式加以讨论.在讨论求积公式的稳定性之前，先分析数值积分问题的敏感性和条件数. 假设f(x)为准确的被积函数， f̃(x )为实际计算时受扰动影响的被积函数，扰动的大小为δ=‖f (x )−f ̃(x )‖∞=max a≤x≤b |f (x )−f ̃(x )|，则扰动对积分计算的影响为： |∫f (x )dx b a −∫f̃(x )dx b a |≤∫|f (x )−f ̃(x )|dx b a ≤(b −a )δ, (7.10)这说明，积分计算结果的误差最多为扰动的(b −a)倍，积分区间的长度(b −a)是绝对条件数的上限. 一般来说，数值积分问题是不太敏感的. 这一点不难理解，因为积分运算本身就是一个平均的过程，它不易受被积函数的某些微小变化的影响.求积公式的稳定性反映计算过程中的扰动是否被放大、以及放大的程度. 具体来说，在计算机械求积公式时，需考虑积分节点的函数值出现误差时，它对结果产生的影响. 假设节点函数值由f (x k )变为f̃(x k )，则数值积分的结果由I n (f )变为I n (f ̃)，两者之差满足： |I n (f )−I n (f̃)|=|∑A k [f (x k )−f ̃(x k )]n k=0|≤∑|A k |nk=0∙|f (x k )−f ̃(x k )|≤(∑|A k |n k=0)ε,(7.11)其中ε=max 0≤k≤n|f (x k )−f ̃(x k )|≤δ. 根据(7.9)式，若同时有A k >0,k =0,…,n ，则不等式(7.11)变为：|I n (f )−I n (f̃)|≤(b −a )ε≤(b −a )δ . (7.12) (7.12)式表明求积公式的结果受扰动影响的程度与积分问题敏感性的结果(7.10)式一致，这是控制数值计算误差能达到的最佳情况. 将它作为一个标准，可定义求积公式的稳定性.定义7.4：若对k =0,1,…,n ，均有A k >0, 则机械求积公式I n (f)=Σk=0n A k f (x k )是稳定的.利用定义7.4很容易直接判断求积公式的稳定性. 在实际情况中，不稳定的求积公式的积分系数绝对值之和Σk=0n |A k |可能远大于b −a ，从而导致函数值的扰动在计算结果上被放大很多.本节介绍了求积公式的基本形式，以及积分余项、代数精度、收敛性和稳定性的概念，其中收敛性是针对一系列公式（积分节点数目逐渐增多）而言的. 后面介绍具体公式时，将考察单个公式的积分余项、代数精度和稳定性，并讨论积分节点数目逐渐增多时的收敛性. 此外，还应注意具体公式中计算函数值的次数，它是度量计算量大小的标准.7.2牛顿-柯特斯公式在积分区间上构造等距节点的多项式插值，对应的插值型求积公式为牛顿-柯特斯公式.7.2.1 柯特斯系数与几个低阶公式假设将积分区间n 等分，步长ℎ=(b −a)/n ，插值节点为x k =a +kℎ,(k =0,…,n)，则可得到等距节点的拉格朗日插值多项式. 根据公式(7.4), (7.5)的推导，得到型如I n (f )=∑A k f (x k )n k=0的求积公式，其中A k =∫l k (x )dx b a =∫(x −x 0)⋯(x −x k−1)(x −x k+1)⋯(x −x n )(x k −x 0)⋯(x k −x k−1)(x k −x k+1 )⋯(x k −x n )dx b a. 这就是n 阶牛顿-柯特斯（Newton-Cotes ）公式.引入变量代换：x =a +tℎ,t ∈[0,n], 则A k =∫∏(t −j )n j=0,j≠k ℎdt n0=∫∏(t −j )n j=0,j≠k b −a dt n 0 . (7.13) 令C k (n )=1n ∫∏(t −j k −j )n j=0,j≠k dt n 0 ,(k =0,⋯,n ), (7.14)它仅与阶数n 有关，而与区间大小无关，则积分系数A k =(b −a )C k (n ).(7.15) 公式(7.14)中的C k (n )常被称为柯特斯系数，可以预先计算出不同n 值对应的柯特斯系数，制成一个表（如表7-1所示），根据它方便地写出各阶牛顿-柯特斯公式.表7-1 柯特斯系数表[本章知识点]: 机械求积公式；积分余项；代数精度；插值型求积公式及其代数精度；中矩形公式、梯形公式；求积公式的稳定性；求积公式的收敛性；牛顿-柯特斯公式及其代数精度；柯特斯系数；辛普森公式、柯特斯公式；梯形公式、辛普森公式的积分余项；复合梯形公式；等距节点求积公式的准确度阶数；复合辛普森公式；步长折半情况下的复合求积公式计算；复合梯形公式的余项展开式；理查森外推法；龙贝格积分算法；自适应积分算法的原理；高斯求积公式；带权积分的高斯求积公式；高斯点的性质以及高斯求积公式的构造；高斯-勒让德公式；向前、向后、中心差分公式；二阶中心差分公式；有限差分公式的准确度阶数；插值型求导公式；数值微分的外推算法.算法背后的历史：“数学王子”高斯卡尔·弗里德里希·高斯（Johann Carl Friedrich Gauss，1777年4月30日—1855年2月23日）是德国著名数学家、物理学家、天文学家、大地测量学家. 高斯的成就遍及数学的各个领域，在数论、非欧几何、微分几何、超几何级数、复变函数论以及椭圆函数论等方面均有开创性贡献. 高斯被公认为是十九世纪最伟大的数学家，与阿基米德、牛顿齐名.高斯于1777年4月30日生于布伦瑞克的一个工匠家庭. 1792年，15岁的高斯进入布伦瑞克学院，在那里独立发现了二项式定理的一般形式、数论上的“二次互反律”、“质数分布定理”和“算术几何平均”. 1795年高斯进入哥廷根大学. 1796年，19岁的高斯取得了一个数学史上极重要的结果，就是论文“正十七边形尺规作图之理论与方法”. 1801年，高斯又证明了"Fermat素数"边数的正多边形可以由尺规法作出. 高斯从1807年起担任哥廷根大学教授兼哥廷根天文台台长直至逝世.高斯的主要贡献●高斯消去法●质数分布定理和最小二乘法●高斯分布（标准正态分布）●高斯求积公式●无穷级数的高斯判敛法●发现了谷神星的运行轨迹●证明代数基本定理●微分几何最重要的创始人高斯的故事一天，德国哥廷根大学，一个19岁的青年吃完晚饭，开始做导师单独布置给他的数学题. 正常情况下，他总是在两个小时内完成这项作业. 像往常一样，前两道题目在两个小时内顺利地完成了. 第三道题写在一张小纸条上，是要求只用圆规和一把没有刻度的直尺做出正17边形，他没多想，埋头做起来. 然而，做着做着，他感到越来越吃力. 困难激起了他的斗志：我一定要把它做出来！天亮时，他终于做出了这道难题. 导师看了他的作业后惊呆了. 他用颤抖的声音对青年说：“这真是你自己做出来的？你知不知道，你解开了一道有两千多年历史的数学悬案？阿基米、牛顿都没有解出来，你竟然一个晚上就解出来了！……我最近正在研究这道难题，昨天不小心把写有这个题目的小纸条夹在了给你的题目里. ”多年以后，这个青年回忆起这一幕时，总是说：“如果有人告诉我，这是一道有两千多年历史的数学难题，我不可能在一个晚上解决它. ”这个青年就是数学王子高斯.高斯在1810年提出了消去法，他的这个方法是为了简化二次型而不是为了矩阵分解. 实际上，求解线性方程组中的高斯消去法是由20世纪40年代由Dwyer 、冯 偌伊曼（von Neumann ）等人提出的. 高斯最早提出，通过积分点的最佳分布可得到更为精确的数值积分方法，并于1814年发表了第一个这样的求积公式，高斯求积公式因此得名.高斯名言● 数学是科学里的皇后，数论是数学中的女王.● 宁可少些，但要好些.● 数学中的一些美丽定理具有这样的特性: 它们极易从事实中归纳出来，但证明却隐藏得极深.● 给我最大快乐的，不是已懂得的知识，而是不断的学习；不是已有的东西，而是不断的获取；不是已达到的高度，而是继续不断的攀登.练习题1. 确定下列求积公式中积分系数或积分节点的待定值，使其代数精度尽量高，并指明所构造的求积公式所具有的代数精度：(1)∫f (x )dx ≈ A −1f (−ℎ)+ A 0f (0)+ A 1f(ℎ)2ℎ−2ℎ(2)∫f (x )dx ≈ 1−1[f (−1)+2f (x 1)+3f(x 2)]3⁄2. 若积分节点x k ,k =0,…,n 给定，要求机械求积公式的积分系数，使得求积公式至少具有n 次代数精度. 请根据定理7.1列出待求解的线性方程组，并判断解的存在性和唯一性.3. 直接验证柯特斯公式：C =b −a [7f (x 0)+32f (x 1)+12f (x 2)+32f (x 3)+7f (x 4)] 具有5次代数精度. 4. 用辛普森公式求积分∫e −x dx 10并估计误差.5. 证明下列等式，它们分别说明了三种矩形求积公式及其余项公式.(1).∫f (x )dx =(b −a )f (a )+ f ′(η)2(b −a)2b a,η∈(a,b); (2).∫f (x )dx =(b −a )f (b )− f ′(η)2(b −a)2b a,η∈(a,b) ;(3).∫f (x )dx =(b −a )f (a +b )+ f ′′(η)(b −a)3b a,η∈(a,b) . 6. 对下列积分，分别用复合梯形公式和复合辛普森公式计算, 其中n 表示计算中使用n +1个区间等分点上的函数值，然后比较两种方法计算结果的准确度.(1).∫x 4+x 2dx,n =810(2).∫√xdx,n =491。