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

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

2、方法是用函数的差商近似函数的导数，即取极限的近似值。下面是与式（7.1）相应的三种差商形式的数值微分公式以及相应的截断误差。向前差商用向前差商（平均变化率）近似导数有：                      （7.2）其中的位置在的前面，因此称为向前差商。同理可得向后差商、中心差商的定义。由泰勒展开 得向前差商的截断误差： 向后差商用向后差商近似导数有： 

3、;                    （7.3）与计算向前差商的方法类似，由泰勒展开得向后差商的截断误差：中心差商用中心差商（平均变化率）近似导数有：                    （7.4）由泰勒展开

4、60;得中心差商的截断误差：            差商的几何意义微积分中的极限定义，表示在处切线的斜率，即图7.1中直线的斜率；差商表示过和两点直线 的斜率，是一条过的割线。可见数值微分是用近似值内接弦的斜率代替准确值切线的斜率。图7.1 微商与差商示意图例7.1 给出下列数据，计算，0.020.040.060.80.105.065.075.0655.055.055解：（5.075.06）/（0.040.02）= 0.5（5.055.07）/（0.080.04）= 0.5（5.055

5、.055）/（0.080.10）= 0.25（(0.10) (0.06)）/（0.100.06）= 18.75设定最佳步长在计算数值导数时，它的误差由截断误差和舍入差两部分组成。用差商或插值公式近似导数产生截断误差，由原始值的数值近似产生舍入误差。在差商计算中，从截断误差的逼近值的角度看，越小，则误差也越小；但是太小的会带来较大的舍入误差。怎样选择最佳步长，使截断误差与舍入误差之和最小呢？一般对计算导数的近似公式进行分析可得到误差的表示式，以中心差商为例，截断误差不超过而舍入误差可用量估计（证明略），其中是函数的原始值的绝对误差限，总误差为当时，总误差达到最小值，即  &#

6、160;                   （*）可以看到用误差的表达式确定步长，难度较大，难以实际操作。通常用事后估计方法选取步长，例如，记为步长等于的差商计算公式，给定误差界，当时，就是合适的步长。例7.2 对函数，取不同的步长计算，观察误差变化规律，从而确定最佳步长。解：误差误差0.100.090.080.070.063.16303.16223.16133.16073.16000.00480.00400.00

7、310.00250.00180.050.040.030.020.013.15903.15883.15833.15753.15500.00080.00060.00010.00070.0032表中数据显示，当步长从0.10减少到0.03时，数值微分误差的绝对值从0.0048减少到0.0001，而随着 的进一步减少，误差的绝对值又有所反弹，表明当步长小于0.03时，舍入误差起了主要作用。在实际计算中是无法得到误差的准确数值的，这时以最小为标准确定步长，本例中取= 0.04。7.1.2 插值型数值微分对于给定的的函数表，建立插值函数，用插值函数的导数近似函数的导数。设为上的节点，给定，以为插值点构造插

8、值多项式，以的各阶导数近似的相应阶的导数，即当时，            （7.5）误差项为：例7.3 给定，并有，计算。解：作过的插值多项式：将代入 得三点端点公式和三点中点公式：    利用泰勒（Taylor）展开进行比较和分析，可得三点公式的截断误差是。类似地，可得到五点中点公式和五点端点公式：     7.1.3 样条插值数值微分把离散点按大小排列成，用关系式构造插值点的样条函数：当则当时，可用计算

9、导数。7.2数值积分在微积分中用牛顿莱布尼兹（Newton-Leibniz）公式计算连续函数的定积分： 但是，当被积函数是以点列的形式给出时，当被积函数的原函数难以得到时，例如，则无法用牛顿莱布尼兹积分公式计算。有时当被积函数的原函数过于复杂时，也不宜套用积分公式计算积分，而应采用数值积分公式计算定积分。在微积分中，定积分是黎曼（Rimann）和的极限，它是分割小区间长度趋于零时的极限，即在数值积分公式中，只能用有限项的和近似上面的极限，通常由函数在离散点函数值的线性组合形式给出。记，在本章中，用表示精确积分值，用表示近似积分值，称为积分节点，称为积分系数。确定中积分系数的过程就是构

