数值积分与数值微分

1、数值分析,第三章 数值积分与数值微分,1 引言,在工程问题和科学实验中，常常需要计算积分。例如：力学和电学中功和功率的计算，电流和电压的平均值和有效值的计算以及一些几何图形的面积、体积和弧长的计算等等。另外微分方程的求解也是以积分计算为基础的。,一 数值求积基本思想,利用积分中值定理： 即以底长b-a，高为 的矩形的面积 恰等于所求曲边梯形的面积i.,这样，只要对平均高度 提供一种算法，便获得一种数值积分方法。如果用两端点的“高度” 与 取算术平均作为平均高度 的近似值 ，这样导出积分近似公式：,对于这些情况，要计算积分的准确值都是十分困难的，这就要 求建立积分的近似计算方法。此外，积分的近似

2、算法又为其它一些数值方法，例如微分方程数值解、积分方程数值解等，提供了必要的基础。,问题是 的具体位置一般是不知道的， 因而难以准确算出 的值。 称 为区间a,b上的平均高度。,一般地，由定积分定义： 积 分为和式的极限,下面通过最简单的情形做进一步分析。如果把区间a,bn等份。节点：,定积分的基本分析步骤是四步：分割、近似、求和、取极限。分割就是把总体（整块梯形面积）分成若干分量（小曲边梯形面积），近似就是在每个分量中用容易计算的量去代表小曲边梯形的面积（这是用矩形面积近似曲边梯形的面积）。求和就是把分量加起来得到总近似值，最后取极限就得到积分的准确值。,以上分析看出，前三步比较容易，最后一

3、步的计算比较困难，但现在只是要求近似值，因而可以省掉求极限这一步，只要经过前三步就可求得积分 近似值，这就是建立数值积分方法的基本步骤。,再把每个小区间的分量相加就得梯形求积公式：,若在每个小区间 中用梯形面积近似曲边梯形面积 ，,由上述知，计算积分的问题归结为计算被积函数 f(x) 在 a,b 上某些点的值 ，这就使积分的计算变得很容易了。现在的问题是如何 提高精确度。从求近似积分的三步看来，前两步分割与近似是提高求积公式精确度的途径。第三步求和无须专门讨论。,先讨论不分割区间a,b 时，哪一种近似得到的求积公式精确度更高。事实上，对于不分割区间 而言，用零次多项式 近似f(x), 然后取积

4、分而得到的矩形公式,直观上看，它比左矩形公式（3）精确些，但总的说来，上述矩形法和梯形法精度都较低，往往不能满足实际计算的要求，因此需要建立精确度更高的求积公式。,分割主要是细分，从理论上讲，分得越细，精确度越高，但必然增加计算量，增大舍入误差。因此，怎样在一种合适的分割下，使用容易计算的分量去近似曲边梯形面积就成为提高求积公式精确度的主要问题。,同样，用一次多项式近似f (x), 取积分而得到梯形公式：,二 代数精度的概念,如果事先选定求积节点，如，以区间a,b的等距节点依次为节点，这时取m=n，求解上述线性方程组（）即可确定系数从而使求积公式至少有m=n次代数精度。具体示例在下面一节中介绍

5、。,由插值余项定理，对于插值型的求积公式其余项为：,三 插值型的求积公式,因此，对于插值型求积公式，由上面余项公式（9）可见，对于次数n 的多项式f(x)，其余项r(f)=0，因而这时求积公式至少具有n 次代数精度。,2 newton-cotes公式,一 cotes系数,由于是多项式积分，cotes 系数计算不会遇到实质性困难。,这就是抛物线公式，又称辛浦生 simpson 公式 。几何意义 就是用抛物线下的面积近似曲线f(x)下的面积。,这就是柯特斯公式（cotes ) 当n 较大时，例如 n=8 时，系数 中出现负数，而且有正有负会使舍入误差增大，数值稳定性较差，因此实际计算并不用 n较大

6、的 公式，而是将区间a,b 分割成若干个小区间，对每个或几个小区间应用n 较小的公式去计算。cotes系数表详见教材.,二 偶阶求积公式的代数精度,另一方面，直接求积分得：,易验证s=i, 即simpson 公式对次数不超过三次的多项式均能准确成立。,易验证此时si （ b a 时） 因此simpson公式实际上具有3次代数精度。一般的有： 定理2：当阶n 为偶数时， newton-cotes公式至少有n+1 次代数精度。,引进变换x=a+ht , 并注意到,证明：只要证明当n为偶数时， newton-cotes公式对,上面已说明 newton-cotes公式的余项为,若 n 为偶数，,三 几

