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

第7章 数值微分和数值积分

7.1 数值微分

7.1.1 差商与数值微分

当函数算

的导数

是以离散点列给出时，当函数的表达式过于复杂时，常用数值微分近似计

。在微积分中，导数表示函数在某点上的瞬时变化率，它是平均变化

率的极限；在几何上可解释为曲线的斜率；在物理上可解释为物体变化的速率。

以下是导数的三种定义形式：

（7.1）

在微积分中，用差商的极限定义导数；在数值计算中返璞归真，导数取用差商（平均变化率）作为其近似值。

最简单的计算数值微分的方法是用函数的差商近似函数的导数，即取极限的近似值。下面是与式（7.1）相应的三种差商形式的数值微分公式以及相应的截断误差。 向前差商

用向前差商（平均变化率）近似导数有：

（7.2）

其中定义。

的位置在的前面，因此称为向前差商。同理可得向后差商、中心差商的

由泰勒展开

1 / 49

得向前差商的截断误差：

向后差商

用向后差商近似导数有：

（7.3）

与计算向前差商的方法类似，由泰勒展开得向后差商的截断误差：

中心差商

用中心差商（平均变化率）近似导数有：

由泰勒展开

（7.4）

得中心差商的截断误差：

差商的几何意义

微积分中的极限定义，表示在处切线的斜

率，即图7.1中直线的斜率；差商

的斜率，是一条过

表示过和

两点直线 的割线。可见数值微分是用近似值内

接弦的斜率代替准确值切线的斜率。

图7.1 微商与差商示意图

例7.1 给出下列数据，计算，

－0.02 0.04 0.06 5.06 0.8 0.10 5.07 5.065 5.05 5.055

解：（5.07－5.06）/（0.04－0.02）= 0.5

（5.05－5.07）/（0.08－0.04）= －0.5

（5.05－5.055）/（0.08－0.10）= 0.25

（(0.10) －(0.06)）/（0.10－0.06）= 18.75

设定最佳步长

在计算数值导数时，它的误差由截断误差和舍入差两部分组成。用差商或插值公式近似导数产生截断误差，由原始值的逼近值的角度看，

的数值近似产生舍入误差。在差商计算中，从截断误差

会带来较大的舍入误差。怎

越小，则误差也越小；但是太小的

样选择最佳步长，使截断误差与舍入误差之和最小呢？

一般对计算导数的近似公式进行分析可得到误差的表示式，以中心差商为例，截断误差不超过

而舍入误差可用量差为

估计（证明略），其中是函数的原始值的绝对误差限，总误

当时，总误差达到最小值，即

（*）

可以看到用误差的表达式确定步长，难度较大，难以实际操作。

通常用事后估计方法选取步长，例如，记为步长等于的差商计算

公式，给定误差界，当时，就是合适的步长。

例7.2 对函数佳步长。 解：

0.10 0.09 0.08 0.07 0.06 ，取不同的步长计算，观察误差变化规律，从而确定最

误差 误差 3.1590 －0.0008 3.1588 －0.0006 3.1583 －0.0001 3.1575 －0.0007 3.1550 －0.0032 3.1630 －0.0048 0.05 3.1622 －0.0040 0.04 3.1613 －0.0031 0.03 3.1607 －0.0025 0.02 3.1600 －0.0018 0.01 表中数据显示，当步长从0.10减少到0.03时，数值微分误差的绝对值从0.0048减少到0.0001，而随着

的进一步减少，误差的绝对值又有所反弹，表明当步长小于0.03时，

舍入误差起了主要作用。

在实际计算中是无法得到误差的准确数值的，这时以定步长，本例中取= 0.04。

最小为标准确

7.1.2 插值型数值微分

对于给定的

的导数。

的函数表，建立插值函数，用插值函数的导数近似函数

设为上的节点，给定，以

的各阶导数近似

，以

的相应阶的导数，即

为

插值点构造插值多项式

当时，

误差项为：

（7.5）

例7.3 给定

。

，并有，计算

解：作过的插值多项式：

将代入 得三点端点公式和三点中点公式：

利用泰勒（Taylor）展开进行比较和分析，可得三点公式的截断误差是 类似地，可得到五点中点公式和五点端点公式：

。

7.1.3 样条插值数值微分

把离散点按大小排列成

的样条函数

：

，用关系式构造插值点

当则当时，可用计算导数。

7.2 数值积分

在微积分中用牛顿－莱布尼兹（Newton-Leibniz）公式计算连续函数

的定积分：

但是，当被积函数是以点列的形式给出时，当被积函数

