证明隐式Euler方法稳定性.doc

第六章　数 值 积 分
6.1　数值积分基本概念
6.1.1　引言
在区间上求定积分
　　　　　　　　　　　　　　 　　(6.1.1)
是一个具有广泛应用的古典问题，从理论上讲，计算定积分可用Newton-Leibniz公式
　　　　　　　　　　　　　　 　　　(6.1.2)
其中F(x)是被积函数f(x)的原函数.但实际上有很多被积函数找不到用解析式子表达的原函数，例如等等，表面看它们并不复杂，但却无法求得F(x).此外，有的积分即使能找到F(x)表达式，但式子非常复杂，计算也很困难.还有的被积函数是列表函数，也无法用(6.1.2)的公式计算.而数值积分则只需计算f(x) 在节点xi(i=0，1,…,n)上的值，计算方便且适合于在计算机上机械地实现.
本章将介绍常用的数值积分公式及其误差估计、求积公式的代数精确度、收敛性和稳定性以及Romberg求积法与外推原理等.
6.1.2　插值求积公式
根据定积分定义，对及都有(极限存在)若不取极限，则积分I(f)可近似表示为
　　　　　　　　 　　(6.1.3)
这里称为求积节点，与f无关，称为求积系数，(6.1.3)称为机械求积公式.
为了得到形如(6.1.3)的求积公式，可在上用Lagrange插值多项式，则得

其中
　　　　　　　　　　 　　　(6.1.4)
这里求积系数由插值基函数积分得到，它与f(x)无关.如果求积公式(6.1.3)中的系数由(6.1.4)给出，则称(6.1.3)为插值求积公式.此时可由插值余项得到
　　　　　　　(6.1.5)
这里ξ∈，(6.1.5)称为插值求积公式余项.
当n=1时，，此时

由(6.1.4)可得

于是
　　　　　　　　　 　　　(6.1.6)
称为梯形公式.从几何上看它是梯形AbB(见图6-1)的面积近似曲线y=f(x)下的曲边梯形面积，公式(6.1.6)的余项为
　　　　　　　　　　　　(6.1.7)
6.1.3　求积公式的代数精确度
当被积函数即f为次数不超过n的代数多项式时，，故由(6.1.5)有，它表明插值求积公式(6.1.3)精确成立.对一般机械求积公式(6.1.3)，同样可以根据公式是否对m次多项式精确成立作为确定公式(6.1.3)中系数及节点的一种方法.在此先给出定义.
定义1.1　一个求积公式(6.1.3)若对精确成立，而对不精确成立，则称求积公式(6.1.3)具有m次代数精确度.
根据定义，当时公式(6.1.3)精确成立，故有等式
　　　　　　　　(6.1.8)
而　　　　　
(6.1.8)是关于系数及节点的方程组，当节点给定时，(6.1.8)取m=n就是关于系数的线性方程组，求此方程组就可求得求积系数.
例如n=1，取，求积公式为

在(6.1.8)中令m=1，可得

解得　　　　　　　　　　　　　　
它就是梯形公式(6.1.6)的系数，它与用公式(6.1.4)算出的结果完全一样.对梯形公式(6.1.6)，当时

故求积公式(6.1.6)的代数精确度为一次.
对于具有(n+1)个节点的插值求积公式(6.1.3)，当时，故公式精确成立，它至少有n次代数精确度.反之，若求积公式(6.1.3)至少有n次代数精确度，则它是插值求积公式，即(6.1.3)的求积系数一定可用(6.1.4)求出.实际上，此时对求积公式(6.1.3)精确成立，若取f(x)为插值基函数，即

由(6.1.3)精确成立，可得

这就是(6.1.4)得到的插值求积公式系数.
定理1.1　求积公式(6.1.3)是插值求积公式的充分必要条件是(6.1.3)至少具有n次代数精确度.
定理表明直接利用代数精确度概念，由(6.1.8)可求得插值求积公式.更一般地，含有被积函数的导数的求积公式也同样可用代数精确度定义建立.如下例所示.
例6.1　求积公式，已知其余项表达式为.试确定系数及，使该求积公式具有尽可能高的代数精确度，并给出代数精确度的次数及求积公式余项.
解　本题虽用到的值，但仍可用代数精确度定义确定参数及.令，分别代入求积公式.令公式两端相等，则
　　　　当
　　　　当
　　　　当
解得，于是有