7、种低阶求积公式的余项,当n=2 即simpson公式的余项。为此，先构造一个次数 的插值多项式h(x) ，使满足：,由于已证明simpson公式具有三次代数精度，因此， simpson公式对于这样构造的三次插值函数h(x) 是准确的：即,（由插值条件）,此即为simpson公式求得的f(x)的积分。,因此simpson公式余项为：,又可以证明：满足上面插值条件的多项式h(x)的插值余项为：,故积分余项为：,由上面截断误差易见，当积分区间a,b较大时，直接使用n-c公式的截断误差增大，积分近似值的精度难以保证。因此，在实际应用中，为了既提高结果的精度，又使算法简便且易在电子计算机上实现，往往采取

8、复合求积的方法，即：先将积分区间分成几个小区间，并在每个小区间上利用低阶newton-cotes公式计算积分近似值。然后对这些近似值求和，从而得所求积分的近似值。由此得到一些具有更大实用价值的数值求积公式，统称为复合求积公式。,四 复合求积法及其收敛性,可导出复合梯形公式：注意 (10),同理可得复化cotes公式,注：,（12）,（11）,定理：若f(x)在积分区间a,b上分别具有二阶、四阶和六阶连续导数，则复化梯形公式、复化simpson公式和复化cotes公式的余项分别为：其中a,b,证明：只证明复合梯形公式的余项，其余由学生类似可做。 由于f (x)在a,b连续，从而在每个小区间 上做

9、 积分 使用梯形公式时的截断误差为：（前面已证）,而 ，,故 #,f(x)在a,b上的连续性及介值定理知：有(a,b),使,又由上面（*）式及定积分定义有：,可见，当h充分小时有：,其余两类公式证明与此完全类似，略。,对于数值求积公式来说，收敛阶越高，近似值 收敛到其值 的速度就越快，在相近的计算工作量（顺便提一下，数 值求积计算工作量的大小，主要取决于计算函数值次数的多少） 下，有可能获得较准确的近似值。,例：利用复合newton-cotes公式计算 的近似值。,解：用两种方法求解。先将积分区间0，1八等分（分点及分 点处的函数值见下表），用复合梯形公式得：（h=1/8）,再将区间0，1四等

10、分， 利用复合simpson公式得：,两种方法都用到表中九个点上的函数值，计算工作量基 本相同。但所得计算结果与积分真值=3.14159265相比 较，复合simpson公式所得近似值 远比复合梯形公式 所得近似值 要精确。因此在实际计算中，较多的应用 复合simpson公式。 为了便于上机计算，常将复合simpson公式改写成：,(19),虽然,我们上面已得到newton-cotes低阶公式的近似值 的误差估计，也可根据精度要求用这些公式确定积分区间 的等份数，即确定步长 h ，但由于余项公式中含有被积函数 f(x) 的高阶导数，在具体计算时往往会遇到困难，因此，在 实际应用时，常常利用误差

11、的事后估计法来估计近似值的 误差，或确定步长 h。此方法的大致做法是：将积分区间逐步分半，每分一次就用同一复合求积公式算出相应的积分近似值，并用前后两次计算来判断误差的大小。,五 误差的事后估计与步长的自动选择,其原理和具体做法如下：,对复化梯形公式（10），由其余项公式（13）或（16）可 见，当f (x)在积分区间上变化不大或积分区间a,b的等分数n 较大（即步长h较小）时，若将a,b的等分数改为2n（即将步长 缩小到原来步长的一半），则新近似值 的余项约为原近似值 余项的1/4，即：（令 ）,先算出 和 ，若 （为计算结果的允许误差）， 则停止计算，并取 为积分近似值。否则，将区间再次分

12、半 后算出新的近似值 ，并检验不等式 是否成立， 直到得到满足精度要求的结果为止。,故将区间逐次分半进行计算的过程中，可以用前后两次计算结果 和 来估计误差和确定步长。具体做法是：,对于复合simpson公式（11），复合cotes公式（12）， 由它们的积分余项（14）、（15）或（17）、（18）可见， 当所涉及的高阶导数在积分区间上变化不大或积分区间的等份数n较大时，,有： 和,则 （21） 及 （22）,因此，也可以像使用复合梯形公式求积分近似值那样，在 将积分区间逐次分半进行计算的过程中，估计新近似值 和 的误差，并判断计算过程是否需要继续进行下去。上述过程很 容易在计算机上实现。,

