本节书摘来自华章计算机《算法基础》一书中的第2章,第2.5节,作者:(美)罗德·斯蒂芬斯(Rod Stephens)著,更多章节内容可以访问云栖社区“华章计算机”公众号查看
2.5 进行数值积分
数值积分,有时也称为求积或数值求积,即用数值计算技术来逼近由函数定义的曲线下面积的方法。通常该函数是一个变量y的函数y=F(x)。这样的结果是一个二维区域,但是,一些应用可能需要计算一个函数z=F(x, y)所限定的表面下的三维区域,甚至还需要计算更高维函数定义的区域。
如果函数是很容易理解的,则可以用微积分来找到确切区域。但是在一些情况下也许无法找到函数的原函数。例如,也许该函数的公式很复杂,或者仅仅有一些物理过程中产生的数据,无从知道函数的方程。在这种情况下,不能使用微积分,但可以使用数值积分。
有几种方法来执行数值积分。最直接的方法是牛顿-柯茨公式,它使用了一系列多项式来逼近函数。在这之中最基本的两种牛顿-柯茨公式是矩形法则和梯形法则。
2.5.1 矩形规则
该矩形规则使用一系列宽度均匀的矩形来逼近曲线下区域的面积。图2-2显示了矩形规则示例程序(可以在这本书的网站上下载),该程序显示了使用矩形规则的情形。 该程序还使用微积分来计算曲线下的确切区域,这样一来可以看到用矩形规则得到的结果与正确的结果相差多少。
下面的伪代码显示了算法应用矩形规则的情形:
该算法简单地将区域划分成一些矩形,这些矩形具有一定的宽度,并且它们的高度等于函数中矩形左边缘的函数值。然后,遍历这些矩形,把它们的面积相加。
2.5.2梯形规则
图2-2显示出了矩形并不一定准确地适合曲线,计算出的总区域会产生一个误差。可以通过使用更多的、窄的矩形来减少误差。在本示例中,矩形的数目从10个增加到了20个,误差大约从6.5%降低到了3.1%。
另一种策略是使用梯形来代替使用矩形逼近曲线。图2-3显示了梯形规则示例程序(可以在这本书的网站上下载)。
下面的伪代码显示了算法应用梯形规则的情形:
此算法和矩形规则算法之间的唯一差别是在于添加了每个切片的面积。本算法采用的公式为梯形的面积公式:面积=宽×平行侧面的长度的平均值。
可以把矩形规则看作是用一个阶梯函数来逼近曲线,这个函数的函数值在每个矩形的边缘处从一个值变化到另一个值。梯形规则用线段来逼近曲线。
牛顿-柯特斯公式的另一个例子是辛普森法则,它使用二次多项式来逼近曲线。还有其他方法使用了更高次的多项式来更好地逼近曲线。
2.5.3 自适应求积
到目前为止所描述的数值积分方法的一个变化是自适应求积,在自适应求积时,程序检测到它的近似方法可能会在某些地方产生较大的误差,并在那些地方细化其方法。
例如,再次看一看图2-3,区域中的曲线接近于直线,梯形非常密切地逼近曲线,然而,在曲线的弯曲急剧的地区,梯形也不适合。
采用自适应求积的程序可以找出有哪些地方梯形不能很好地逼近曲线,并在这些地方使用了更多的梯形。
如图2-4所示,该自适应中点求积示例程序采用了自适应求积梯形规则。当计算切片的面积时,这个程序首先使用单一的梯形来逼近该区域。然后,程序将其分成两部分,并且使用两个较小的梯形来计算他们的面积。如果较大梯形的面积、较小梯形的面积和之间的差超过一定比例时,程序再将其分成两部分,并以相同的方式计算这些部分的面积。
下面的伪代码描述了这个算法:
如果运行自适应中点求积程序,开始时只有两个最初的切片,该程序将其划分为如图2-4所示的24个切片,并且估计出曲线下的面积,这样有0.035%的误差。如果使用具有均匀宽度的24个切片的梯形规则程序来计算,则该程序具有0.072%的误差,大约两倍于自适应程序产生的误差。虽然这两个程序使用的片数相同,但自适应程序能更有效地定位它们。
自适应梯形求积示例程序使用不同的方法来决定何时把切片分成子切片。它可以在该切片的起始的x值处计算该函数的二阶导数,并把此区间分成二阶导数值加一的份数。例如,如果二阶导数是2,程序将其分成三段。 (片数的计算公式选得有些武断,可以用一个不同的公式得到更好的结果。)
注 如果对积分有点生疏,函数的导数可以给出它在任何给定的点的斜率。它的二阶导数给出斜率变化速率,或曲线弯曲得有多快。较大的二阶导数值指出该曲线是弯曲得相对比较紧,所以自适应梯形求积程序采用更多的切片。
当然,如果不能计算曲线的二阶导数,这一技术将无法正常工作。而所采用的自适应中点求积程序的技术似乎在任何情况下工作得相当好,并且值得信赖。
自适应技术在许多算法中都有用到,因为它们能产生更好的结果,而不会在不需要的区域浪费精力。如图2-5所示,该自适应网格求积程序采用自适应技术来估算阴影区域的面积。此区域是纵向和横向的椭圆的结合再减去椭圆内的三个圆的面积 。
这个程序将整个图像变成一个方格,并确定方格内点的网格。图2-5的程序使用了四排和分列网格,对于网格中的每个点,该程序确定该点位于内部或阴影区域之外。
如果方格中点没有位于阴影区域内,该程序假定方格不在该区域内,并忽略它。
如果方格中的每一个点都位于阴影区域内,该程序认为方格完全位于该区域内,并将该方格的面积加入该区域的待测面积。
如果方格中有一部分点位于阴影区域内部,而另外一部分位于区域之外,该程序将方格细分成更小的方格并用同样的方法来计算这些更小的方格的面积。
在图2-5中,自适应网格求积程序已经制定了它认为合适的方格,所以可以看到它们。这样就可以看到,该程序考虑得更多的是将方格位于阴影区域的附近而不是位于区域内或区域外。总体而言,本示例中程序将方格分成1921个,并且主要集中在被整合区域的边缘。
2.5.4 蒙特卡罗积分
蒙特卡罗积分是一种数值积分的方法。在这种方法中,程序在一个均匀的区域内产生一系列的伪随机点,并确定是否每个点都位于目标区域内。当上述步骤完成后,程序使用区域内的点占所有点的百分比来估计该区域的总面积。
例如,假设算法产生的点是一个20×20平方单位的集合,所以它具有400平方单位的面积,并且伪随机点中有37%是在区域内。该区域有大约0.37×400=148平方单位的面积。
如图2-6所示的蒙特卡罗求积示例程序使用这一技术来估计在自适应网格求积程序中所使用的相同形状的区域面积:
相比起那些更有条理的方法(例如梯形规则和自适应求积),蒙特卡罗积分一般是比较容易出错,但有时它更容易计算,因为不要求很了解待测量的模型的性质。只需向那个模型抛出点,看看有多少打中它。
注 本章介绍如何使用伪随机数计算面积,但更普遍地是,可以使用类似的技术来解决很多问题。在蒙特卡罗模拟中,通过挑选伪随机数,看看有多少比例满足一定的标准来估计符合标准值的总数。