10、造数值积分公式的过程。怎样判断数值积分公式的效果？代数精度是衡量数值积分公式优劣的重要标准之一。定义7.1 (代数精度)记上以为积分节点的数值积分公式为若满足而，则称具有阶代数精度。由此可知当具有阶代数精度时，对任意的阶多项式都有 。7.2.1 插值型数值积分对给定的被积函数在上的点列作拉格朗日插值多项式，以近似计算，即记，则有数值积分误差，也就是对插值误差的积分值 或对一般的函数，但若是一个不高于次的多项式，由于，而有。因此，阶插值多项式型式的数值积分公式至少有阶代数精度。例7.4 建立上节点为的数值积分公式。解：由得得以数值积分公式：7.2.2 牛顿-柯特斯（Newton-Cot

11、es）积分把积分区间分成等分，记步长为，取等分点作为数值积分节点，构造拉格朗日插值多项式，取由此得到的数值积分称为牛顿柯特斯积分。下面可以看到，牛顿柯特斯积分系数和积分节点以及积分区间无直接关系，系数固定而易于计算。梯形积分以和为插值节点构造线性函数，有那么，提取公因子后，得到牛顿柯特斯积分的组合系数：，它们已与积分区间没有任何关系了。记           （7.6）称为梯形积分公式。它的几何意义是用梯形面积近似代替积分值（图7.2）。图7.2 梯形积分怎样确定梯形积分公式的代数精度？我们可以

12、取验证取时，有即：取时，有即：取时，有得梯形求积公式具有一阶代数精度。梯形积分公式的误差：由  得  因为在上不变号，由积分中值定理得到梯形求积公式的截断误差：   （7.7）辛普森（Simpson）积分对区间作二等分，记。以和为插值节点构造二次插值函数，那么，有计算得到积分组合系数：，。   （7.8）S (f)称为辛普森或抛物线积分公式。它的几何意义是用过三点的抛物线面积近似代替积分的曲边面积（图7.3）。图7.3 抛物线积分面积分别将代入到和中，可以得到表明辛普森公式对于次数不超过三次的多项式准确成立，具有三阶代数精度

13、。因此可设一个三次多项式满足条件：，计算得到误差为：  于是有但   故辛普森求积公式的截断误差：（7.9）牛顿-柯特斯积分系数等分区间，取等分点为积分节点，其中。以为插值节点构造插值函数。 其中令，代入上式得（7.10）这里称为牛顿-柯特斯系数可见在取等距节点时，积分系数与积分节点和积分区间无直接关系，只与插值的节点总数有关，而在例7.3中的积分系数是待定系数，这就简化了数值积分公式，而不必对每一组插值节点xi都要计算一组相应的积分系数ai。在公式（7.10）中取=1，可算出梯形积分系数；取=2，可算出辛普森积分系数。在表7.1中列出从1到6的牛顿-柯

14、特斯系数。从表中可以看出牛顿-柯特斯系数具有对称性。表7.11          2        3      4    5  67.2.3 求积公式的收敛性与稳定性定义7.2 若，则称求积公式是收敛的。定义中的包含了，通常都要求计算积分的求积公式是收敛的。稳定性是研究计算公式当有误差时，的误差是否增长，现设，误差记为。定义7.3 对任给，只要，就有则称求积公式是稳定的。定理7.1 若求积公式的系数，则求积公式是稳定的。证明 由于,，故有于是对