13、先看复合梯形公式（10），计算 时，需计算n+1个点 （它们是积分区间a,b的n等分的分点）上的函数值，当 不 满足精度要求时，根据上面提供的方案，就应将各小区间分半， 计算出新近似值 。若仍用（10）计算 ，就需求出2n+1个点（a,b的2n等分点）上的函数值。,3 rombeng算法,一 复合梯形的递推化,上面介绍的步长变化的计算方案，虽然提供了估计误差 与选取步长的简便方法，但是还没有考虑到避免在同一节点上重复计算函数值的问题，故有进一步改进的余地。,而实际上，在这2n+1个分点中，包含有n+1个n分点，对应的函数值在计算 时已算出，现重新来计算这些点上的函数值，显然是极不合理的。,为避

14、免这种重复计算，我们来分析新近似值 与原有近 似值 之间的关系。由复合梯形公式（10）知：,注意到在2n分点 （k=1,2,2n-1)中， 当k取偶数时， 即为n分点，k为奇数时， 才是新增加的分 点，将新增加的分点处的函数值从求和记号中分离出来，就有:,由递推复化梯形公式（23）可见，在已计算出 基础上再计算 时，只要计算n个新分点上的函数值就行了，这与直接利用复合梯形公式（10）相比，计算工作量几乎节省一半。,例2：由（23）重新计算 的近似值，误差不 超过 。,(23),再将各个小区间二等分，出现两个新分点x=1/4, x=3/4,由 （23）得：,解：在积分区间逐次分半的过程中顺序计算

15、积分近似值 ， 并用是否满足 （= ）来判断计 算过程是否需要计算下去。,这样不断将各小区间二分 下去可利用递推公式（22） 依次算出 计算结果如右表，因为： 故 =3.14159202是满足 要求的近似值。,romberg算法是在积分区间逐次分半的过程中，对用复合 梯形法产生的近似值进行加权平均，以获得准确度较高的近似 值的一种方法，具有公式简练，使用方便，结果较可靠的优点。,上面介绍的递推公式（23）或（24），虽然具有结构简单， 易在计算机上实现等优点，但是由它产生的梯形序列 ，其 收敛速度却非常缓慢。如上例中用（23）计算 的近似值时，要一直算到 才获得误差不超过 的近似值。 因此，用

