数值计算-第7章数值微分和数值积分

1、第7章数值微分和数值积分7.1数值微分7.1.1差商与数值微分当函数是以离散点列给出时，当函数的表达式过于复杂时，常用数值微分近似计算畑的导数广（初。在微积分中，导数表示函数在某点上的瞬时变化率，它是平均变化率的极限；在几何上可解释为曲线的斜率；在物理上可解释为物体变化的速率。以下是导数（力的三种定义形式：广=1摆公牛空KJhs2h（7.1）在微积分中，用差商的极限定义导数；在数值计算中返璞归真，导数取用差商（平均变化率）作为其近似值。最简单的计算数值微分的方法是用函数的差商近似函数的导数，即取极限的近似值。下面是与式（7.1）相应的三种差商形式的数值微分公式以及相应的截断误差。向前差商用向前

2、差商（平均变化率）近似导数有：h（7.2）其中+k的位置在坯的前面，因此称为向前差商。同理可得向后差商、中心差商的定义。由泰勒展开f（x0+h）=+h得向前差商的截断误差：向后差商用向后差商近似导数有：沟（7.3）与计算向前差商的方法类似，由泰勒展开得向后差商的截断误差盹）寸-皿匕（）=0h）中心差商中心差商用中心差商（平均变化率）近似导数有7.4)7.4)由泰勒展开由泰勒展开了（观+血）=了（珥）+幼认）+务广（观）+令厂得中心差商的截断误差：应（对=/-（咼）-应（对=/-（咼）-川DFD备八+ZU）2h差商的几何意义微积分中的极限定义，表示了微积分中的极限定义，表示了在“叼处切线的斜畑+

3、却-畑+却-了）率，即图7.1中直线尸的斜率；差商h表示过九E和（观+曲/筑+划）两点直线的斜率，是一条过珥的割线。可见数值微分是用近似值内接弦的斜率代替准确值切线的斜率。图7.1微商与差商示意图例7.1给出下列数据，计算艸叫艸网艸呵，艸网,0.020.040.060.80.105.065.075.0655.055.0550.020.040.060.80.105.065.075.0655.055.055/YQ02)刘/X0.06)-/解：八丿(5.075.06)/(0.040.02)=0.5/X0.06)-/X0.10)-r(o.os)-5.055.07)/(0.080.04)=0.5/X0.

4、10)-r(o.os)-5.055.055)/(0.080.10)=0.25(/(0.10)(0.06)/(0.100.06)=18.75:厂L宀戸.上土1丄乙设定最佳步长在计算数值导数时，它的误差由截断误差和舍入差两部分组成。用差商或插值公式近似导数产生截断误差，由原始值X的数值近似产生舍入误差。在差商计算中，从截断误差的逼近值的角度看，I糾越小，则误差也越小；但是太小的I糾会带来较大的舍入误差。怎样选择最佳步长，使截断误差与舍入误差之和最小呢？一般对计算导数的近似公式进行分析可得到误差的表示式，以中心差商为例，截断误差不超过而舍入误差可用量曲估计（证明略），其中曰是函数朋的原始值的绝对误差

5、限，总误差为当当时，总误差达到最小值，即*)可以看到用误差的表达式确定步长，难度较大，难以实际操作。通常用事后估计方法选取步长用,例如，记D醐D|通常用事后估计方法选取步长用,例如，记D醐D|为步长等于h2的差商计算公式，给定误差界当时，2就是合适的步长。例7.2对函数，取不同的步长计算，观察误差变化规律，从而确定最佳步长。解：误差误差0.103.16300.00480.053.15900.00080.093.1622-0.00400.043.15880.00060.083.16130.00310.033.15830.00010.073.16070.00250.023.15750.00070.

6、063.16000.00180.013.15500.0032表中数据显示，当步长应从0.10减少到0.03时，数值微分误差的绝对值从0.0048减少到0.0001，而随着h的进一步减少，误差的绝对值又有所反弹，表明当步长曲小于0.03时,舍入误差起了主要作用。在实际计算中是无法得到误差的准确数值的，这时以最小为标准确在实际计算中是无法得到误差的准确数值的，这时以最小为标准确定步长，本例中取曲=0.04。7.1.2插值型数值微分对于给定的畑的函数表，建立插值函数皿初，用插值函数圮对的导数近似函数冷的导数。设H用为小上的节点，给定毎弘）厂0丄7，以尽2）为插值点构造插值多项式，以g的各阶导数近似m