15、，只要，就有故求积公式是稳定的。7.3复化数值积分由插值的龙格现象可知，高阶牛顿-柯特斯积分不能保证等距数值积分系列的收敛性，同时可证（略）高阶牛顿-柯特斯积分的计算是不稳定的。因此，实际计算中常用低阶复化梯形等积分公式。7.3.1 复化梯形积分把积分区间分割成若干小区间，在每个小区间上用梯形积分公式，再将这些小区间上的数值积分累加起来，称为复化梯形公式。复化梯形公式用若干个小梯形面积逼近积分比用一个大梯形公式效果显然更好，如图7.4所示。这种作法使我们想起定积分定义，即它为被积函数无限分割的代数和。这也正是计算定积分最朴素的算法。图7.4 复化梯形公式积分视图复化梯形积分计算公式对作等距分割

16、，有，于是在上，则有 记等分的复化梯形公式为，有 （7.11）复化梯形公式截断误差由，根据均值定理，当时，存在，有，于是   （7.12） 由此看到复化梯形公式的截断误差按照或者的速度下降，事实上，可以证明，只要在上有界并黎曼可积，当分点无限增多时，复化梯形公式收敛到积分。记 ，则有对于任给的误差控制小量，有   或   就有，式中表示取其最大整数。7.3.2 复化辛普森积分把积分区间分成偶数等分，记，其中是节点总数，是积分子区间的总数。记，在每个区间上用辛普森数值积分公式计算，则得到复化辛普森公式，记为

17、。复化辛普森积分计算公式而，称    （7.13）为复化辛普森积分公式，它是在上采用辛普森积分公式叠加而得。下面用图7.5显示复化辛普森积分计算公式中节点与系数的关系，取，在每个积分区间上提出因子后，三个节点的系数分别是1，4，1；将4个积分区间的系数按节点的位置累加，可以清楚地看到，首尾节点的系数是1，奇数点的系数是4，偶数点的系数是2。141141141141142424241图7.5  复化辛普森积分系数复化辛普森公式的截断误差设，在上的误差为因此，即        &

18、#160; （7.14）与复化梯形公式类似，误差的截断误差按照或者的速度下降。可以证明，只要在上有界并黎曼可积，当分点无限增多时，复化辛普森公式收敛到积分。记，则有对任给的误差控制小量，只要    或    就有。例7.5 求，计算中要求有5位有效数字。用复化梯形和复化辛普森求积公式的分点应取多少？解：      由复化梯形误差公式得到：计算出，复化梯形公式至少要在.00等分n = 68。由复化辛普森误差公式，有在复化辛普森公式中取 或。   &

19、#160;   7.3.3 复化积分的自动控制误差算法复化积分的误差公式表明，截断误差随分点的增大而减小，对于给定的误差量，用估计函数导数的界的方法可计算出。用误差公式计算满足精度的分点数，像是在做一道计算导数上界的微积分习题（如例7.5所示）。但是在实际运算中，一般难以估计出函数的各阶导数界，也就无法确定分点数。在计算中常用误差的事后估计方法，即用估计误差。T2n (f )的计算公式对定积分，取分点，计算得取分点，计算得这里，。可以看到，的值是与新增分点的组合。取分点，计算得这里，。同理，计算时只要在的基础上计算新增分点，的值再做组合，如图7.6所示。图7.6 与一般地，

20、每次总对前一次的小区间分半，分点加密一倍，并可充分利用老分点上的函数值，每次只需计算新增分点的和。对上等分， ，则有记上的中点为则  （7.15）其中。或        其中。类似地，可得积分节点为，的辛普森求积公式的关系式：        （7.16）其中：由误差公式：由于，分别为及个点上的均值，可视，于是 上式表明的误差大约是误差的4倍。或       

21、   （7.17）由此得到启发，对任给的误差控制量，要，只需即可，而用作为控制手段简单直接，序列在计算机上也不难实现。复化积分的算法描述从数值积分的误差公式可以看到，截断误差随分点的增长而减少，控制计算的精度也就是确定分点数。在计算中不用数值积分的误差公式确定分点数的理论模式，而用作为控制，通过增加分点自动满足精度的方法称为数值积分公式的自动积分法。即在计算中构造序列，直到或时停止计算，由分点数自动控制积分值的误差，并取。下面描述复化数值积分公式的自动控制误差算法，详细程序和算例请看本章7.6节。1输入：误差控制精度e = eps；初始分点值 。2计算分点的复化梯