16、这种方法计算更复杂的高精度要求的近似值，显然是 费时、费力甚至是不可能的 。,二 romberg算法,如何提高收敛速度，节约计算工作量，自然是我们所最 关心的问题。,故（2）式成立,(由上节（10）,（由上节（11）,（4）,由近似式（22） 类似推导可得：,同理，由上节近似式（21）,类似推导可得：,（3）,即将simpson序列 按（3）作线性组合就可产生收敛,速度更快的新序列,例3: 利用公式(2).(3).(4)“加工”例2中得到的 近似值,解:依公式(2),(3),(4)按上图中计算顺序可得结果如下表: 其中k为二次分数.,这种加速的方法称为,romberg算法。其“加工”过程如下图

17、，,其中圆圈中号码表示计算顺序。, ,注意 教材中介绍的richardson外推法，为便于上机计算，引用记 号 来表示各近似值，其中k仍代表积分区间的二次分数，而 下标m则指出了近似值 所在序列的性质。如m=0在梯形序列中， m=1在simpson序列中，m=2在cotes序列中， 引入上面记号 后，romberg算法所用到的各个计算公式可统一化为：,可见，“加工”效果十分显著，而“加工”计算量因只需 做少量四则运算没有涉及到求函数值，故可忽略不计。,0 3,k,1 3.1 3.133333,2 3.131176 3.141596 3.1412118,3 3.138988 3.141593 3

18、.141594 3.141588,三 romberg算法计算 公式 的简化,由此可逐行构造 出一个三角形数表-称为t数表,.,这是因为： 时，梯形积分公式余项可写成：,此即为richardson外推法。其中 由（5）即知上面两 个极限成立。可见加速后余项数量级下降很快。,2、 可用二维数组来存放并参加运算，也可用一维数组。,四 几点说明,3、 对于积分限为无穷的积分 ，可利用变量代换化成有限区 间的积分然后再进行计算。例如：,4、若被积函数有奇异点（间断点）存在于积分区间内，则 可得积分 区间分成小部分，使间断点在子区间的端点处。 也可用变量代换法处理。,令,则,代入得：,一 gauss型求积

19、公式的引进,4 gauss型求积公式,但是，如果我们对公式（*）中系数 和节点 都不加 限制，那么就可以适当选取 使所得公式的代数精 度 m1。 事实上，若要求（*）对f（x）=1，x， 都准确成立，,证明： 设p(x)是任意次数不超过n的多项式。则 p(x)w(x) 的次 数不超过2n+1。因此，如果 是gauss点，则求 积分公式（1）对 p(x)w(x) 能准确成立。即有：,因此，可以想象对机械求积公式（1），只要适当选择2n+2个待定参数 ，是可以使其代表精度达 到 2n+1的。可以证明，gauss型求积公式是插值型的。,对于任意给定的次数不超过2n+1的多项式f（x），用 w（x）除

20、 f（x），记商为p(x)，余式为q（x），p（x） 与q（x）都是次数不超过n的多项式。 f(x)=p(x)w(x)+q(x), 由,而求积公式（1）是插值型的，代数精度至少为m=n，因此 对q（x）能准确成立：,但,故（2）成立。,由（2）得：,（2）、再利用gauss点确定求积系数,利用解方程组求gauss点 和权 以确定gauss公式，需 解2n+2个未知数的方程组，工作量太大。比较简便的做法是：,（1）、先利用 上的n+1次正交多项式确定gauss点,gausslegendre公式余项为：,再取 的两个零点 构造求积公式：,令它对f（x）=1，x均准确成立，即,从而得两点gaussl

21、egendre公式：,对于一般区间 上的积分，也可以利用教材中表格的数据,写出gauss型求积公式。其原理与方法是：先作变量替换，令,则将 上积分化为 上的积分：,记 ，则上式化为：,利用表中数据，对于给定的n=1，2，3，4，可以写出gauss型 公式：,即,代入（a）得：,其中系数 和节点 可查表得出，由变量替换公式,易见，由于求积公式（c）对变量t不高于2n+1,的多项式准确成立，从而求积公式（d）对自变量x的不高于2n+1,的多项式也准确成立，即（d）是gauss型求积公式。,例：利用四点gauss求积公式计算 的近似值。,解：由教材中表示gauss型求积公式（d）得：,其中a=0，b

22、=1，,= -0.86113631，,= -0.33998104，,将上述各数据代入上公式中有：,优点：在此例计算过程中，只涉及到 四个点上的函数值，可 见gauss型求积公式具有计算工作量小，所得近似值精确度 高的优点，是一种高精度的求积公式。,gauss型求积公式的缺点是：当n改变大小时，系数和节点几 乎都在改变。虽然可以通过其他资料查到较大n的系数和节点， 但应用时却十分不便。同时，余项却涉及到高阶导数（被积 函数的），要利用它们来控制精度也十分困难。,为克服这些缺点，在实际计算中较多地采用复合求积 的方法。例如，先把积分区间 分成m个等长的小 区间 然后，在每个小区间上使用同 一低阶（

23、如二点的，三点的， ）高斯型求积公式算 出积分近似值，再相加即将积分 的近似值：,其中， 由查表可得。同时在实际计算中，,还常用相邻两次计算结果 和 的关系式,相当于相对误差）即算出 后，观察 是否 成立，以判定是否终止计算。请同学们据此编程计算。,来控制运算（当 时，,证明：以 为节点构造次数 的多项式h(x) 使满足条件： 由第二章hermite插值知，h(x)为hermite插值多项式。 由于gauss公式具有2n+1次代数精度，则它对h(x)能准 确成立，即,三 gauss公式的余项,再注意到函数 在 上保号，应用积分中值定理即可得结论。,对比newton-cotes公式，gauss不

24、但具有高精度，而且是数值稳 定的。gauss公式之所以能够保证数值稳定性，是由于其求积函数 具有稳定性。,四 gauss公式的稳定性,，由于实际问题中，通常,不一定提供准确的数据 ，而只是给出含有误差的数据 （如：由于计算机字长的限制而产生的舍入误差）。因而实际 求得的积分值为：,那么，原始数据的误差对积分值的影响能否加以控制呢？ 上面已经证明gauss公式的求积系数 则,又由于gauss公式对f(x)=1准确成立,即,下面讨论求积过程的数值稳定性问题。,利用公式 计算积分时,可以仿照处理普通积分的方法讨论带权的积分。如，求积分公式,五 带权的gauss公式,确定了gauss点后，再利用gau

25、ss公式2n+1次代数精确度关系定 值得指出的是：运用正交多项式的零点构造gauss求积公式，这 种方法只针对某些特殊的权函数才有效。 构造gauss公式的一般方法是利用代数精度的概念用待定系数法 求解。 现举一例加以说明： 设要构造下列形式的gauss公式：,从而所构造的gauss求积公式为：,解：取f(x)=1.则对上求积公式，左端,六 综合例题,解此方程组得：,f(x)=x. 则 ，左端=,令左端=右端,左端=,再取,故，当取,时，求积公式具有3次代数精度。,例2 若用复化梯形公式计算积分,问积分区间要多少等分,才能保证有6位有效数字？,解：由复化梯形公式截断误差（13）式知,由于该积分

26、有一位整数，所以要求使近似积分有6位有效数字，只需取n满足：,（有效数字定义见第一章，表示不 超过某一位数字的一半）,即,即,因此，至少要将0，1区间212等分。若将同一问题改为复化 simpson公式，则由复化simpson公式截断误差（14）式同样可得：,由此可见，simpson公式复化型比复化梯形公式计算量少得多。,的近似值，要求误差,例3 用romberg求积法计算,解：此时积分限为a=0,b=1.而 (本例主要说明romberg过程）,如此继续算得：,由于,这个实例表明,romberg求积法计算过程不便于手工计算，但由于计算程序具有规律性，不必存储求积系数和节点，而且精度较高，因此适

27、合于在电子计算机上进行计算。,用romberg方法计算时，是把区间逐次分半的，因此有时称该法 为逐次分半加速法。,例4 构造三个节点的gauss-legendre求积公式，并给出余项估计式。 解：由于三次legendre多项式为：,其三个零点分别为：,令它对,准确成立,则三点gauss-legendre求积公式为：,余项为：,（n+1节点数）,例如，若要计算,的近似值，则由上积分公式得：,上述积分准确值为：,若利用三点simpson求积公式。则,可见在节点数目相同的情况下，gauss求积公式的精度是相当 高的。,例5 给出计算积分,的两点计算 公式，,使得对f(x)为三次多项式时精确成立。,解

28、：设,取 为二次多项式，对w（x） 上 式应精确成立：,显然 则,但,因而 即,不妨令 且 于是,令积分公式对f（x）=1，x准确成立 得：,解之得a=b=,故所求积分公式为：,显然，由上述过程知，积分公式对,精确成立。 可验证： 从而积分公式对,也准确成立。,令,（奇函数积分）,而,因此积分公式对,准确成立。,例6 求下列求积公式的代数精度,解：设,为任意实数,则,而,即积分公式对任意3次多项式准确成立。 又取,，则,而右端,故其代数精度为3次。 例7 建立下述形式的求积公式并确定它的代数精度：,解：由于有4个待定系数,一般应对三次多项式精确,成立。可取,得：,故所求积分公式为：,故所得积分

29、公式具有3次代数精度。,例8 导出下述形式的求积公式：,解：它有4个待定系数，应该对三次多项式准确成立。取,得：,易验证，它对 不精确成立，因此具有3次代数精度。,数值微分就是用离散方法近似地求出函数在某点的导数值。按照taylor展开原理可得,5 数值微分,其中h为一增量。上面几个公式是很实用的，下面我们再讨论一些常用方法。,3.5.1 插值型求导公式,设f（x）是定义在a，b上的函数，并给定区间a，b上的函数， 并给定区间a，b上的n+1个节点 出的函数值 这样,我们可以建立函数 的n次插值多项式 多项式的求导是容易的,称 (3.5.1) 为插值型求导公式。,应当指出，即使 和 的值相差不多，导数的近似值 与导数的值 仍然可能相差很大。因而在使用求导公式 （3.5.1）时，应注意误差的分析。,依据插值余项定理，求导公式（3.5.1）的余项为,式中,在上述余项公式中，由于 是 x 的未知函数，我们无法对右端的 第二项作出进一步的说明。因此，对于任意给出的点 x，求导公式 的余项是很难估计的。,然而，如果我们限定求节点上的导数值，那么有余项公式,（3.5.2）,下面我们考察节点处的导数值。为简化讨论，假定所给的节点是等距的， h是步长。,1.两点公式 当n=1时，由（3.5.2）得带余项的两点公式,（3.5.3）,（3.5.4）