的原函数难以得到时，例如，则无法用牛顿－莱布尼兹积分公式计算。有时当被积函数的原函数过于复杂时，也不宜套用积分公式计算积分，而应采用数值积分公式计算定积分。

在微积分中，定积分是黎曼（Rimann）和的极限，它是分割小区间长度趋于零时的极限，即

在数值积分公式中，只能用有限项的和近似上面的极限，通常由函数在离散点函数值的线性组合形式给出。

记

表示近似积分值，

，在本章中，用

称为积分节点，称为积分系数。确定

表示精确积分值，用

中积分系数

的

过程就是构造数值积分公式的过程。

怎样判断数值积分公式的效果？代数精度是衡量数值积分公式优劣的重要标准之一。

定义7.1 (代数精度)记上以为积分节点的数值积分公式

为

若

具有

满足

阶代数精度。

而，则称

由此可知当具有阶代数精度时，对任意的阶多项式都有 。

7.2.1 插值型数值积分

对给定的被积函数在

，以

上的点列

，即

作拉格朗日插值多项式

近似计算

记，则有

数值积分误差，也就是对插值误差的积分值

或

对一般的函数而有

，但若是一个不高于次的多项式，由于，

。因此，阶插值多项式型式的数值积分公式至少有阶代数精度。

例7.4 建立上节点为的数值积分公式。

解：由得

得以数值积分公式：

7.2.2 牛顿-柯特斯（Newton-Cote’s）积分

把积分区间分成等分，记步长为，取等分点，取

作为数值积分节点，构造拉格朗日插值多项式

由此得到的数值积分称为牛顿—柯特斯积分。下面可以看到，牛顿－柯特斯积分系数和积分节点以及积分区间无直接关系，系数固定而易于计算。 梯形积分

以和为插值节点构造线性函数，有

那么，

提取公因子后，得到牛顿－柯特斯积分的组合系数：，，它们

已与积分区间没有任何关系了。

记 （7.6）

称为梯形积分公式。它的几何意义是用梯形面积近似代替积分值（图7.2）。

图7.2 梯形积分

怎样确定梯形积分公式的代数精度？

我们可以取验证

取时，有

即：

取时，有

即：

取时，有

得梯形求积公式具有一阶代数精度。 梯形积分公式的误差：

由

得

因为在上不变号，由积分中值定理得到梯形求积公式的截断误差：

辛普森（Simpson）积分

（7.7）

对区间作二等分，记

和

。以

为插值节点构造二次插值函数

，那

么，有

计算得到积分组合系数：，，。

（7.8）

S (f)称为辛普森或抛物线积分公式。它的几何意义是用过三点的抛物线面积近似代替积分的曲边面积（图7.3）。

图7.3 抛物线积分面积

分别将代入到和中，可以得到

表明辛普森公式对于次数不超过三次的多项式准确成立，此可设

一个三次多项式满足条件：

具有三阶代数精度。因

计算得到误差为：

，

于是有

但

故辛普森求积公式的截断误差：

牛顿-柯特斯积分系数

（7.9）

以

等分区间，取等分点为积分节点，

为插值节点构造插值函数

。

，其中。

其中

令，代入上式得

（7.10）

这里称为牛顿-柯特斯系数

可见在取等距节点时，积分系数与积分节点和积分区间无直接关系，只与插值的节

点总数有关，而在例7.3中的积分系数是待定系数，这就简化了数值积分公式，而不必对每一组插值节点xi都要计算一组相应的积分系数ai。在公式（7.10）中取=1，可算出梯形积分系数；取=2，可算出辛普森积分系数。在表7.1中列出从1到6的牛顿-柯特斯系数。从表中可以看出牛顿-柯特斯系数具有对称性。

表7.1

1 2 3 4 5 6 7.2.3 求积公式的收敛性与稳定性

定义7.2 若的。

，则称求积公式是收敛

定义中的的。

包含了，通常都要求计算积分的求积公式是收敛

稳定性是研究计算公式增长，现设

，误差记为

当有误差时，。

的误差是否

定义7.3 对任给，只要，就有

则称求积公式是稳定的。

定理7.1 若求积公式稳定的。

的系数，则求积公式是

证明 由于,，故有

于是对，只要，就有

故求积公式是稳定的。

7.3 复化数值积分

由插值的龙格现象可知，高阶牛顿-柯特斯积分不能保证等距数值积分系列的收敛性，同时可证（略）高阶牛顿-柯特斯积分的计算是不稳定的。因此，实际计算中常用低阶复化梯形等积分公式。

7.3.1 复化梯形积分