22、形积分T1=T2+100      /迭代计算中T1和T2分别表示和3while | T1T2|eT1 = T2H=Hn           /计算新增节点的值T2= (T1+H )/2H = h/2，n =2n    /将区间一分为二end while4输出积分值T2。在自动控制误差算法中初始分点值不宜过小，以防假收敛。7.3.4 龙贝格（Romberg）积分由前面得到的关系式（7.17），可以将 作为的修正值补充

23、到，即 （7.18）其结果是将梯形求积公式组合成辛普森求积公式，截断误差由提高到。这种手段称为外推算法。外推算法在不增加计算量的前题下提高了误差的精度，是计算方法中一种常用手法。不妨对再做一次组合。由得到          （7.19）复化辛普森公式组成复化柯特斯公式，其截断误差是。同理对柯特斯公式进行组合： 得到具有7次代数精度和截断误差是 的龙贝格公式：还可以继续对做上去。为了便于在计算机上实现龙贝格算法，将统一用表示，列标 分别表示梯形、辛普森、柯特斯积分，行标表示分点数或步长j。龙贝格计

24、算公式：对每一个从2做到，一直做到小于给定控制精度停止计算。龙贝格算法龙贝格算法按表7.2元素的行序进行运算， 在计算中每个元素只用到上一行和本行的元素。对上面的算法进一步优化，对每k行可将计算定义在两行元素之间，令；在每计算一行元素后，要将。表7.2 龙贝格算法计算元素顺序表            1输入区间端点 ，精度控制值，循环次数，定义函数，取；2；3for    =2  to          4输出。7.4重积分计

25、算在微积分中计算二重积分是用化为累次积分的方法进行的。计算二重数值积分也是计算累次数值积分的过程。为了简化问题，我们仅讨论矩形域上的二重积分。有很多非矩形域上的二重积分可作变换将其转换到矩形域上。                （7.20）其中：是常数，在上连续。像在微积分中一样，将二重积分化为累次积分：      （7.21）或     

26、0;      二重积分的复化梯形公式对区间和分别选取正整数，在轴和轴上分别有步用复经梯形公式计算，计算中将当作常数，有 （7.22）再将当作常数，在方向上计算式（7.23）中每一项的积分，有=则+=积分区域的4个角点的系数是1/4，4个边界的系数是1/2，内部节点的系数是1。误差：和在积分区间内。例子7.6 用复化梯形公式计算二重积分，取。解：如表7.3所示：表7.3 数值表的数值如表7.4所示：表7.4 数值表=0.873601的准确值是0.886176。二重复化辛普森求积公式对区间和分别等分和等分，在轴和轴上分别有步长 

27、;均为偶数类似于二重复化梯形公式推导，得记  误差： 和在积分区间内。按例7.6的化分区间，的值如表7.5所示：表7.57.5  高斯（Gauss）型积分公式介绍对插值型积分公式在牛顿-柯特斯积分公式中要求节点是等距的，其优点是计算积分系数的公式规则相同，缺点是制约了求积公式的代数精度。可以证明：当节点个数为偶数时，求积公式具有阶的代数精度；当节点个数为奇数时，求积公式具有阶的代数精度。如果我们不预先指定求积节点的位置，和权系数都作为待定的常数，能否适当地确定它们，以提高积分公式的代数精度？回答是肯定的。个待定常参数，需要个方程来确定，取一个函数组：，这一组函数构

28、成了次多项式的基，任一小于等于次的多项式，都可以用这组函数的线性组合来表示。如果某一积分公式，对这组函数都能精确积分，则此积分公式就有次代数精度。例7.7 计算求积系数和求积节点，使得至少具有3阶代数精度。解：按照求积公式的代数精度定义，分别令，得方程组：解方程组得：  ， 求积公式：按例7.7的方式，构造更高阶的代数精度的求积公式，生成求积系数和求积节点的方程组并无困难，而求解该方程组则无一定的章法可循。一般地，通过计算正交多项式的零点作为求积节点。当取积分节点为正交多项式的零点时，则节点个数是的求积公式具有阶的代数精度。并称积分节点为正交多项式的零点的数值积分公式为高斯型积分公式