7、的相应阶的导数，即了（x）=厶（兀）=工（小/（阳）2-0八兀）=4。）=S?；（x）/（珂）2-0时，f%）f%）疼总&于召）=0,1,卫3-07.5）误差项为：例7.3给定Kgw例7.3给定Kgw，并有宀rf，并有，计算解：作过爼畑）厂0的插值多项式:将吗代入，（方得三点端点公式和三点中点公式:八坯）=（-”（观）+4了（珂）-了区）广（珂）=存-血）+了区）f显）=舟（孑（观）-4孑凶）+（吃）利用泰勒（Taylor）展开进行比较和分析，可得三点公式的截断误差是。曲。类似地，可得到五点中点公式和五点端点公式：八坯）=吕了-泌）-卩临-血）+8/（x0+紛-畑+M）12/2兵观-2血氏+2

8、h/Vo）=总L/%）+4盯（观+）-36/（x0+2h）+1叮0+弘）-（坯+恥）+牛严奘砂。+4划7.1.3样条插值数值微分把离散点按大小排列成以从57，用粋关系式构造插值点(召，于(召，于(吗)=0,1,2,，挖的样条函数纨力:1+2忑-再+i1阳+1-聲丿帆-吗+12$.,X-xi+l必+伏一吗）m.1再一召+1f/+1+2一吗+iL聲-聲+1和-“疋兀风+1+0一聲+1）强+i严+1_%则八时，可用计算导数。则八时，可用计算导数。7.2数值积分在微积分中用牛顿一莱布尼兹(Newton-Leibniz)公式计算连续函数的定积分:ff(x)dx=F(b)-F(a)但是，当被积函数是以点列

9、用的形式给出时，当被积函数的原函数月(力难以得到时，例如止,则无法用牛顿一莱布尼兹积分公式计算。有时当被积函数的原函数过于复杂时，也不宜套用积分公式计算积分，而应采用数值积分公式计算定积分。在微积分中，定积分是黎曼(Rimann)和的极限，它是分割小区间长度趋于零时的极限，即在数值积分公式中，只能用有限项的和近似上面的极限，通常由函数在离散点函数值的线性组合形式给出。KJ）=（召）记7，在本章中，用心）表示精确积分值，用2表示近似积分值，区称为积分节点，陽称为积分系数。确定2中积分系数陽的过程就是构造数值积分公式的过程。怎样判断数值积分公式的效果？代数精度是衡量数值积分公式优劣的重要标准之一。

10、定义7.1（代数精度）记也旬上以卫为积分节点的数值积分公式若2满足恥wexn丄期而，则称a具有嗨阶代数精度。由此可知当2具有禅阶代数精度时，对任意的卿阶多项式都有wy。7.2.1插值型数值积分对给定的被积函数在也旬上的点列（召，于（吗）=0,1,2,，挖作拉格朗日插值多项式，以近似计算如，则有对给定的被积函数在也旬上的点列（召，于（吗）=0,1,2,，挖作拉格朗日插值多项式，以近似计算如，则有数值积分误差，也就是对插值误差的积分值或&C/）=f于坯血,心-吗）必对一般的函数瓦H,但若冷是一个不高于挖次的多项式，由于严7而有爲（了）=。因此，阶插值多项式型式的数值积分公式至少有阶代数精度。例7.

11、4建立闵上节点为右2的数值积分公式。碍=（畫）必解：由得得以数值积分公式：C/)=t-3/(0)+16/(0.5)+5/(2)722牛顿-柯特斯（Newton-Cotes）积分把积分区间也切分成挖等分，记步长为挖,取等分点盹丸丄卫）作为数值积分节点，构造拉格朗日插值多项式，取由此得到的数值积分称为牛顿柯特斯积分。下面可以看到，牛顿柯特斯积分系数和积分节点以及积分区间无直接关系，系数固定而易于计算。梯形积分以e畑和了）为插值节点构造线性函数“，有f（x）dx瘗f厶（对必那么，厶（对必=f仏（x）f（心）+k（x）f（珂）必I*(x)dx=丄dx=(b-a)=-a)中提取公因子少一后，得到牛顿一柯

12、特斯积分的组合系数：2，2，它们已与积分区间没有任何关系了。了（加辛了+八毋T（J）=字了）+了记上（7.6）称乳刀为梯形积分公式。它的几何意义是用梯形面积近似代替积分值（图7.2）。图7.2梯形积分怎样确定梯形积分公式的代数精度？我们可以取)=1,验证取畑时，有(/)=(dx=baF=乎C/0)+/乎(1+1)=占p即：取/W=X时，有即：取/W=X时，有丁=乎C/S)+了)戶(+坊=与1即：取2时，有/（/）=f急=兰子H乎C/+/）=Tf得梯形求积公式具有一阶代数精度。梯形积分公式的误差：由Eg=f冷#(x-L3)(xd)必得因为匕-泳宀在禺旬上不变号，由积分中值定理得到梯形求积公式的截

13、断误差:7.7)辛普森（Simpson）积分对区间讪作二等分，记廿7=3对区间讪作二等分，记廿7=3卸H。以仮）,+切2+卯2）和为插值节点构造二次插值函数g），那么，有=J(x-(g+b)/2)(x-b)(a-(a+h)/2)(ah)=+&-总)=0-总)即=竿卜）+4彳|+伯（7.8）Sf）称为辛普森或抛物线积分公式。它的几何意义是用过三点的抛物线面积近似代替积分的曲边面积（图7.3）。图7.3抛物线积分面积分别将=代入到/和占中，可以得到心=辛1+4/畔+子=心）表明辛普森公式对于次数不超过三次的多项式准确成立，应力具有三阶代数精度。因此可设一个三次多项式满足条件：a+ba+b计算得到误