把积分区间分割成若干小区间，在每个小区间上用梯形积分公式，再将这些小

区间上的数值积分累加起来，称为复化梯形公式。复化梯形公式用若干个小梯形面积逼近积分

比用一个大梯形公式效果显然更好，如图7.4所示。这种作法使我们想起定积分

定义，即它为被积函数无限分割的代数和。这也正是计算定积分最朴素的算法。

图7.4 复化梯形公式积分视图

复化梯形积分计算公式

对作等距分割，有，，于是

在上，，则有

记等分的复化梯形公式为，有

复化梯形公式截断误差

（7.11）

由，根据均值定理，当时，存在

，有，于是

（7.12）

由此看到复化梯形公式的截断误差按照只要

在

或者的速度下降，事实上，可以证明，

上有界并黎曼可积，当分点无限增多时，复化梯形公式收敛到积分。

记 ，则有

对于任给的误差控制小量，有

或

就有，式中表示取其最大整数。

7.3.2 复化辛普森积分

把积分区间分成偶数等分总数。

，记，其中是节点总数，是积分子区间的

记，，在每个区间。

上用辛普森数值积分公

式计算，则得到复化辛普森公式，记为 复化辛普森积分计算公式

而，称

（7.13）

为复化辛普森积分公式，它是在上采用辛普森积分公式叠加而得。下

，在每个积分区间

面用图7.5显示复化辛普森积分计算公式中节点与系数的关系，取

上提出因子后，三个节点的系数分别是1，4，1；将4个积分区间的系数按节点的位置

累加，可以清楚地看到，首尾节点的系数是1，奇数点的系数是4，偶数点的系数是2。

1 4 1 4 1 1 4 1 1

1 1 4 2 4 2 4 2 4 4 1 1 图7.5 复化辛普森积分系数

复化辛普森公式的截断误差

设，在上的误差为

因此，

即 （7.14）

与复化梯形公式类似，误差的截断误差按照

在

或者的速度下降。可以证明，只要

上有界并黎曼可积，当分点无限增多时，复化辛普森公式收敛到积分

。

记，则有

对任给的误差控制小量，只要

或

就有。

例7.5 求

公式的分点应取多少？

，计算中要求有5位有效数字。用复化梯形和复化辛普森求积

解：

由复化梯形误差公式得到：

计算出，复化梯形公式至少要在.00等分n = 68。

由复化辛普森误差公式，有

在复化辛普森公式中取 或。

7.3.3 复化积分的自动控制误差算法

复化积分的误差公式表明，截断误差随分点的增大而减小，对于给定的误差量，用估计函数导数的界的方法可计算出。用误差公式计算满足精度的分点数，像是在做一道计算导数

上界的微积分习题（如例7.5所示）。但是在实际运算中，一般难以估计

出函数的各阶导数界，也就无法确定分点数。在计算中常用误差的事后估计方法，即用

估计误差

。

T2n (f )的计算公式

对定积分，取分点，计算得

取分点

，计算得

这里，。可以看到，的值是与新增分点的组合。

取分点，计算得

这里，。

同理，计算如图7.6所示。

时只要在的基础上计算新增分点，的值再做组合，

图7.6 与

一般地，每次总对前一次的小区间分半，分点加密一倍，并可充分利用老分点上的函数值，每次只需计算新增分点的和。

对上等分， ，则有

记上的中点为则

（7.15）

其中。

或

其中。

类似地，可得积分节点为，的辛普森求积公式的关系式：

（7.16）

其中： 由误差公式：

由于可视

，于是

，分别为及个点上的均值，

上式表明的误差大约是误差的4倍。

或 （7.17）

由此得到启发，对任给的误差控制量

即可，而用

在计算机上也不难实现。

，要，只需

作为控制手段简单直接，序列

复化积分的算法描述

从数值积分的误差公式可以看到，截断误差随分点的增长而减少，控制计算的精度也就是确定分点数。在计算中不用数值积分的误差公式确定分点数的理论模式，而用

作为控制，通过增加分点自动满足精度的方法称为数值积分公式的自

动积分法。即在计算中构造序列

计算，由分点数自动控制积分值的误差，并取

，直到

。

或时停止

下面描述复化数值积分公式的自动控制误差算法，详细程序和算例请看本章7.6节。 1．输入：误差控制精度e = eps；初始分点值

。

2．计算分点的复化梯形积分

T1=T2+100 //迭代计算中T1和T2分别表示 3．while | T1－T2|＞e T1 = T2

和

H=Hn //计算新增节点的值 T2= (T1+H )/2

