计算方法 第5章 数值积分与数值微分

第五章 数值积分与数值微分在高等数学中我们学过定积分的计算方法，若找到被积函数在区间上的一个原函数，利用Newton-Leibniz公式可以轻易得计算出积分值，但在实际问题中，往往会遇到一些困难。有些函数虽然能找到原函数,但表达式过于复杂，例如的原函数为有些函数找不到初等函数形式的原函数，例如积分有些情况下，函数值是用表格形式给出的，例如：对于以上这些积分问题，解决的方法就是使用数值积分方法。其实数值积分方法不仅可以解决上述问题，最为重要的优点是对任意被积函数任意积分区间的积分问题都可以采用统一的数值积分公式，非常便于计算机编程实现。 对于微分问题，虽然对每一个初等函数都可以求出其导数，但是不同函数其求导方法依赖于各自不同的求导公式，没有简单、统一的处理方法，而数值微分法却可以对不同的函数使用统一的数值微分公式或数值微分算法。 本章首先介绍一些数值积分公式，最后再简单的介绍数值微分问题。5.1 数值积分公式 1.数值积分的基本思想 我们知道定积分的几何意义就是所围成的曲边形面积，而数值积分的基本思想是利用函数在区间上某些点处函数值的线性组合来计算其定积分的近似值，把计算定积分这一复杂问题转换为仅仅涉及到函数值的计算问题，而无需考虑函数本身的结构以及函数值的真实来源，这样就很便于计算机编程实现。 为求的数值解，我们考虑在积分区间上用的插值函数近似，则有 (5.1)称(5.1)式为函数在区间上的求积公式，为求积系数，插值函数的插值节点在这里称为求积节点，如果记 (5.2)称(5.2)式为求积公式(5.1)的截断误差，实质上求积公式的截断误差就是插值余项所致，求积公式的精确程度可由截断误差的绝对值来度量。从(5.1)式还容易看出求积系数仅仅与节点的选取有关，而不依赖被积函数，因此求积公式适用于所有函数，当然其截断误差的大小与节点的选择和函数本身的特性均有关系。 2.求积公式的代数精确度 定义5.1如果求积公式，对于是精确成立的，而对于不再精确成立，则称该求积公式的代数精确度为。根据定义5.1，如果求积公式的代数精确度为，则求积系数应满足 ……………. 如果将求积公式中的系数作为未知量，则我们已经列出了确定求积系数的线性方程组，选出节点，取，求解此线性方程组，即可以通过求解此方程组确定出求积系数，得到其求积公式且求积公式的代数精确度至少为。 例5.1试确定中的待定参数，使该积分公式的代数精确度尽可能的高，并指出所求积分公式的代数精确度。 解：1)求待定系数，利用求积公式，得到关于的方程组 解之得 2)确定求积公式的代数精确度，按此方程组确定的待定系数，求积公式的代数精确度至少为2，现在依次验证时的求积公式是否精确成立，直到出现首个不成立为止。 则直到均精确成立，但对于不再精确成立，表明所求的求积公式的代数精确度为5。用解线性方程组的方法构造求积公式很不方便，通常是利用以为插值节点的插值多项式来构造出常用的求积公式，在编程中直接使用这些常用求积公式，以下我们将构造这些求积公式。3．梯形求积公式对于，取，用近似，则有称积分公式 (5.3)为梯形求积公式。梯形求积公式具有直观的几何意义：图5.1梯形求积公式的几何意义关于梯形求积公式的截断误差，有如下结论。定理5.1若函数在上有连续的二阶导数，则梯形求积公式的截断误差为 (5.4) 证明：由于梯形求积公式的误差是由插值余项而引起的，其余项公式为其中与有关，再注意到在上有连续的二阶导数，记，则在区间内且不变符号，由广义积分中值定理有 ，对上式右端求积分可得例5.2验证梯形求积公式的代数精确度为1。 解： 对于精确成立，但对于不成立，所以代数精度为1，从梯形公式的代数精确度太小也可以看出梯形求积公式的精度不高。 4.Simpson求积公式 对于，将区间划分为二等分,令，取，用近似，则有 称积分公式 (5.5)为Simpson求积公式，Simpson求积公式的几何意义：图5.2Simpson求积公式的几何意义关于Simpson求积公式的截断误差，仿照前边定理5.1的证明过程可证明如下定理。定理5.2若函数在区间上有连续的4阶导数，则Simpson求积公式的截断误差为 (5.6)例5.3求Simpson求积公式的代数精确度。解：对于，考察 所以Simpson求积公式的代数精确度为3，从代数精确度、截断误差和实际计算效果(稍后的例5.4比较了3个求积公式的计算结果)看，Simpson求积公式都优于梯形求积公式。 5.Cotes求积公式 对于，将区间划分为4等分，令，取，用近似，则有称积分公式 (5.7)为Cotes求积公式，其截断误差为 (5.8)可以验证Cotes求积公式的代数精度为5。 6．Newton-Cotes求积公式 现在我们讨论一般情况下的求积公式，先将积分区间划分为等分，令，积分节点为，用近似，则有 (5.9)称此公式为Newton-Cotes求积公式，其中 (5.10)令，则对于有称为Cotes系数。可以证明当n为偶数时，Newton-Cotes型求积公式的代数精确度为n+1，当n为奇数时，其代数精确度为n。Cotes系数表从Cotes系数表可以看出，梯形求积公式、Simpson求积公式均是Newton-Cotes积分公式的特例。另一方面从Cotes系数的表达式还可以看出，Cotes系数与积分区间无关，与被积函数也无关，仅与区间的等份个数有关。Cotes系数具有如下性质。由于Cotes系数是由插值公式推导而来，而插值公式在高阶情况下具有不稳定性，因此Cotes公式必然也是不稳定的，从Cotes系数表中的出现负值就可以看出问题，通常情况下，Newton-Cotes求积公式中最常用的是梯形求积公式、Simpson求积公式和4阶Cotes求积公式，通常不用高阶的Newton-Cotes求积公式。 例5.4分别用梯形求积公式、Simpson求积公式和4阶Cotes求积公式计算定积分的近似值。 解：梯形求积公式： Simpson求积公式： Cotes求积公式： 和精确值比较，梯形求积公式、Simpson求积公式和4阶Cotes求积公式所求结果分别有2、4、6位有效数字。5.2 复化求积公式 Newton-Cotes求积公式的两个问题高阶求积公式稳定性难以保证；从截断误差公式可看出，低阶求积公式对于积分区间较大时同样精度也难以保证。解决方法： 针对前述的两个问题，注意到前节中各截断误差公式中均有因子，如果将区间缩小定会提高积分的精度，复化求积公式就是将大的积分区间分解为小的积分区间，在小的积分区间上使用低阶求积公式。复化梯形公式将积分区间划分为等分，则每个子区间长度均为，称为步长，此时分点为，在每个子区间上施行梯形公式，得到称 (5.11)为复化梯形公式。由定积分的定义，显然有此极限也说明，当N充分大时，复化梯形公式会充分接近积分的精确值。复化梯形公式的截断误差由以下定理给出。定理5.3若函数在上有连续的二阶导数，则复化梯形公式截断误差为 (5.12)其截断误差的上限为 (5.13) 证明：复化梯形公式的截断误差就是每个小区间上的截断误差之和，根据(5.4)式有，由于在上有连续的二阶导数，根据高等数学中的介值定理，存在使从而有，截断误差限估计式(5.13)式是显然成立的。复化Simpson公式类似于复化梯形公式的构造方法，再构造复化Simpson公式，将积分区间划分为等分，则每个子区间长度均为，在每个子区间上施行Simpson公式(此时有个分点)，得到称 (5.14)为复化Simpson公式，对于复化Simpson公式同样也有这也表明当N充分大时，复化Simpson公式也能很好的近似定积分的精确值。关于复化Simpson公式截断误差，有比复化梯形公式更好的结论。 定理5.4若函数在上有连续的4阶导数，则复化Simpson公式截断误差为 (5.15)其截断误差上限 (5.16)复化Cotes公式类似于复化梯形公式和复化Simpson公式的构造过程可以得到复化Cotes公式 (5.17)其中，称该式为复化Cotes公式，容易看出复化Cotes公式共有4N+1个分点。同样根据定积分的定义也有 定理5.5若函数在上有连续的6阶导数，则复化Cotes公式截断误差为 (5.18)且截断误差的上限为 (5.19) 例5.5试取9个节点分别用复化梯形公式、复化Simpson公式、复化Cotes公式计算定积分。 解：首先建立如下表格根据此表可以计算出的结果如下： 所求积分的精确值为，则的误差分别为，，，则有效数字分别为2、7、6位。此例也说明同样计算9次函数值，计算量相同，使用不同的公式所得结果的精度相差很多。 在复化积分公式中，编者推荐使用复化Simpson公式，此公式简单易记且精度高。 例5.6如果计算积分的精度要求为，试问分别用复化梯形公式、复化Simpson公式、复化Cotes公式，积分区间至少应划分多少等分？ 解：，则 复化梯形公式： 复化Simpson公式： 复化Cotes公式： 因此可以得到如下结论，对于所给的定积分，到达指定的精度，复化梯形公式需要674等分积分区间，也就是将要计算675个分点的函数值；复化Simpson公式需要7等分积分区间，总共需要计算个分点的函数值；复化Cotes公式需要2等分积分区间，总共需要计算个分点的函数值。5.3 Romberg求积法复化公式的一个缺点是必须先给出步长，从截断误差公式可以看出，复化求积公式的计算精度严重依赖于步长，而如何选取步长也没有给出一个可行的方法。本节所给出的求积公式采用变步长的方法，把步长逐次分半，直至达到要求的精度为止，这就是本节所要讨论的内容—Romberg求积法。变步长梯形公式变步长梯形公式的思想是，首先取步长是整个积分区间，使用梯形公式计算积分值，逐次将步长减半，使用复化梯形公式计算积分值，当积分值的计算精度到达要求时停止计算，否则为步长继续减半。在算法实施过程中，我们要解决三个问题。1)在计算出后，如何计算，即导出与的递推公式。2)如何避免重复计算相同节点函数值？3)在计算出后，如何估计其误差，即算法的终止准则是什么？首先推导变步长梯形公式的递推公式，在计算时，需要计算个点上的函数值。在求得后再计算时，需将每个子区间再二等分，新增加个分点。为了避免重复计算，计算时，已经计算出的个点上的函数值应保留下来，而只计算新增加的个分点上的函数值，为此把表示成两部分之和。即这样就可以得到积分值递推公式 (5.20)有了递推公式同时也就解决了重复计算问题。以下我们根据和的截断误差，推导出算法的终止准则，按照(5.12)式有这里，如果在积分区间内的变化不大，可以假设，综合以上两式则有该式表明，用作为定积分的近似值，其误差大约为，因此终止准则可以用和的差来度量定积分计算值的精度。即对于给定的精度，算法终止准则为 (5.21)例5.7用变步长梯形公式求定积分，要求终止准则为。 解：首先用梯形公式得到，再利用递推公式得到、、、…直到满足精度要求为止，计算结果如下，，，，精确值为，。Romberg(龙贝格)求积公式我们先利用前边推导终止准则时得到的近似公式 (5.22)对例5.7中的和数据进行计算该值与精确值的误差为，达到的四位有效数字，实际上好于的精度。 基于这一事实，我们猜想把(5.22)这一近似公式作为新的积分近似值应该优于，正是这一猜想产生了非常有效的Romberg积分法。首先我们可以验证 (5.23)仿照对梯形公式的推导过程，对做同样的处理，则有如果在积分区间内的变化不大，可以假设，则有通过验算还可以得出 (5.24)继续对做同样处理可推出还可以得到 (5.25)可以想象要优于，用计算积分其精确将可能更高，如此递推所得的递推公式称为Romberg递推公式。Romberg递推公式的算法设计通常将该计算步骤设计成以下表格来递推计算积分值便于编程的Romberg递推公式为 ，注意在此递推格式中内循环的是同步变化的，直到完成为止。在编程时可以将递推公式中的作以下等价变换其中后一项就是本次递推的修正项，也可以作为终止精度的判断项。如同第4章计算差商表一样，最新一行近似积分值的存放方式按倒序更方便，即每次迭代对角线元素为首个元素，算法结束时首元素就是所求积分值，设计算法如下：算法5.1Romberg求积法[初始化]算法中将最后一行按反序存入向量，令初始步长，给定终止精度，，计算首个梯形值；令，；[更新梯形值]；[递推]令，然后对循环3.1)3.2)[判断终止准则]如果，算法结束，为满足精度要求的积分近似值。，转1)。[算法结束]例5.8用Romberg递推公式求定积分，终止精度为。 解：同样先用梯形公式得到，再利用Romberg递推公式得到以下结果20.78333330.03333330.002064840.78552940.00219610.000131280.78539640.00013300.0000017建议：对于已知函数表达式的情况计算定积分推荐使用Romberg算法，对于给出的仅仅是离散函数值的情况推荐使用复化Simpson公式。Richardson(理查森)外推加速法前面根据猜想的方式得到了Romberg求积公式并验证了有效性，现在我们探讨其理论依据，Romberg求积法本质上是用误差补偿的加速算法，由收敛较慢的序列构造成收敛很快的序列，此方法即是Richardson外推加速法。记是以为步长的梯形值，可以证明当被积函数有任意阶导数时，有 (5.26)其中为积分的理论值，与无关，该式说明，在被积函数和积分区间确定的情况下，Romberg递推公式中的每一步积分值可以看作是步长的函数。 从(5.26)式可以看出，，即用复化梯形公式计算积分，其计算误差阶为。如果在式中用来代替，则得到 (5.27)用4乘(5.27)式减去(5.26)式再除则消除了项，记为，则得 (5.28)显然对理论积分值的精度达到。前边我们已经验证过是以为步长的复化Simpson公式，如果在中再用来代替又可以得到此时我们再消除其中的项，用乘，减去(5.28)式再除将其记为，则得到下式显然对理论积分值的精度达到，就是步长为的复化Cotes公式。依次类推就得到Romberg求积法。我们从递推过程可以看出，每次递推其步长将折半，对积分理论值的精度又提高两阶，这就是Romberg求积法高效的原因所在。其实Richardson外推加速法不仅可以对梯形积分值序列加速，只要是满足以上特征的所有序列均可加速，在下节数值微分中还会用到Richardson外推加速策略，现在我们看另外一个例子，体会一下Richardson外推加速的用途。例5.9已知有如下恒等式试依据此式的值按，用Richardson外推法求的近似值。 解：设，则有使用Richardson外推法外推一次得到再外推一次可得到将计算结果列于如下表中精确值为，误差为。5.4 数值微分 所谓数值微分就是仅利用函数在一些点处的函数值求出其导数的近似值，其优点是可以对不同的函数使用同一算法求各自对应的导数，便于计算机编程实现。本节将讨论数值微分的差商法、Richardson外推法、插值法等。差商代替导数为计算的导数，取一较小的步长，考虑的一阶Taylor展开从该式可以解出 (5.29)同样如果考虑的一阶Taylor展开则有 (5.30)如果再考虑和的Taylor展开两式相减则有 (5.31)以上(5.29)、(5.30)、(5.31)三个公式分别称为向前差商、向后差商、中心差商公式，从无穷小项的阶数可以看出中心差商近似导数更为精确。如果计算是精确的，显然愈小差商愈能逼近导数值，但从误差分析的角度讨论，太小将意味着两个相近的数相减，势必导致有效数字的损失，这一点我们在第一章的例1.2已经看到过，现在再看一个例子。例5.10用中心差商计算在处的导数。解：使用中心差商公式(5.31)，有取以下不同的步长，按5位有效数字进行计算，结果见如下表格(精确值0.3535533906)按5位有效数字进行计算在时误差最小，达3位有效数字，误差也不算大，反而是时仅1位有效数字。读者可以自己编程试算一下，如果用双精度进行计算，取，计算将达到8位有效数字，当后误差将逐步增大；如果用单精度进行计算，时精度最高计算达到6位有效数字，当后将无一位有效数字。从误差分析的角度，本例的计算应该按如下方式用计算已达4为有效数字，用5位有效数字计算，达到4位有效数字已算相当精确。Richardson外推法计算微分使用差商近似导数并不是一个让人满意的选择，但是对中心差商进行外推却可以得到很好的算法，我们考察中心差商公式：根据Taylor展开公式，有两式相减并除以h得这里常数与无关，显然符合Richardson外推形式，从而有如下外推公式 (5.32)外推格式表如下所示注意到对的计算，当太大也就是步长太小时将导致两个非常相近的数相减，造成有效数字的损失，因此使用外推要适可而止，一般控制步长，从以下例5.11也可看出外推法中的步长不须太小就可以得到较为精确的近似解。例5.11运用外推法计算在处的导数。 解：计算结果见下表(计算保留5位数字)精确值为e=2.718281828…。我们可以仿照外推法计算积分的算法5.1设计如下算法，算法中将外推表的最后一行按反序存入向量，因此算法结束时就是导数值。 算法5.2外推法计算导数[初始化]取初始步长，给定终止精度，，计算首个导数值；令，；[更新导数]；[递推]令，然后对循环3.1)3.2)[判断终止准则]如果，算法结束，为满足精度要求的近似导数值。，转1)。[算法结束]插值函数计算微分利用插值函数计算导数也是数值微分的一种方法。其思想很简单，就是对于给定的的函数值表，建立的插值函数，令，在某些情况下，用这种方法计算导数是可以接受的选择。考虑函数在插值区间上有个插值节点，其插值函数为，则对此等式两边求导令则有其中介于之间，显然用插值函数的导数近似表示在处的导数其误差由上式中的余项决