14、差为：=x（xhab计算得到误差为：=x（xha0定义中的挖Tw的。max(无ja-xJ0包含了1，通常都要求计算积分的求积公式是收敛稳定性是研究计算公式当曲有误差血时,W的误差是否增长，现设心叮，误差记为斜张mi)定义7.3定义7.3对任给小，只要汨皿r，就有14(/)厶匚=工昭(吗)则称求积公式是稳定的。定理7.1若求积公式的系数碍叫=0,1,，则求积公式是稳定的。证明由于详gTW证明由于详gTW，故有14-4(7)1=1习辺C/E)-Q卜占工%=呢弋)Vs0,3=-，就有于是对，就有故求积公式是稳定的。7.3复化数值积分由插值的龙格现象可知，高阶牛顿-柯特斯积分不能保证等距数值积分系列的

15、收敛性，同时可证（略）高阶牛顿-柯特斯积分的计算是不稳定的。因此，实际计算中常用低阶复化梯形等积分公式。7.3.1复化梯形积分把积分区间分割成若干小区间，在每个小区间V川上用梯形积分公式，再将这些小区间上的数值积分累加起来，称为复化梯形公式。复化梯形公式用若干个小梯形面积逼近积ff（x）dx分血比用一个大梯形公式效果显然更好，如图7.4所示。这种作法使我们想起定积分定义，即它为被积函数无限分割的代数和。这也正是计算定积分最朴素的算法。图7.4复化梯形公式积分视图复化梯形积分计算公式切作等距分割，有，于是了必唔广了（电在g上门（曲冷心+小）-心*上，则有M-l心)=jjjM-l”+号E+”叫-三

16、广）寻12记挖等分的复化梯形公式为爼或叽旳，有T(h)=TK(J)=h”+他+询+”)复化梯形公式截断误差7.11)，根据均值定理，当时，存在兵“，有，于是7.12)丄由此看到复化梯形公式的截断误差按照妒或者屏的速度下降，事实上，可以证明，只要畑在d上有界并黎曼可积，当分点无限增多时，复化梯形公式收敛到积分!(./)=f(x)dx，则有，则有对于任给的误差控制小量0，有M2s式中M20，只要设了亡/向，在恋，吩设了亡/向，在恋，吩+2】上的误差为设了亡/向，在恋，吩设了亡/向，在恋，吩+2】上的误差为就有遛，计算中要求有5位有效数字。用复化梯形和复化辛普森求积心，计算中要求有5位有效数字。用复

17、化梯形和复化辛普森求积例7.5求公式的分点应取多少？解.了7/0)=严0)7irwn/C4)wiox0，要，只需W）-T5%即可，而用鬼一兀I作为控制手段简单直接，序列在计算机上也不难实现。复化积分的算法描述从数值积分的误差公式可以看到，截断误差随分点的增长而减少，控制计算的精度也就是确定分点数。在计算中不用数值积分的误差公式确定分点数挖的理论模式，而用I爲-石（刀|eT1=T2M-1禺弋工乳）H=Hn计算新增节点的值H=HnT2=(Tl+H)/2H=h/2，n=2n/将区间一分为二endwhile4输出积分值T2。在自动控制误差算法中初始分点值不宜过小，以防假收敛。734龙贝格（Romber