再令，此时，而上式右端为，两端不等，则求积公式对不精确成立，故它的代数精确度为二次.
为求余项可将代入求积公式

当，代入上式得
，即
所以余项.
6.1.4 　求积公式的收敛性与稳定性
定义1.2　若[，则称求积公式(6.1.3)是收敛的.
定义中n→∞包含了通常都要求用于计算积分(6.1.1)的求积公式(6.1.3)是收敛的.本章后面给出的求积公式都必须先证明其收敛性.
稳定性是研究计算和式

当有误差时，的误差是否增长.现设，误差为.
定义1.3　对任给，只要，就有

则称求积公式(6.1.3)是(数值)稳定的.
定义表明只要被积函数f(x)的误差充分小，积分和式的误差限就可任意小，则(6.1.3)就是数值稳定的.
定理1.2　若求积公式(6.1.3)的系数则求积公式(6.1.3)是稳定的.
证明　由于，故有

于是对，只要，就有

故求积公式(6.1.3)是稳定的.
讲解：
数值积分就是将求积分转化为求和（即6.1.3式）这样不管被积函数多么复杂，它都能在计算机上机械实现。把（6.1.3）式称为机械求积公式，为求积节点，为求积系数，建立求积公式有两种途径，一是利用 的插值多项式积分得到，二是根据代数精确度概念，通过解方程得到及。特别当节点给定时，方程（6.1.8）是关于的线性方程组，它是容易求解的。
定理1.1给出了插值求积公式与代数精确度之间的关系。
求积公式收敛性简单说就是当和式收敛于积分值。而稳定性是研究计算和式 的误差积累，即当有误差时，只要误差充分小则和式误差也任意小，这就是稳定的。定理1.2表明只要求积公式（6.1.3）的系数，则求积公式就是稳定的。
6.2　梯形公式与Simpson求积公式
6.2.1　Newton-Cotes公式与Simpson公式
在插值求积公式中，若求积节点，此时可将求积公式写成
　　　　　　　　 　　　(6.2.1)
称为Newton-Cotes求积公式，其中系数称作Cotes系数.按(6.1.4)式引入变换，则有
　　　　　　(6.2.2)
由于是多项式积分，计算不会有困难.例如n=1时，，.这时求积公式就是我们熟悉的梯形公式(6.1.6).n=1～8时的系数见表6-1.从表中看到n=8时出现负数，稳定性没保证,所以一般只用n≤4的公式.
表　6-1

当n=2时，由(6.2.2)可得，求积公式为
　　　　　　　　(6.2.3)
称为Simpson求积公式.
　　对梯形公式(6.1.6)，已知它的代数精确度为一次，且它的余项已由(6.1.7)给出，若记，则

由于，在上不变号，由积分中值定理得知，使
　　　　　　　(6.2.4)
这就是梯形公式(6.1.6)的截断误差.
　　对Simpson公式(6.2.3)可以证明它的代数精确度是三次，根据定理1.1，显然(6.2.3)对精确成立，再对，左端为，右端为

故(6.2.3)对也精确成立.而对，公式(6.2.3)不精确成立.故求积公式(6.2.3)的代数精确度是三次，即(6.2.3)对任何精确成立.若令P(x)满足条件

这里,于是由(6.2.3)有

根据上章例5.6中(5.5.12)式有

上式两边积分，并记，则得

由于在区间上(不变号，故由积分中值定理得

于是有
　　　　　　　　　　　　(6.2.5)
这就是Simpson求积公式(6.2.3)的余项，即截断误差.
讲解：
当求积节点(k=0,1,…n)为等距节点，直接由插值求积得到的求积公式就转为Newton-Cotes求积公式。但使用时通常只用n=1，2，4三种情况，它们分别称为梯形公式，Simpson公式和 Cotes公式，梯形公式

它的局部截断误差直接由插值余项积分得到，由（6.2.4）给出，即

对于n=2的情形给出的Simpson求积公式（6.2.3），它具有3次代数精确，也就是它对任何次数不超3次的多项式精确成立，为使求积公式余项能反映这一性质，我们不用2次的Lagrange插值近似被积函数f(x)，而用带导数的3点3次Hermite插值，即造，满足条件：

于是由Hermite插值余项表达式，可得

对上式两端从到b积分，并利用积分中值定理就可得到Smpson求积公式的余项（6.25）对于n=4的Newton-Cotes公式根据（6.2.2）可算出，　　　，于是有

这里，余项
6.2.2　复合梯形公式与复合Simpson公式
直接用梯形公式(6.1.6)及Simpson公式(6.2.3)计算积分I(f)误差较大，为达到精度要求，通常可将分为n个小区间，在小区间上应用梯形公式及Simpson公式即可达到要求.为此取分点,在每个小区间上用梯形公式(6.1.6)，则得

或
　　　　　　　　 　　　(6.2.6)
称为复合梯形公式.根据定积分定义可知

故复合梯形公式(6.2.6)是收敛的，且(6.2.6)的求积系数，故它也是稳定的.
(6.2.6)的截断误差可由(6.2.4)得到

若，根据连续函数性质，,使

于是得

它表明复合梯形公式(6.2.6)的截断误差阶为
如果在每个小区间上使用Simpson公式(6.2.3)，则得
(6.2.8)
称为复合Simpson公式，它的余项由(6.2.5)可得

即　　　　　 (6.2.9)
它表明.此外，还可证明

故复合Simpson公式(6.2.8)是收敛的，并且，故公式也是稳定的.
例6.2　用n=8的复合梯形公式及n=4的复合Simpson公式，计算积分，并估计误差.
　　解　只要将区间［0，1］分为8等分，用公式(6.2.6)时取n=8，h=0.125，对复合Simpson公式取n=4,h=0.25.计算各分点的函数值.

由公式(6.2.6)及(6.2.8)得

为了估计误差，要求的高阶导数，由于

所以

故

由(6.2.7)得

对Simpson公式，由(6.2.9)得

例6.3 计算积分，若用复合梯形公式，问区间［0,1］应分多少等分才能使截断误差不超过?若改用复合Simpson公式，要达到同样精确度，区间［0，1］应分多少等分?
解 本题只要根据及的余项表达式(6.2.7)及(6.2.9)即可求出其截断误差应满足的精度.由于，对复合梯形公式

即.取n=213,即将区间［0,1］分为213等分时，用复合梯形公式计算误差不超过.而用复合Simpson公式，则要求

即.取n=4，即将区间分为8等分，用n=4的复合Simpson公式即可达到精确度.
讲解：
由于插值求积公式只能用n=1,2,4,用它算积分往往达不到精度要求，为提高求积精度我们可以采用将求积区间分为n个等距的小区间在每个小区间上应用梯形公式或Simpson公式，这就得到复合梯形公式及复合Simpson公式。它们的余项由每个小区间积分余项相加即可得到，结果分别由（6.2.7）及（6.2.9）给出。这两个求积公式是收敛和稳定的。如何利用公式计算积分近似值并估计误差可见例6.2及例6.3。
6.3　外推原理与Romberg求积
6.3.1　复合梯形公式递推化与节点加密
在计算机上用等距节点求积公式时，若精度不够可以逐步加密节点.设将区间分为n等分，节点，在区间上梯形公式为

若节点加密一倍，区间长为，记中点为在同一区间上的复合梯形公式为

于是
　　　 (6.3.1)
它表明是在的基础上再加新节点的函数值之和乘新区间长，而不必用(6.2.6)重新计算，这时有误差估计式

若，则得
　　　　　　　　　　　　(6.3.2)
它表明用，其误差近似.这也是在计算机上估计梯形公式误差的近似表达式.若(给定精度)，则.
若在区间［a,b］中做2n等分时，在上用Simpson公式计算，则由(6.2.8)可知

它恰好是(6.3.2)中I(f)的近似值，即

它表明用(6.3.2)计算I(f)，其精度已由提高到如果再将区间分半，使分为4个小区间，长度为，则可由(6.3.1)计算出及，利用复合公式余项(6.2.9)得

如果，则有
　　　　　 　(6.3.3)
从而有复合Simpson公式的误差估计

如果用(6.3.3)近似，即
　　　　　　　　　　　　　 　　(6.3.4)
则精度可达到.类似做法还可继续下去.这样对区间逐次分半，利用公式(6.3.1)逐次递推.再由(6.3.2)，(6.3.3)逐次构造出精度愈来愈高的计算积分I(f)的公式，这就是Romberg求积的基本思想.
6.3.2　外推法与Romberg求积公式
仍从梯形公式出发，区间中的节点如前所述.此时复合梯形公式可表示为

当分为2n等分，区间长变为时，记，由于

此处将看作h的函数，将按的幂展开，可得到以下结果.
定理3.1　设f在上的各阶导数存在，则复合梯形公式可展成
　　　　　　　　　(6.3.5)
其中为不依赖h的常数.
定理证明见［2］.由(6.3.5)，显然有.在(6.3.5)中若用代替则得
　　 　　　　(6.3.6)
用4乘以(6.3.6)减去(6.3.5)除以3，则得

若记
　　　　　　　　　　　　　(6.3.7)
显然
　　　　　　　　　　　　　　(6.3.8)
具有精度.实际上就是复合Simpson公式中的为提高精度，可由及中消去，得

用精度为，如此逐次做下去，可得到
　　　　　　 　　　　　(6.3.9)
用时，精度为，这种将步长h逐次减半，使逼近，以便精度逐次提高的方法称为外推法，它对于可展成h的幂级数的计算公式的加速收敛是很有效的，这里只将外推法用于计算积分.
下面若用表示将区间二分k次得到的复合梯形公式，此时分为等分，步长，当k=0,1,…逐次得到即为等分的复合梯形公式，加速一次得序列即为Simpson公式序列.加速m次则得，由(6.3.9)可将它表示为
　　　　　　　　　(6.3.10)
称为Romberg求积公式.计算从k=0，即h=b-出发记，逐次二分得到T表(见表6-2).
当k增加时，先由(6.3.1)根据算出，再由(6.3.10)对m=1,2,…,k计算.当f充分光滑时可证明
　　　　　　　　　　(T表任一列)
　　　　　　　　　　　(T表对角线)
计算到(精度要求)为止.
例6.4 用Romberg求积公式求的近似值，使其具有6位有效数字.
解 本题直接用梯形递推公式(6.3.1)及Romberg求积公式(6.3.10)，按T表依次计算

其余计算结果见T表.

由于，故计算停止，I≈0.946 083 1即为所求.
例6.5 证明等式
试依据的值，用外推算法求π的近似值.
解 本题可利用Taylor展开式用外推原理求π的近似值.可令.由Taylor公式展开得

若记，则其误差为.
由外推法
　　　　，其误差为
　　　　，其误差为
根据以上公式计算结果如下表所示.

π=3.141 58即为所求.
讲解：
　　在等距节点的情况下，通过对求积区间的逐次分半，由梯形公式出可逐次提高求积公式精度，这就是Romberg求积的基本思路，由于梯形公式余项只有精度，即，但当节点加密时可组合成其精度达到，如果再由与组合成则可使误差精度达到，于是依赖于x，若在上各阶导数存在，将展开，可将展成的幂级数形式，即（6.3.5）表达式，记的计算精度，可利用外推原理逐次消去式（6.3.5）右端只要将步长h逐次分半，利用及组合消去，重复同一过程最后可得到递推公式（6.3.9），此时.说明用其误差阶为，这里表示m次加速。计算时用序列表示区间分半次数，即具体计算公式为（6.3.10），就是Romberg求积方法。Romberg求积是外推原理在数值积分中的应用，外推法也用于数值微分，微分议程数值解，实际上只要计算的问题可类似（6.3.5）展成h的幂级数可利用外推法，如例6.5计算的值，利用的Taylor展开式就可利用外推法。
6.4 Gauss型求积公式
6.4.1 最高代数精确度求积公式
考虑带权积分

其中为权函数，它的求积公式为
　　　　　　　　　　(6.4.1)
其中为求积系数，不依赖于f,为求积节点，在式(6.4.1)中作为待定参数，可选择使(6.4.1)对精确成立，从而得到关于,的(2n+2)个参数的非线性方程组，由式(6.4.1)得
　　　　　　　　　　　(6.4.2)
例如，当n=0，ρ(x)=1时，求积公式为
当=1，得，得，于是，可得求积公式
　　　　　　　　　　　　　 　　(6.4.3)
称为中点求积公式，它的代数精确度为一次.
例6.6 当n=1时求积分公式的系数及节点，使它具有最高代数精确度.
解 由代数精确度定义，公式对精确成立，由(6.4.2)得

它是关于与的非线性方程组，由前两式有

由后两式得

故有

化简得

令，上式可改写为

解得，从而求得

从例题看到直接解方程组(6.4.2)计算太复杂，n≥2时一般都不易求解.但若先确定求积节点，则由(6.4.2)求出系数就容易了.下面先证明求积公式(6.4.1)的代数精确度最高为(2n+1)次.若令，则，而

说明(6.4.1)对(2n+2)次多项式不精确成立，故它的最高代数精确度为(2n+1)次.具有(2n+1)次代数精确度的求积公式称为Gauss型求积公式，相应的求积节点称为Gauss点.
定理4.1 插值求积公式(6.4.1)的节点是Gauss点的充分必要条件是区间上以这组节点为零点的多项式

与任何次数不超过n的多项式带权正交，即
(6.4.4)
证明略，见［4］.
根据此定理可知，Gauss型求积公式的节点就是在上带权正交多项式的零点，这就避免了根据代数精确度定义求解非线性方程组的困难.在给出求积节点后，求积系数就可直接由解(6.4.2)的前(n+1)个方程组得到.而公式(6.4.1)的余项可通过f(x)的Hermite插值多项式得到，设，满足插值条件

于是有

两端乘权函数，并从到b积分，则得
　　　　　　　　　　　　(6.4.5)
其中
　　　　　　　　　　　 　(6.4.6)
(6.4.5)右端第一项是因为(6.4.1)对任何(2n+1)次多项式精确成立，且，则可得到，而可利用积分中值定理得到的.下面还可证明Gauss型求积公式(6.4.1)的稳定性和收敛性.
4.2 Gauss型求积公式(6.4.1)的系数.
证明 由于(6.4.1)对任何不大于(2n+1)次的多项式精确成立，若取，其中为n次插值基函数

于是，故有

证毕.
再由定理1.2得如下推论.
推论 求积公式(6.4.1)是稳定的.
定理4.3 设，则Gauss型求积公式(6.4.1)是收敛的，即

定理证明可见［4］.
6.4.2 Gauss-Legendre求积公式
若，区间为［-1，1］的求积公式
　　　　　　　　　　　　 　　　(6.4.7)
其中节点是Legendre多项式

的零点，则(6.4.7)称为Gauss-Legendre求积公式，简称Gauss求积公式.其系数

这里是最高项系数为1的Legendre多项式.余项可由(6.4.6)得到
　　　　　　　　　 　(6.4.8)
例如当n=1时，，则，，.它比Simpson公式(三点)的余项的系数绝对值小，Gauss求积公式(6.4.7)的节点与系数可见表6-3.
表6-3

例6.7 用四点(n=3)的Gauss求积公式计算

　　解 先将区间变换为［-1，1］，令

其中

(精确解=0.467 401…)
这结果与用n=8的Romberg求积相当.
6.4.3　Gauss-Chebyshev求积公式
区间为［-1，1］，权函数的Gauss型求积公式，其节点是Chebyshev多项式的零点，即，而，于是得到
　　　　　　　　 　　　　(6.4.9)
称为Gauss-Chebyshev求积公式，公式的余项为
　　　　　　　　　　　　(6.4.10)
这种求积公式可用于计算奇异积分.
例6.8　用三点和四点Gauss-Chebyshev求积公式计算积分，并估计误差.
解　这里，由Gauss-Chebyshev求积公式(6.4.9)可得

当n=2时，，求得

代入上式得

估计误差可用余项表达式(6.4.10)，因，故

当n=3时，，求得

误差　　　　　
讲解：
Gauss型求积公式是上带权的求积公式（6.4.1），它具有最高代数精确度2n+1，实际上由于求积系数及节点都是待定系数，它共有2n+2个，可使（6.4.1）对任何2n+1次多项式精确成立。具有2n+1次代数精确度的求积公式节点就是Gauss点，它可由定理4.1确定，实际上它就是在上带权正交多项式的零点。得到求积节点以后，同样可利用（6.4.1）对精确成立，得到关于的线性方程组

解此方程组得到的求积公式系数，它是稳定的，也是收敛的，具有较高的精度。通常使用的具体公式是Gauss-Legendre求积公式（简称Gauss求积公式），它是区间为，权函数为 的公式（6.4.7），其节点系数均可直接由表6－3得到，余项由表达式（6.4.8）给出，当n=1时可得，比n=2（三点）simpson公式好，当n=2时可得比n=4（五点）的Cotes公式好，而计算量却减少。另一个Gauss型求积公式时Gauss-Chebyshev求积公式，它由（6.4.9）给出，它除了精度高，还可计算反常积分，如例6.8。
【本章小结】
1.建立求积公式并求其余项表达式是本章重点，通常有两种方法，一是利用插值多项式积分，二是根据求积公式精确度定义。
（1）利用被积函数的插值多项式建立求积公式，在求积节点给定式可直接由公式（6.1.4）求得系数　，对只有一个节点的求根公式可利用函数的Taylor展开得到，其余项表达式也可利用插值（含带导数插值）的余项积分得到。
（2）利用求积公式代数精确度定义可得到含待定参数（系数或节点）的方程组，直接解方程组即可求得系数和节点，但对节点未知的情形方程组为非线性，只能对m较小情形才能计算，一般都是在节点求出后，方程组变成(k=0,1,…,n)的线性方程组，才能较方便求出.至于求公式的余项表达式在已知求积公式代数精确度为m,则可得到求积公式余项表达式为
，由于余项为，为确定系数K，可令，而，由此可得 如例6.1所示。
2.用已知求积公式计算积分并估计误差是本章最常见和必须掌握的内容。
（1）等距节点求积公式最基本的是梯形公式，Simpson公式及其余项和相应的复合求积公式及其余项，要能用这些公式计算积分并估计误差，也能根据误差要求确定区间的等分数。
（2）对Romberg求积方法，要掌握提高精确度的外推原理及利用外推原理计算积分的Romberg求积法。
（3）Gauss型积分是具有最高代数精确度的求积。共求积节点是相应正交多项式的零点，这类积分节点系数都可直接查表，只要能用Gauss公式及Gauss-chebyshev求积公式计算积分即可。在区间不是时计算要通过区间变换，将变换为.
【本章重点】
1． Euler法，后退Euler法，改进Euler法，梯形法的基本公式和构造，并用它计算初值问题数值解。
2． Runge-Kutta法基本思想，二阶R－K法推导及经典四阶R－K法的使用。
3． 单步法的局部截断误差定义及局部截断误差主项与方法阶的确定。
4． 单步法的绝对稳定性及常用方法绝对稳定区间的确定。
5． 线性多步法的构造及局部截断误差的推导，主要是用Taylor展开方法，并能用Adams公式，Hamming公式计算初值问题的解。　　
【课前思考】
1. 怎样导出Euler法，后退Euler法，改进Euler法和梯形法的公式？它们中那些是显式公式？那些是隐式公式？
2. 何谓方法的局部截断误差？方法的阶？试求Euler法和梯形法的局部截断误差和阶。
3. 什么是单步法的绝对稳定性？试给出Euler法，改进Euler法，四阶R－K方法及梯形法的绝对稳定区间。
4. 怎样推导二阶R－K方法的公式。
5. 什么是线性多步法？它有什么优缺点。怎么用Taylor展开推导Adams显式和隐式公式并求出它的局部截断误差主部。（以四阶为例证明）
30