H = h/2，n =2n //将区间一分为二 end while

4．输出积分值T2。

在自动控制误差算法中初始分点值不宜过小，以防假收敛。

7.3.4 龙贝格（Romberg）积分

由前面得到的关系式（7.17），可以将 到 即

，

作为的修正值补充

（7.18）

其结果是将梯形求积公式组合成辛普森求积公式，截断误差由提高到。

这种手段称为外推算法。外推算法在不增加计算量的前题下提高了误差的精度，是计算方法中一种常用手法。

不妨对再做一次组合。由

得到

（7.19）

复化辛普森公式组成复化柯特斯公式，其截断误差是组合：

。同理对柯特斯公式进行

得到具有7次代数精度和截断误差是 的龙贝格公式：

还可以继续对

做上去。

为了便于在计算机上实现龙贝格算法，将统一用表示，列标 或步长

分别表示梯形、辛普森、柯特斯积分，行标表示分点数

j

。

龙贝格计算公式：

对每一个 龙贝格算法

从2做到，一直做到小于给定控制精度停止计算。

龙贝格算法按表7.2元素的行序进行运算， 在计算中每个元素只用到

上一行和本行的元素。对上面的算法进一步优化，对每k行可将计算定义在两行元素之间，令

；在每计算一行元素后，要将

表7.2 龙贝格算法计算元素顺序表

。

1．输入区间端点

；

，精度控制值，循环次数，定义函数，取

2．；

3．for {

=2 to

}

4．输出

。

7.4 重积分计算

在微积分中计算二重积分是用化为累次积分的方法进行的。计算二重数值积分也是计算累次数值积分的过程。为了简化问题，我们仅讨论矩形域上的二重积分。有很多非矩形域上的二重积分可作变换将其转换到矩形域上。

（7.20）

其中：次积分：

是常数，在上连续。像在微积分中一样，将二重积分化为累

（7.21）

或

二重积分的复化梯形公式

对区间和分别选取正整数，在轴和轴上分别有步

用复经梯形公式计算，计算中将当作常数，有

（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。

二重复化辛普森求积公式

对区间和分别等分和等分，在轴和轴上分别有步长

均为偶数

类似于二重复化梯形公式推导，得

记

误差：

和在积分区间内。

按例7.6的化分区间，的值如表7.5所示：

表7.5

7.5 高斯（Gauss）型积分公式介绍

对插值型积分公式

在牛顿-柯特斯积分公式中要求节点是等距的，其优点是计算积分系数的公式规则相同，缺点是制约了求积公式的代数精度。可以证明：当节点个数为偶数时，求积公式具有阶的代数精度；当节点个数为奇数时，求积公式具有阶的代数精度。

如果我们不预先指定求积节点的位置，和权系数都作为待定的常数，能否适当个待定常参数，需要

个

地确定它们，以提高积分公式的代数精度？回答是肯定的。方程来确定，取一个函数组：任一小于等于

，这一组函数构成了次多项式的基，

次的多项式，都可以用这组函数的线性组合来表示。如果某一积分公式，

次代数精度。

对这组函数都能精确积分，则此积分公式就有

例7.7 计算求积系数具有3阶代数精度。

和求积节点，使得至少

解：按照求积公式的代数精度定义，分别令，得方程组：

解方程组得：

求积公式：

，

按例7.7的方式，构造更高阶的代数精度的求积公式，生成求积系数和求积节点的方程组并无困难，而求解该方程组则无一定的章法可循。一般地，通过计算正交多项式的零点作为求积节点。

当取积分节点为正交多项式的零点时，则节点个数是的求积公式具有精度。并称积分节点为正交多项式的零点的数值积分公式为高斯型积分公式。 为了一般性，考虑积分

阶的代数

其中称为权函数。当时，即是普通的积分。

对于不同的权函数选定的节点也不相同。

如何构造高斯型积分公式呢？

对给定的及权函数，由施密特（Schmidt）正交化过程作出正交多项式

的个零点

，这个个零点就

；解出正交多项式

是积分节点；以这些节点构造插值多项式，计算积分系数

其中是拉格朗日插值基函数。

高斯型求积公式为

高斯型积分公式的优点是它的代数精度高，特别是对无穷区间或瑕积分更有效，但计算正交多项式的零点即积分节点有一定的工作量，好在数值学家们已算出一些特定的函数的积分节点和积分系数，在计算中我们可以查表直接得到这些数值。 本章并不构造各种高斯型积分公式，有关的详细内容请参考有关的教材。

下面给出上，取权函数的高斯型积分。