18、g）积分爲-抵C/）T（八由前面得到的关系式（7.17），可以将5作为$的修正值补充到心），7.18)7.18)其结果是将梯形求积公式组合成辛普森求积公式，截断误差由曲提高到（沪）。这种手段称为外推算法。外推算法在不增加计算量的前题下提高了误差的精度，是计算方法中一种常用手法。不妨对比心S再做一次组合。皿=皿=-2880(ba)厶7C/）一普P）3讣得到心）-砥诂毘心）-凡）心孕-存C/)7C/)7.19)心孕-存C/)7C/)7.19)复化辛普森公式组成复化柯特斯公式，其截断误差是曲。同理对柯特斯公式进行组合：心）7心）=住得到具有7次代数精度和截断误差是曲的龙贝格公式:尽=|%刀-右GC/

19、）还可以继续对W做上去。为了便于在计算机上实现龙贝格算法，将九昭统一用血h表示，列标g分别表示梯形、辛普森、柯特斯积分，行标疋表示分点数存或步长h/21j。龙贝格计算公式：对每一个匕从2做到疋，一直做到I鬆厂鬆31小于给定控制精度停止计算。龙贝格算法龙贝格算法按表7.2元素的行序进行运算，知息“氐在计算中每个元素只用到上一行和本行的元素。对上面的算法进一步优化，对每k行可将计算定义在两行元素之间,令；在每计算一行元素后，要将表7.2表7.21.输入区间端点恥，精度控制值，循环次数虫，定义函数，取23.fork=2toM瓦=(兔祖+血_1工了3+-W)/2/hk=初沪forJ=2tok傀厂堆冲+

20、（限-鬆心）“铲】-1）泊隔-也-心退岀循坏4.输出也。7.4重积分计算在微积分中计算二重积分是用化为累次积分的方法进行的。计算二重数值积分也是计算累次数值积分的过程。为了简化问题，我们仅讨论矩形域上的二重积分。有很多非矩形域上的二重积分可作变换将其转换到矩形域上。7.20)7.20)其中：是常数，了Z在。上连续。像在微积分中一样，将二重积分化为累次积分：7.21)二重积分的复化梯形公式对区间“切和叩分别选取正整数询，在兀轴和厂轴上分别有步h=匕衣=匕m旳用复经梯形公式计算EE，计算中将当作常数，有M-lJfy)dy花M-lJfy)dy花（“）+訂（兀必）+工了（乜）7.22)再将当作常数，在

21、*方向上计算式（7.23）中每一项的积分，有必厂:抑0必）必二M-1|”（和曲必厂:抑0必）必二M-1|”（和曲+”（耳,曲+工/屛兀）M-1I”区必）+”（耳必）+&（西必）则fM-l工乳再Jo）M-l-1工血必j57（心丹）+M-1工刃兀必）j-1J-lJ2-0J-0积分区域的4个角点的系数是1/4，4个边界的系数是1/2，内部节点的系数是1。误差：妤詁（卩叭+F等丘）和窗在积分区间内。例子7.6用复化梯形公式计算重积八和窗在积分区间内。例子7.6用复化梯形公式计算重积八ffsin（y+刃矽必取血=k=0.25解:了Z）如表7.3所示:表7.3了匕旳数值表X1.001.251.501.75

22、2.000.000.8414710.9489850.9974950.9839860.9092970.250.8735750.9668270.9999660.9709320.881530.500.9489850.9974950.9839860.9092970.7780730.750.9999660.9709320.331530.7373190.5472651.000.9092970.7780730.5984720.38110.14112的数值如表7.4所示:表7.4%数值表X0123401/41/21/21/21/411/21111门21/21111/231/21111/241/41/21/21

23、/21/4+y)dydx=0.25x0.25|(0.841471+0.909297+0.909297+0.14112)+1(0.948985+0.997495+0.983986)+1(0.778073+0.598472+0.381661)+j(0.873575+0.948985+0.999966)+10.88153+0.778073+0.547265)+(0.966827+0.999966+0.970932+0.997495+0.983986+0.909297+0.970932+0.88153+0.737319)=0.873601的准确值是0.886176。二重复化辛普森求积公式对区间和匕引分

24、别欣等分和挖等分,在对区间和匕引分别欣等分和挖等分,在轴和，轴上分别有步长均为偶数均为偶数类似于二重复化梯形公式推导，得口=他场，如=1,424,24,1误差：和窗在积分区间内。按例7.6的化分区间，叫的值如表7.5所示:表7.5*012340:2:75高斯（Gauss）型积分公式介绍对插值型积分公式在牛顿-柯特斯积分公式中要求节点是等距的，其优点是计算积分系数的公式规则相同，缺点是制约了求积公式的代数精度。可以证明：当节点个数挖为偶数时，求积公式具有挖T阶的代数精度；当节点个数为奇数时，求积公式具有阶的代数精度。如果我们不预先指定求积节点阳的位置，吗和权系数碍都作为待定的常数，能否适当地确定