29、。为了一般性，考虑积分其中称为权函数。当时，即是普通的积分。对于不同的权函数选定的节点也不相同。如何构造高斯型积分公式呢？对给定的及权函数，由施密特（Schmidt）正交化过程作出正交多项式；解出正交多项式的个零点，这个个零点就是积分节点；以这些节点构造插值多项式，计算积分系数其中是拉格朗日插值基函数。高斯型求积公式为高斯型积分公式的优点是它的代数精度高，特别是对无穷区间或瑕积分更有效，但计算正交多项式的零点即积分节点有一定的工作量，好在数值学家们已算出一些特定的函数的积分节点和积分系数，在计算中我们可以查表直接得到这些数值。本章并不构造各种高斯型积分公式，有关的详细内容请参考有关的教材。下面

30、给出上，取权函数的高斯型积分。取1，1上权函数的正交多项式为勒让德（Legendre）多项式：高斯-勒让德相应的积分节点和积分系数表如下：2 =0.5773503， =0.57735031.0000000，1.00000004=0.8611363，=0.3399810 =0.3399810， =0.86113630.3478548，0.65214520.6521452，0.34785485x1= x5=0.9061798x2= x4=0.538469=0.00.2369269，0.23692690.4786287，0.47862870.5688889要计算一般区间上的积分，

31、只需作变量代换则有，其中，这样 仍可用高斯积分求积，即 例7.8 应用两点高斯-勒让德积分公式计算。解：I = (0.5773503)2cos (0.5773503) + (0.5773503)2 cos (0.5773503)=0.558608例7.9 应用两点高斯-勒让德积分公式计算。解：令，得到积分例7.10 证明不存在使求积公式的代数精度超过次。证明：只要能找到一个次多项式，使求积公式两边不相等即可。用反证法，假定存在求积系数和节点及使求积公式对任何次多项式精确成立。现取代入求积公式左端得而公式右端，故右端与左端不相等，与假设矛盾，说明不存在次代数精度的求积公式，故高斯型求积

32、公式是具有最高代数精度的求积公式。7.6程序示例程序7.1  复化梯形公式的自动控制误差算法。算法描述输入的值，定义；   ；T1=T2+100while |T1T2|T1 = T2；H = h        ！或用复化梯形公式计算T2T2 = (T1 + H) /2；h= h/2；n = 2n；end while输出T2。程序源码/  Purpose：梯形公式的自动控制误差算法      /# include &l

33、t;stdio.h># include <math.h># define f (x)  (sin (x)         / f(x)由define定义# define m  100# define a 1.0# define b 2.0# define epsilon 0.00001void main( )  int n = m；   int i；   double T n，H_n，T1，T2；   dou

34、ble h = (ba) / n；   T_n = (f (a) + f(b) /2；   for (i =1；in；i+       T_n + = f(a+h i)；       T_n= h；   T2 = T_n；   T1 = T2+100；   do     T1 = T2；    for (i=0，H_n

35、=0；i =1；in；i+          H_n+=f (a+hi+h/2)；      H_n= h；      T2 = (T1 + H_n) /2h = h /2；n = n2；        while (fabs (T1T2)epsilon)       

36、 printf (“T=% 1f n”，T2)；          / End of File计算实例对于，在区间上验证梯形公式的自动控制误差公式。程序输入输出对于不同，修改程序# define f(x)项，本例，初始，区间。）        T = 0.956447程序7.2  龙贝格积分算法计算公式和算法描述第7.3.4节所述。程序源码/       Purpose：龙贝格

37、算法        /# include<stdio.h># include<math.h># define f (x)             (sin (x)# define N_H             20# define MAXREPT        10# define a                 1.0# define b