取[－1，1]上权函数的正交多项式为勒让德（Legendre）多项式：

高斯-勒让德相应的积分节点和积分系数表如下：

2 =－0.5773503，=－0.8611363， =0.5773503 1.0000000，1.0000000 =－0.3399810 0.3478548，0.6521452 0.6521452，0.3478548 4 =0.3399810， =0.8611363 0.2369269，0.2369269 －x1= x5=0.9061798 5 －x2= x4=0.538469 0.56888

0.4786287，0.4786287

=0.0 要计算一般区间上的积分，只需作变量代换则有

，其中

可用高斯积分求积，即

，这样 仍

例7.8 应用两点高斯-勒让德积分公式计算 解：

。

I = (－0.5773503)2cos (－0.5773503) + (0.5773503)2 cos (0.5773503) =0.558608

例7.9 应用两点高斯-勒让德积分公式计算 解：

。

令，得到积分

例7.10 证明不存在使求积公式

的代数精度超过次。

证明：只要能找到一个在求积系数和节点现取

及

次多项式，使求积公式两边不相等即可。用反证法，假定存

使求积公式对任何

次多项式

精确成立。

代入求积公式左端得

而公式右端不存在

，故右端与左端不相等，与假设矛盾，说明

次代数精度的求积公式，故高斯型求积公式是具有最高代数精度的求积公式。

7.6 程序示例

程序7.1 复化梯形公式法。

算法描述

的自动控制误差算

输入的值，定义；

；

T1=T2+100

；

while |T1－T2|＞ε T1 = T2；

H = h ！或用复化梯形公式计算T2

T2 = (T1 + H) /2；h= h/2；n = 2n； end while 输出T2。 程序源码

///////////////////////////////////////////////

// Purpose：梯形公式的自动控制误差算法 // ////////////////////////////////////////////// # include # include

# define f (x) (sin (x)) // f(x)由define定义 # define m 100 # define a 1.0 # define b 2.0

# define epsilon 0.00001 void main( ) { int n = m； int i；

double T n，H_n，T1，T2； double h = (b－a) / n；

T_n = (f (a) + f(b)) /2；

for (i =1；i＜n；i+

T_n + = f(a+h i)；

T_n= h； T2 = T_n； T1 = T2+100； do

{ T1 = T2；

for (i=0，H_n=0；i =1；i＜n；i+

H_n+=f (a+hi+h/2)；

H_n= h； T2 = (T1 + H_n) /2

h = h /2；n = n2；

} while (fabs (T1－T2)＞epsilon) printf (“T=% 1f \\ n”，T2)； }

// －－－－－－－－End of File－－－－－－－－－ 计算实例

对于

程序输入输出

，在区间上验证梯形公式的自动控制误差公式。

对于不同

。）

，修改程序# define f(x)项，本例，初始，区间

T = 0.9547 程序7.2 龙贝格积分算法

计算公式和算法描述第7.3.4节所述。 程序源码

////////////////////////////////////// // Purpose：龙贝格算法 // ////////////////////////////////////// # include # include

# define f (x) (sin (x)) # define N_H 20 # define MAXREPT 10 # define a 1.0 # define b 2.0 # define epsilon 0.00001

double compute T (double aa，double bb，long int n) // 复化梯形公式 {

int i；double sum，h = (bb－aa)/n；sum = 0；

for (i=1；i＜n；i+

sum + = f (aa + ih)；

sum + = (f (aa) + f (bb)) / 2；

return ( hsum)； }

void main ( ) { int i；

long int n = N_H，m = 0； double T[2] [MAXREPT +1]； T[1][0] = compute T (a，b，n)；

n = 2；

for (m = 1；m＜MAXREPT；m++) }

for ( i = 0；i＜m；i++ ) { T[0][i]=T[1][i];}

T[1][0]=computeT(a，b，n)；

n=2；

for ( i =1；i＜m + + ＝

T[1][i]=T[1][i－1]+(T[1][i－1]－T[0][i－1]) / (pow(2，2m)－1)； if ((T[0][m－1]－T[1][m])＜epslon)

{ printf (“T=% lf \\ n”，T[1][m])； // 输出 return； }

}

printf (“Return no solved… \\ n”)； }

// －－－－－－－－－－End of File－－－－－－－－－ 计算实例

利用龙贝格积分法计算 程序输入输出

。

（对于不同初始

，修改程序# define

。）

项，本例，区间，

对于本程序的给定，输出结果。

0.9547

友情提示：方案范本是经验性极强的领域，本范文无法思考和涵盖全面，供参考!最好找专业人士起草或审核后使用。