25、它们，以提高积分公式的代数精度？回答是肯定的。加个待定常参数，需要加个方程来确定，取一个函数组：兀，,“门，这一组函数构成了抵T次多项式的基，任一小于等于MT次的多项式，都可以用这组函数的线性组合来表示。如果某一积分公式，对这组函数都能精确积分，则此积分公式就有次代数精度。例7.7计算求积系数例7.7计算求积系数】心和求积节点心兀至少具有3阶代数精度。解：按照求积公式的代数精度定义，分别令=，得方程组:11+1=fdx=2解方程组得：求积公式：JJ（x伽=7（-0.57735）+/（0.57735）按例7.7的方式，构造更高阶的代数精度的求积公式，生成求积系数和求积节点的方程组并无困难，而求解

26、该方程组则无一定的章法可循。一般地，通过计算正交多项式的零点作为求积节点。当取积分节点为正交多项式的零点时，则节点个数是挖的求积公式具有一1阶的代数精度。并称积分节点为正交多项式的零点的数值积分公式为高斯型积分公式。为了一般性，考虑积分其中盼）（叭小0）称为权函数。当炉时，即是普通的积分。对于不同的权函数琢选定的节点也不相同。如何构造高斯型积分公式呢？对给定的及权函数防），由施密特（Schmidt）正交化过程作出正交多项式爲（X）尼（X）,丘（X）；解出正交多项式恥）的挖个零点W%，这个挖个零点就是积分节点；以这些节点构造插值多项式，计算积分系数碍=丄血）叭对血g=i2，挖其中址是拉格朗日插值

27、基函数。高斯型求积公式为ghc/）=Svw高斯型积分公式的优点是它的代数精度高，特别是对无穷区间或瑕积分更有效，但计算正交多项式的零点即积分节点有一定的工作量，好在数值学家们已算出一些特定的函数的积分节点和积分系数，在计算中我们可以查表直接得到这些数值。本章并不构造各种高斯型积分公式，有关的详细内容请参考有关的教材。下面给出小十上，取权函数盼)的高斯型积分。取一1,1上权函数貶)1的正交多项式为勒让德(Legendre)多项式:厶=亠加-1门用71乙也！ax2=0.5773503,=0.57735031.0000000,1.0000000=0.8611363,2=0.5773503,=0.57

28、735031.0000000,1.0000000=0.8611363,=0.33998100.3478548,0.65214524=0.3399810,=0.86113630.6521452,0.3478548x1=x5=0.90617980.2369269,0.23692695x2=x4=0.5384690.4786287,0.4786287=0.00.5688889高斯-勒让德相应的积分节点和积分系数表如下:上的积分E_a+ba则有要计算一般区间“切其中只需作变量代换，这样5仍可用高斯积分求积，即7=fw2cosadii例7.8应用两点高斯-勒让德积分公式计算山解：I=(0.5773503

29、)2cos(0.5773503)+(0.5773503)2cos(0.5773503)=0.558608例7.9例7.9应用两点高斯-勒让德积分公式计算解：，得到积分dt八弘4+（，得到积分dt八弘4+（“）=2(/()+/)=0.786885例7.10证明不存在%曲血使求积公式的代数精度超过T次。证明：只要能找到一个加次多项式，使求积公式两边不相等即可。用反证法，假定存在求积系数和节点及使求积公式对任何加次多项式畑精确成立。现取5代入求积公式左端得(x)f(x)dx=w(x)(x)dx0工=工=工監武e=而公式右端zj，故右端与左端不相等，与假设矛盾，说明不存在加次代数精度的求积公式，故高斯

30、型求积公式是具有最高代数精度的求积公式。7.6程序示例程序7.1复化梯形公式的自动控制误差算程序7.1复化梯形公式的自动控制误差算法。算法描述输入汉恥的值，定义畑；n=m,h=(bT2=TK=h”何+*0+逸)+扣3);T1=T2+100whileIT1T2IT1=T2；M-lA(xi+l/2)H=!或用复化梯形公式计算T2T2=(T1+H)/2；h=h/2；n=2n；endwhile输出T2。程序源码/Purpose：梯形公式的自动控制误差算法/includeincludedefinef(x)(sin(x)/f(x)由define定义definem100definea1.0defineb2.0defineepsilon0.00001voidma