第7讲(2014-11-12).doc

第4章 数值积分与数值微分上次课要点：1 数值积分的基本概念与基本方法一、求积公式的一般形式其中仅依赖于节点而与无关。更一般的形式：由节点处的函数值、一阶导数值（甚至高阶导数值）的线性组合。二、有关求积公式的几个概念1代数精度定义1 如果某个求积公式对于次数不超过的多项式都能精确成立，但对于次多项式不精确成立，则称该求积公式具有次代数精度。2求积公式的收敛性定义2 在求积公式中，若其中，则称求积公式是收敛的。3求积公式的稳定性一个求积公式是数值稳定的，是指当输入数据的误差很小时，输出数据的误差也很小。三、插值型求积公式插值型求积公式:式中。公式余项:定理1 具有个插值节点的插值型求积公式其中，为Lagrange插值基函数，至少具有次代数精度。 定理2 若插值型求积公式中系数，则此求积公式是数值稳定的。 2 牛顿科特斯公式：等分，步长：，节点：牛顿科特斯公式:其中称为科特斯系数。可以证明，当为偶数时，牛顿科特斯公式具有次代数精度。（1）梯形公式（2）辛普森公式（3）科特斯公式这里。当时，科特斯系数有正有负，数值稳定性没有保证。三、几种低阶求积公式的余项1、梯形公式的余项2、辛普森公式的余项3、科特斯公式的余项。3 复化求积公式一、复化梯形公式余项：。二、复化辛普森公式 余项：三、埃尔米特积分法余项：。本次课继续。4 龙贝格求积公式(Romberg Integration)实践中，固定复化求积公式中的步长会导致：（1）步长取得太大，精度难以保证；（2）步长太小，则计算量太大，有时不必要。一种可行的解决思路：采用变步长方法。即在步长逐步分半的过程中，反复使用复化求积公式，直到所得的积分值满足精度要求为止。一、梯形法的递推化将区间进行等分，设节点为，则梯形法求积公式为再将每一个子区间二等分，记相应的中点为，则此时梯形法求积公式为于是有递推公式。注意：上述递推公式中不变。二、龙贝格算法的思路上述算法过于依赖梯形公式，效率较低。问题：能否利用和的某种组合，得到精度更高的求积公式？理查德外推法：用精度较低的近似公式，组合成精度较高的近似公式，在数值计算的许多问题中都有应用。理查德在1927年发现的外推法，当时并没有用来做数值积分问题。而数学家龙贝格在1955年将它应用到数值积分上，取得很好的成效。下面考察梯形法、辛普森法和科特斯法，当步长减半时误差的变化情况。可以看出，当减半时，梯形法：误差大致减至原有误差的；辛普森法：误差大致减至原有误差的；科特斯法：误差大致减至原有误差的。由上述讨论，可得移项整理可得说明当很小时，也很小。因此，实践中常给定计算精度，用作为停止计算的准则。进一步，由于可以期望，用计算积分值可能得到更好的结果，而不需要增加计算量。定理4 梯形法二分后的两个积分值的线性组合等于辛普森法积分值，即.用同样方法，依据辛普森公式可导出科特斯公式依据科特斯公式可导出龙贝格公式在变步长过程中，应用上述公式，就能将粗糙的梯形公式逐步加工成精度较高的辛普森公式、科特斯公式和龙贝格公式。三、理查德外推法与龙贝格算法将区间等分后，设子区间的长度为，用梯形法计算的积分值记为。定理5 设，则有（1）其中，系数与无关。利用公式（1），可得（2）记 ，则 （消去了项） （3）其中也与无关。由（3）得（4）记，则又可消去项，得如此进行下去，每加速一次，误差的量级便提高2阶。一般地，记，则有递推公式上述方法称为理查德外推加速法。以表示二分次后的梯形值，表示的次加速值，则有龙贝格算法：步1 初始化：计算，置（记录区间的二分次数），；步2（二分） 计算积分值：；步3（加速） 求加速值：计算，；步4 精度检验：对指定的精度，若，则终止计算，并取作为所求的结果；否则置，转步2。计算次序：第1次循环二分：加速：第2次循环二分：加速：第3次循环二分：加速：第4次循环二分：加速：5 自适应积分方法一般而言，用复合求积方法求积时：1）被积函数变化剧烈将区间分得细一些；2）被积函数变化平缓不需要分得很细。问题：如果被积函数在积分区间上，有的部分变化剧烈，有的部分变化平缓，则将整个积分区间一起细分是非理性的。解决方法：根据各个子区间满足精度的情况，采用不同的步长这就是自适应积分方法。下面以复合辛普森公式为例，介绍自适应积分方法。设给定精度要求，计算积分的近似值。先取步长，记应用辛普森公式，有（1）将区间二等分，步长，在每一个小区间上用辛普森公式，得，记则有即（2）如果在区间上变化不大，可假设，此时，由（1）、（2）两式可得由（2）式，得比较上述两式，得这里。如果有 （3）则可期望得到此时可取作为的近似。 考虑到前面假设不一定合理，（3）式改为 （4）成立时，认为也成立。如果（4）式不成立，则分别对区间和重复上述过程，此时步长。由于区间长度减半，对每一个子区间而言，（4）式右端的应换成。对于已经满足精度要求的子区间，停止细分，记录该子区间上的近似积分值；对于不满足精度要求的子区间，继续二分。8 数值微分研究内容：用数值方法求函数在离散点处的导数值。解决方法：用某些点处函数值的线性组合，近似指定点处的导数值。问题分为两类：（1）求一点处的导数值；（2）求等距分布的多个点处的导数值。一、泰勒展开法函数展开法是基于函数的泰勒展开式求导数的近似值的方法。1、显式求导公式（用于求一个点处的一阶、二阶导数）（1）（2）由（1），保留等式右边的前两项，得忽略掉，得一阶向前差商公式 （3）公式（3）用于求区间左端点的近似导数。同理，由（2）可得一阶向后差商公式 （4）公式（4）用于求区间右端点的近似导数。由，得忽略掉，可得一阶中心差商公式 （5）公式（5）用于求区间内点的导数值。类似地，得从而可得二阶中心差商公式（6）公式（6）用于求区间内点的二阶导数值。从截断误差看，越小，公式误差越小。但上述公式的数值稳定性不好。以一阶中心差商公式为例，设为理论值，为计算值，则输入数据的舍入误差为，考察一阶中心差商公式输出数据的误差分析：与是两个相近数，它们相减，会损失有效数字的位数，导致舍入误差较大；在此基础上再除以，又将上述舍入误差放大了倍，雪上加霜；上述公式在与异号时尤甚（虽然与随变化，但具有保号性）。综上所述，不能太小，也不能太大。2、隐式求导公式已知：节点：函数值：边界条件：和求：。由泰勒展开式（1）（2），整理后，得而 ，可得。用表示的近似值，得近似等式即这是由个未知数、个方程组成的线性方程组。考虑到两个边界条件是已知的，再令，便得（大约）具有精度的一阶导数隐格式方程组上述方程组主对角占优，解存在且唯一。用追赶法求解是数值稳定的。再看二阶导数。由于而，于是用表示的近似值，得近似等式即再加上两个边值条件，并令便得（大约）具有精度的二阶导数隐格式方程组上述方程组同样是主对角占优的，解存在且唯一，用追赶法求解数值稳定。二、插值型求导公式过已知点，建立插值多项式，用的导数作为的导数的近似值，即。（一）用Lagrange插值多项式求近似导数常用的有三点求导公式和五点求导公式。1、三点公式已知三个点，则的导函数因此，各点处的一阶导数近似公式：2、五点公式已知五个点用三点类似的方法可得（二）用三次样条插值求导（求多个点处的导数值）（略）三、利用数值积分求导设，则上式右边的积分可以采用不同的求积公式，从而得到不同的求导公式。1、中矩形法（用于求一点处的导数）对用中矩形公式则得即。2、辛普森法（用于求多点处的导数）对用辛普森公式则得用表示的近似值，得这是由个未知数、个方程组成的线性方程组。若已知，则可得到前面相同的一阶导数隐格式。四、数值微分的外推算法外推法是一种逐步提高公式精度的方法，使截断误差逐步减小。由泰勒展开式 （1） （2）得 （3）其中只与有关，与无关。令，则（4） （5）由得（6）可见，用精度更高，余项为。令，可以继续外推。一般地，令，可得，误差：。当较大时，截断误差很小。考虑到舍入误差，一般不能太大。五、复数法求近似导数William Squire , George Trapp. Using Complex Variables to Estimate Derivatives of Real Functions. SIAM Review, Vol.40, No.1, PP110-112, March 1998.算法要点：将按复数展开，表示虚数单位。则虚部令 。该式计算的截断误差为。没有两个相近数相减产生的有效数字位数的损失。上述算法与中心差商法进行数值比较 （1） （2）例1 ，公式（1）公式（2）0.1D-010.1D-020.1D-030.1D-040.1D-050.1D-060.1D-070.1D-080.1D-090.1D-0100.1D-0110.1D-0120.1D-0130.1D-0140.1D-0150.18602018344501897D+020.18600824790342307D+020.18600812854818738D+020.18600812735480865 D+020.18600812734248517 D+020.18600812735081185 D+020.18600812723423843 D+020.18600812723423843 D+020.18600814222224926 D+020.18600815332447951 D+020.18600898599174798 D+020.18601231666082185 D+020.18585133432225120 D+020.18596235662471372 D+020.20539125955565396 D+020.18599607128036329D+020.18600800678177631 D+020.18600812613698936 D+020.18600812733054151 D+020.18600812734247702 D+020.18600812734259637 D+020.18600812734259757 D+020.18600812734259759 D+020.18600812734259759 D+020.18600812734259759 D+020.18600812734259759 D+020.18600812734259759 D+020.18600812734259759 D+020.18600812734259759 D+020.18600812734259759 D+02可以看出，公式（1）在时比较准确。例2 ，公式（1）公式（2）0.100000E-10.100000E-20.100000E-3