NIntegrate 积分规则

介绍

积分规则通常使用加权和对一个区域上的积分进行估值. 对于 NIntegrate 的用法而言,一个积分规则对象既进行积分估计,也进行误差估计,作为积分估计准确性的测量手段.

下面的一般说明与诸如 "GaussKronrodRule""MultidimensionalRule" 一类的加权总和类型积分规则相关. 适用于 "LevinRule" 等其它类型的规则另有单独的讨论.

积分规则在一组点上对被积函数进行采样,这些点称为采样点. 在文献中这些也称为横坐标. 与每一个采样点 相对应,有一个权重数 . 积分规则根据加权和 估计积分 . 积分规则是一个泛函,也就是说,它在区间 (或一个更一般的区域)上将函数映射到实数.

若一个规则适用于区域 ,则记为 ,其中 是被积函数.

对于以下讨论的规则,对积分进行估值的取样点选取于区间 、单位立方体 或者是位于中心的单位立方体 ,其中 是积分的维数. 所以如果 是这样一个区域,则使用 而不是 . 当这些规则应用于其它区域时,其横坐标和估计值需要相应地进行缩放.

如果 ,则认为积分规则 对于函数 是准确的.

一个积分规则 在一个函数 中的应用将被视为 的一个积分,例如, 积时,得到 ."

对于一个一维积分规则,如果被积多项式的次数都小于或等于 ,并且对于至少一个次数为 的多项式会失败,则称这个一维积分规则的次数为 .

对于一个多维积分规则,如果被积单项式的次数都小于或等于 ,并且对于至少一个次数为 的单项式会失败,则称这个多维积分规则的次数为 ,也就是说,对于所有形如 的单项式,积分规则是准确的,其中 为维数,,且 .

一个次数为 的空规则对于所有次数 的单项式的积分为零,且对于对于至少一个次数为 的单项式会失败. 每一个空规则可以被视作一个基本积分规则与一个次数较低的适当积分规则的差.

如果次数为 的规则 的样本点集包含一个次数为 的规则的 样本点集,且 ,则称 嵌入 ,并记作 .

如果次数为 的积分规则是具有共同推导和性能但不同次数的一个规则族的成员,则被记作 ,其中 可能被选出对规则族进行标识. (例如,4次梯形规则可以被称为 .)

如果在一个族中的每一个规则都是该族中另一个规则的嵌入,则称该族的规则是渐进的. (对于任意 ,存在 ,满足 ).

如果被积函数在区间的两个端点不进行计算,则称这个积分规则是开放型的. 如果在区间端点计算被积函数则称之为封闭型.

一个 NIntegrate 积分规则对象对于积分估计有一个积分规则,对于误差估计有一个或几个空规则. 积分规则和空规则的样本点是一致的. 根据上下文应该判断出积分规则规则指的是一个 NIntegrate 积分规则对象,还是通常数学意义上的一个积分规则.

积分规则规范

下述所有的积分规则,除了 "MonteCarloRule" 以外,都将用于NIntegrate自适应策略. 在NIntegrate 中,包括原始自适应的所有蒙特卡洛策略都使用 "MonteCarloRule".

改变一个积分策略的积分规则组件会导致不同的积分算法.

NIntegrate 中的自适应策略所用的积分规则(见 "全局自适应策略""局部自适应策略")是通过 Method 子选项指定的.

这个例子使用的是策略为 "GlobalAdaptive" 的积分规则:
这个例子使用的是与上例相同的积分规则,但积分策略不同("LocalAdaptive"):

如果 NIntegrate 给出的方法选项只有一个不为 "MonteCarloRule" 的积分规则规范,则该规则与 "GlobalAdaptive" 策略结合使用. 下面的两个输入是等价的.

这个积分只给定了积分规则. 默认情况下,使用 "GlobalAdaptive" 策略:
对于这个积分,积分策略和积分规则都指定了:

对于 "MonteCarloRule",自动使用自适应蒙特卡洛策略. 下面两个命令是等价的.

对于这个蒙特卡洛积分,只指定了 "MonteCarloRule"
对于这个蒙特卡洛积分, 蒙特卡洛策略和 "MonteCarloRule" 都指定了:

"TrapezoidalRule"

积分估计的梯形规则是最简单也是最古老的规则之一(可能曾被巴比伦人使用,但一定曾被古希腊数学家使用):

复合梯形规则是一个形如下式的 Riemann 和:

其中 .

如果选项 Method 的值为 "TrapezoidalRule",则复合梯形规则被用来估计由积分策略形成的每一个子区间.

一个 "TrapezoidalRule" 积分:
选项名称
默认值
"Points"5粗略梯形点数
"RombergQuadrature"True是否使用 Romberg 求积法
"SymbolicProcessing"Automatic进行符号预处理的秒数

"TrapezoidalRule" 选项.

梯形规则及其复合(多面板)扩展并不十分准确. (复合梯形规则对于线性函数是准确的,并且如果被积函数的二阶导数连续,则至少以 的速度收敛 [DavRab84].)多面板梯形规则的准确度可以通过使用 Romberg 求积法 提高.

由于 的横坐标是 的子集,差值 可用作积分估计 的误差估计,且不用额外计算被积函数就可求出.

选项 "Points"->k 可用来指定粗略点的个数. "TrapezoidalRule" 所用的总点数为 .

这里证实了样本点与 (1) 中相同:

"TrapezoidalRule" 也可以用于多维积分.

这是一个 "TrapezoidalRule" 的多维积分. 其准确结果是

评注:NIntegrate 既有梯形规则也有梯形策略;见教程积分策略中的" 'Trapezoidal' 策略". 所有 NIntegrate 的内部执行的规则都有后缀 -Rule. 因此 "TrapezoidalRule" 用于指定梯形积分规则,"Trapezoidal" 用于指定梯形策略.

Romberg 求积法

Romberg 求积法的思想是使用 的线性组合消去 的截断近似误差的阶数相同的项.

根据 EulerMaclaurin 公式 [DavRab84],有

其中

故可写作

如果从第二个方程中减去第一个方程4次,可以消去上述方程中的 项. 结果是

这个例子表明,使用 Romberg 求积法的梯形规则,其性能优于标准梯形规则. 并且前者的结果更逼近于准确结果
这是一个使用梯形规则但不用 Romberg 求积法的积分:

"TrapezoidalRule" 采样点和权重

下面计算的是给定精度的梯形样本点、权重和误差权重:
这里是 Romberg 求积法如何推导权重和误差权重:

"NewtonCotesRule"

NewtonCotes 积分公式是有相同间隔采样点的插值型公式.

NIntegrate 的 NewtonCotes 求积法可通过Method 选项值 "NewtonCotesRule" 来指定:
选项名称
默认值
"Points"3粗略 NewtonCotes 点数
"Type"ClosedNewtonCotes 规则的类型
"SymbolicProcessing"Automatic进行符号预处理的秒数

"NewtonCotesRule" 选项.

将积分区间 通过点

分成 个相同长度的子区间. 则插值型积分公式由

给出,其中

很大时,NewtonCotes 点系数也很大并且符号混合.

由于相互抵消,可能导致大量有效数位的丢失,因此使用高阶 NewtonCotes 规则时必须谨慎.

"NewtonCotesRule" 采样点和权重

下面计算给定精度的 NewtonCotes 采样点、权重和误差权重:

"GaussBerntsenEspelidRule"

高斯求积法使用最佳采样点(通过多项式插值法得到)形成在这些点上的被积函数值的加权和. 在这些采样点的子集上,可以使用一个较低阶的求积法规则. 两个规则的差可用来进行误差估计. Berntsen 和 Espelid 通过去除奇数个采样点的高斯规则的中点推导得到了误差估计规则.

NIntegrate 中的高斯求积法可以通过设定选项 Method 的值为 "GaussBerntsenEspelidRule" 来实现:
选项名称
默认值
"Points"Automatic高斯点的个数
"SymbolicProcessing"Automatic进行符号预处理的秒数

"GaussBerntsenEspelidRule" 选项.

对被积函数 个点的高斯规则 对于次数为 的多项式是准确的(也就是说,若 为次数为 的多项式,则 ).

由于不计算被积函数在区间两个端点处的值,因此高斯规则是开放型的. (Lobatto 规则, ClenshawCurtis 规则梯形规则 是封闭型的,因为它们计算被积函数在区间端点处的值.)

这里定义了差商泛函 [Ehrich2000]

对于采样点为 的高斯规则 ,Berntsen 和 Espelid 已推导得到了下面的误差估计函数(error estimate functional)(see [Ehrich2000])

([Ehrich2000] 中的原始公式针对的是区间 上的采样点. 上述公式针对的是区间 上的采样点.)

这个例子显示出 "GaussBerntsenEspelidRule" 的选项"Points" 取不同值时的 NIntegrate 所用的采样点个数:

"GaussBerntsenEspelidRule" 采样点和权重

下面计算了给定粗略点数和精确度时的高斯横坐标、权重和 BernsenEspelid 误差权重:

BerntsenEspelid 误差权重在下面被执行.

这里执行的是差商(divided differences):
这里计算 的横坐标和权重:
计算 BerntsenEspelid 误差权重:

"GaussKronrodRule"

高斯求积法使用最佳采样点(通过多项式插值得到)来形成在这些点处被积函数值的加权和. 高斯规则的 Kronrod 扩展是在高斯点之间增加新采样点而形成一个更高阶的规则,而该规则重新使用高斯规则被积函数的计算方法.

NIntegrate 中的高斯Kronrod 求积法可以通过设定Method 的选项值为"GaussKronrodRule" 来实现:
选项名称
默认值
"Points"Automatic与 Kronrod 点一起被扩展的高斯点的个数
"SymbolicProcessing"Automatic进行符号预处理的秒数

"GaussKronrodRule" 选项.

对被积函数 个点的高斯规则 对于次数为 的多项式是准确的,也就是说,若 为次数为 的多项式,则 ).

由于不计算被积函数在区间两个端点处的值,因此高斯Kronrod 规则是开放型的.

具有 个点的高斯规则 的 Kronrod 扩展 增加 个点. 扩展的规则对于次数为 为偶数)或 为奇数)的多项式是准确的. 与高斯规则相关的权重在其 Kronrod 扩展中会改变.

由于 的横坐标是 的一个子集,差值 可用于表示积分估计 的一个误差估计,可以在不用额外计算被积函数的情况下被计算.

这个例子显示 "GaussKronrodRule" 的选项 "Points" 取不同值时NIntegrate 所用的采样点的个数:

关于高斯规则的 Kronrod 扩展的执行介绍,请见[PiesBrand74].

"GaussKronrodRule" 采样点和权重

下面计算了给定粗略点数和精确度时的高斯Kronrod 横坐标、权重以及误差权重:

下面的计算显示出高斯-Kronrod 积分规则的次数(见).

计算高斯-Kronrod 积分规则的次数:
定义一个函数:

下面的命令执行了积分估计 和误差估计 的积分规则加权和,其中 是横坐标, 是权重, 是误差权重.

这是根据该规则计算的 的积分和误差估计:
积分估计正好是准确结果:

该误差估计非零,原因是嵌入的高斯规则对于次数 的多项式是准确的. 如果对这样次数的多项式进行积分,则误差估计为零.

定义一个函数:
按照该规则计算 得到的积分和误差估计:
这是使用 Integrate 得到的准确结果:

"LobattoKronrodRule"

Lobatto 积分规则是一种预先对横坐标赋值的高斯型规则. 它使用积分区间端点和区间内的最佳采样点来形成在这些点上被积函数值的加权和. Lobatto 规则的 Kronrod 扩展是在Lobatto 规则点之间增加新采样点而形成一个更高阶的规则,而该规则重新使用 Lobatto 规则被积函数的计算方法.

如果 Method 选项的值为 "LobattoKronrodRule",则 NIntegrate 使用 Lobatto 规则的 Kronrod 扩展:
选项名称
默认值
"Points"5与 Kronrod 点一起被扩展的高斯Lobatto 点的个数
"SymbolicProcessing"Automatic进行符号预处理的秒数

"LobattoKronrodRule" 选项.

对被积函数 个点的 Lobatto 规则 对于次数为 的多项式是准确的,也就是说,若 为次数为 的多项式,则 ).

具有 个点的 Lobatto 规则 的 Kronrod 扩展 增加了 个点,扩展了的规则对于次数为 为偶数)或 为奇数)的多项式是准确的. 与 Lobatto 规则相关的权重在其 Kronrod 扩展中会变化.

"GaussKronrodRule" 相同,选项 "GaussPoints" 设定高斯点的个数. 如果 "LobattoKronrodRule""Points"->n 激活,则规则点的总个数为 .

这个例子显示出 NIntegrate"LobattoKronrodRule" 的选项 "Points" 取不同值时的采样点的个数:

由于 Lobatto 规则 是封闭型的,需要计算区间端点的被积函数值. 如果在这些端点处有奇点,则将被 NIntegrate 忽略.

关于 Lobatto 规则的 Kronrod 扩展的执行介绍,请见[PiesBrand74].

"LobattoKronrodRule" 采样点和权重

下面计算了给定粗略点数和精确度时的LobattoKronrod 横坐标、权重以及误差权重:

下面的计算显示了 LobattoKronrod 积分规则的次数( 见).

计算 LobattoKronrod 积分规则的次数:
定义一个函数:

下面的命令执行了积分估计 和误差估计 的积分规则加权和,其中 是横坐标, 是权重, 是误差权重.

这是根据该规则计算的 的积分和误差估计:
上面的积分估计正好是准确结果:

上面的误差估计非零,原因是嵌入的 Lobatto 规则对于次数 的多项式是准确的. 如果对这样次数的多项式进行积分,则误差估计为零.

定义一个函数:
这是根据该规则计算的 的积分和误差估计:
使用 Integrate 得到的准确结果:

"ClenshawCurtisRule"

ClenshawCurtis 规则使用的采样点是通过被积函数的切比雪夫多项式近似得到的.

NIntegrate 的 ClenshawCurtis 求积法可通过设置Method 选项值为 "ClenshawCurtisRule" 得到:
选项名称
默认值
"Points"5粗略 ClenshawCurtis 点的个数
"SymbolicProcessing"Automatic进行符号预处理的秒数

"ClenshawCurtisRule" 选项.

理论上一个具有 个采样点的 ClenshawCurtis规则对于次数等于或小于 的多项式是准确的. 然而实际上ClenshawCurtis 规则能达到高斯规则 [Evans93][OHaraSmith68] 的准确度. ClenshawCurtis 公式的误差分析请见 [OHaraSmith68].

经典 ClenshawCurtis 规则的采样点是切比雪夫多项式的零点 .实用 ClenshawCurtis 规则的采样点是切比雪夫多项式的极值点. 经典 ClenshawCurtis 规则不是渐进的,但实用 ClenshawCurtis 规则是 [DavRab84][KrUeb98].

代表函数 个采样点的一个实用 ClenshawCurtis 规则.

渐进性质意味着 的采样点是 的采样点的子集. 因此差值 可用于表示积分估计 的一个误差估计,可以在不用额外计算被积函数的情况下被计算.

NIntegrate 选项 Method->{"ClenshawCurtisRule","Points"->k} 使用的是具有 个点的 的实用 ClenshawCurtis 规则:
这个例子显示出 NIntegrate"ClenshawCurtisRule" 的选项 "Points" 取不同值时的采样点的个数:

"ClenshawCurtisRule" 采样点和权重

下面是给定粗略点数和精确度时的采样点和ClenshawCurtis 规则的权重:
这是计算 的采样点的另一种方法:
定义一个函数:
这是根据该规则计算的 的积分和误差估计:
Integrate 得到的准确结果:

"MultipanelRule"

"MultipanelRule" 在两个或两个以上的相邻区间上将一维积分规则的应用合并至一个规则中. 在任何一个相邻区间上的原始规则的应用被称为一个面板.

这是一个使用 "MultipanelRule" 的积分的例子:
选项名称
默认值
Method"NewtonCotesRule"为一个单面板提供横坐标、权重和误差权重的积分规则设定
"Panels"5面板个数
"SymbolicProcessing"Automatic进行符号处理需要的秒数

"MultipanelRule" 选项.

将单位区间 用点 分为 个子区间.

如果有规则

它可在区间 上变形为一个规则,

,且 , . 则基于 面板积分规则可被显式写为

是封闭的,即 存在 作为采样点,则 ,且 的采样点个数可减至 . (这通过执行 "MultipanelRule" 完成.)

更多关于多面板规则的理论(也被称为复合规则),请见 [KrUeb98] 和 [DavRab84].

"MultipanelRule" 采样点和权重

"MultipanelRule" 的采样点和权重可以通过这个命令得到:
这是高斯Kronrod规则的横坐标和权重:
多面板规则横坐标可以通过 Rescale 得到:
这里表明如何由原始权重推导得到多面板规则权重:

"CartesianRule"

一个 维 Cartesian 规则的采样点是 个一维规则的采样点的 Cartesian 乘积. 与 Cartesian 规则的采样点相关的权重是对应于其坐标的一维规则权重的乘积.

NIntegrate 的 Cartesian 乘积积分可通过设定 Method 选项值为 "CartesianRule" 实现:
选项名称
默认值
Method"GaussKronrodRule"形成 Cartesian 乘积规则的一个规则或规则列表
"SymbolicProcessing"Automatic进行符号预处理需要的秒数

"CartesianRule" 选项.

举例来说,假设有下面的公式:

这些公式分别对次数为 的多项式是准确的. 则不难看出 个点的公式

对于次数为 中的多项式是准确的. 注意与横坐标 相关的权重是 .

个一维规则(其中 个规则有 个采样点 和权重 ),其一般 Cartesian 乘积公式是

显然 (2) 可以写作

其中 ,并且对于每一个整数 .

这是 Cartesian 乘积规则积分的图示. 在 轴上使用 "TrapezoidalRule";在 轴上使用 "GaussKronrodRule"

Cartesian 规则适用于维数较低的情形(),原因在于对于较高的维数,它们受制于"组合爆炸". 例如 个相同的一维规则,其中每个规则有 个采样点,则对应的五维 Cartesian 乘积将有 个采样点.

如果积分为多维,且 Method 选项被给出一个一维规则或一个一维规则列表,则 NIntegrate 使用 Cartesian 乘积规则.

这个例子使用 GaussKronrodRule 设定了 Cartesian 乘积规则积分:
这里使用一个一维积分规则列表进行 Cartesian 乘积规则积分的设定:
这又是一个使用一维积分规则列表进行 Cartesian 乘积规则积分的设定的例子:

更多关于 Cartesian 规则的讨论可以在 [Stroud71] 中找到.

"CartesianRule" 采样点和权重

"CartesianRule" 规则的采样点和权重可以通过命令 NIntegrate`CartesianRuleData 得到.

NIntegrate`CartesianRuleData 使每个规则的横坐标和权重保持分开. 否则,正如(2)中所看到的,对于较高的维数,结果可能过大.

NIntegrate`CartesianRuleData 的结果可以通过这个函数放到(2)的形式中:

"MultidimensionalRule"

组成立方体 的一个完全对称积分规则由具有下述性质的的点的集合组成:(1) 一个集合中的所有点可以通过置换及/或改变集合中任意固定点坐标的符号生成;(2) 一个集合中的所有点都具有与之相关的相同权重.

NIntegrate 的完全对称多维积分(完全对称求体积法)可以通过设置 Method 选项值为 "MultidimensionalRule" 实现:

满足上述性质的一个完全对称积分规则的点集被称为一个轨道. 轨道的一个点 ,若其坐标满足不等式 ,则被称为一个生成元. (见 [KrUeb98][GenzMalik83].)

选项名称
默认值
"Generators"5完全对称规则的生成元个数
"SymbolicProcessing"Automatic进行符号预处理需要的秒数

"MultidimensionalRule" 选项.

如果积分规则有 个分别记为 的轨道,且其中的第 个轨道 的权重为 ,则积分估计可通过下述公式进行计算:

一个次数为 的空规则对于所有次数 的单项式的积分为零,且对于对于至少一个次数为 的单项式会失败. 每一个空规则可以被视作一个基本积分规则与一个次数较低的适当积分规则的差.

NIntegrate"MultidimensionalRule" 对象基本上是三个积分对象的接口,将一个积分规则和一个或几个空规则合并. 下面的表格对生成元的个数和阶数作一总结. 具有6个和9个生成元的规则对象使用三个空规则,其中每一个是两个空规则的线性组合. 使用空规则的线性组合是为了防止相位误差. 关于如何使用空规则的更多细节请见[BerntEspGenz91].

NIntegrate 的完全对称多维积分中生成元的个数和阶数:

生成元的个数积分规则阶数每一个空规则的阶在何处被介绍
575[GenzMalik80]
675,3,1[GenzMalik83][BerntEspGenz91]
997,5,3[GenzMalik83][BerntEspGenz91]

NIntegrate 使用的采样点的个数,其中完全对称多维积分用于形如 , 的积分:

"MultidimensionalRule" 采样点和权重

这一小节给出一个用完全对称多维积分规则计算积分估计的例子.

这是生成元个数的参数:
这个函数取一个生成元点并创建它的轨道:
给定生成元个数的生成元和权重:
计算每一个生成元的轨道:
定义一个函数:
应用多维规则:
这是准确结果:
这里作出轨道中点的图形基元:
这里显示出轨道如何不同:
将所有规则的点集合在一起:

"LevinRule"

Levin 型规则通过用多项式和被积函数振荡部分的乘积对不定积分进行近似来估计振荡函数的积分.

这是利用 "LevinRule" 计算的一个高振荡积分:

默认时,"LevinRule" 自动检测被积函数的振荡部分,并应用下述配置方法来形成一个积分估计. 可以使用选项来指定振荡部分,并控制数值算法的细节行为,包括是否自适应地切换到一个备选的非振荡积分规则.

选项名称
默认值
"Points"Automatic粗略的配置点数
"EndpointLimitTerms"Automatic在非数值端点的多项式外推阶
"TimeConstraint"10求解配置方程的最大时间
"MethodSwitching"True是否自动切换到非振荡规则
"OscillationThreshold"10应用 "LevinRule" 所需超过的近似振荡数
MethodAutomatic备选的非振荡规则
"LevinFunctions"Automatic包括在内核中的振荡函数
"MaxOrder"50内核的最大微分阶数
ExcludedForms{}被内核排除在外的形式
"AdditiveTerm"Automatic被积函数内附加的非振荡部分
"Amplitude"Automatic被积函数中内核的系数
"DifferentialMatrix"Automatic内核满足的线性常微分方程的矩阵形式
"Kernel"Automatic显式振荡内核

"LevinRule" 选项.

Levin 型配置方法

范围

一维数值积分的基本 Levin 型方法 [Levin96] 适用于如下形式的振荡积分:

其中,振荡部分 必须满足阶数为 的线性常微分方程(线性 ODE),以矩阵形式写作:

其中 中的一个. 函数 的向量是振荡内核,矩阵 微分矩阵.

举例来说,如果 ,振荡内核可以是

微分矩阵为

Levin 型方法也可以处理更一般的情形:

通常,上式比较简单的情况是 . 函数 的向量是振幅.

绝大多数振荡型特殊数学函数,如 ExpBesselJAiryAi 等,满足线性常微分方程,因此可以利用 Levin 型方法进行积分. 此外,此类函数的任何乘积、求和或整数幂,或者是它们与任意非快速振荡函数的复合,仍满足线性常微分方程. 另请参见 DifferentialRootReduce.

"LevinRule" 方法可以应用于一种更一般的积分类别:

其中,附加项 通过一种传统的求积法,例如 "GaussKronrodRule",单独处理.

当振幅 、附加项 以及微分矩阵 自身非快速振荡时,Levin 型方法非常有效 [Levin97].

配置方法

在 Levin 型配置方法中,被积函数 的一个不定积分以 的形式求出,其中 是一些未知函数:

对等式两边关于 求导,并利用 ,我们发现所有的 可以消去,并且函数 满足非振荡微分方程:

利用一种配置方法,这些方程中的 可以求解,在该配置方法中,我们将每个 近似表示成一组 个预定义的简单基函数 的线性组合:

"LevinRule" 中,选取切比雪夫多项式 ChebyshevT[k,x] 作为基函数 . 将其代入关于 的微分方程,得到 个配置系数 个线性方程:

由于 均为已知函数,我们可以在积分范围内部的 个配置节点 处计算这些线性方程. 这将生成关于 个系数 个线性方程. 一旦得到这些方程的数值解,我们将得到 的近似形式,并能够利用微积分基本定理来直接计算积分,如:

NIntegrate 使用与 "GaussKronrodRule" 对应的横坐标作为配置节点 .

"LevinRule" 利用 个配置节点计算积分估计,其中 "Points" 方法选项的值. 同时,一个较低阶数的积分估计也利用每隔一点的配置节点 () 计算得到,这两个积分估计之差作为对算法误差的估计.

在每个积分子域上利用23个配置节点计算积分估计,利用11个配置节点形成用作误差估计的一个较低阶数的积分估计:

指定振荡内核

默认时,"LevinRule" 自动选择构成被积函数尽可能多的部分作为振荡内核. 该行为可以通过选项进行细节修正.

如果被积函数包含的振荡函数在积分区域上不是高度振荡的,将其排除在振荡内核以外可能会更有效. 例如,Sin[x] 在区域 {x,0,5} 上不是高度振荡的.

"Kernel" 方法选项可用于指定被积函数的振荡部分.

使用内核为 Sin[1000x]"LevinRule". 因子 Sin[x] 没有包括在振荡内核中,而是包括在振幅中:

"ExcludedForms" 方法选项可用于指定不应包括在自动检测的振荡内核中的被积函数部分.

这是指定上述同一内核的另一种方法,利用的是 "ExcludedForms" 选项:

"LevinFunctions" 方法选项可用于指定哪种类型的函数可以包括在自动检测到的振荡内核中. 设置"LevinFunctions"->{h1,h2,} 仅检测具体函数 hi.

仅将 "LevinRule" 用于 ExpSin. BesselJ 被忽略:

设置 "LevinFunctions"->{grp1,grp2,} 检测已命名组 grpi 中的函数.

使用 "LevinRule",仅检测与三角函数和指数函数相关的函数. BesselJ 不被包括在振荡内核中:

"LevinFunctions" 方法选项识别的已命名组.

缺省时,"LevinFunctions" 的设置除了包括 "Other""InverseTrig""Elliptic" 外,还包括所有已命名组. 该设置对应于定义域某些部分振荡的初等及特殊函数.

NIntegrate`LevinIntegrand 对象

利用内部函数 NIntegrate`LevinIntegrandReduceNIntegrate 探测的振荡内核、微分矩阵、振幅与附加项进行详细检查是可能的. NIntegrate`LevinIntegrandReduce 返回一个 NIntegrate`LevinIntegrand 对象,该对象表示一个可以应用 Levin 型配置方法的完全处理的被积函数.

这是一个处理后的 NIntegrate`LevinIntegrand,被积函数 +x x 分解成了一个振荡内核和其它部分:

NIntegrate`LevinIntegrand 对象 li 的属性可以利用 li["prop"] 形式访问.

这是上述被积函数的振荡内核
这是同一被积函数的微分矩阵 . First 被启用的原因是在多维情形中,在积分的每一维度上,都有一个不同的微分矩阵:
这里是如何验证该振荡内核的 Levin 微分方程
属性的完全列表可以通过 li["Properties"] 访问:

属性 "Rules" 给出用于 Levin 型配置算法的 NIntegrate`LevinIntegrand 的主要属性.

除了上述被积函数的振荡内核和微分矩阵外,显示附加项 和振幅

利用该属性我们可以迅速查看不同 "LevinRule" 方法选项关于振荡内核检测的效果,例如 "AdditiveTerm" 选项.

这是同一被积函数的两种不同分解. 缺省时,Sin[x]Cos[x] 项均被包括在内:
在设置为 "AdditiveTerm"->Sin[x] 时,Sin[x] 项不被包括在检测到的内核中:

我们可以看到如何使用 "LevinFunctions" 选项将函数组包括在被积函数中或从被积函数中排除.

这里是如何查看 "LevinFunctions" 方法选项的不同设置下的效果. 缺省时 Log[x] 从检测到的内核中被排除,因为它不是振荡函数:
在设置 "LevinFunctions"->All 下,Log[x] 被包括在检测到的内核中:

完全指定被积函数分解的全部四部分(附加项、振幅、内核、微分矩阵)是可能的. 这可以用于对相似被积函数进行 NIntegrate 的有效多个调用.

这是一个关于如何使用 "LevinRule" 的例子,被积函数分解完全是手动指定的:

注意在这种情况下,NIntegrate 的第一个参数式完全被忽略了.

下面的积分与上面完全等价,NIntegrate 的第一个参数被忽略:

"MonteCarloRule"

蒙特卡洛规则通过形成一个在随机(半随机)采样点处被积函数值的统一加权和进行积分估计.

这个例子使用具有1000个采样点的 "MonteCarloRule"
选项名称
默认值
"Points"100采样点个数
"PointGenerator"Random采样点坐标生成器
"AxisSelector"Automatic全局自适应蒙特卡洛积分被使用时分裂轴的选择算法
"SymbolicProcessing"Automatic进行符号预处理需要的秒数

"MonteCarloRule" 选项.

在 [KrUeb98] 的蒙特卡洛方法中, 维积分 被解释成下述期望(平均)值:

其中 是函数 相对于在 上的均匀分布,被解释为一个随机变量时的平均值,即概率密度为 的分布. 代表区域 上的特征函数; 代表 的体积.

期望值 的原始蒙特卡洛估计通过选取密度为 个独立随机向量 得到(也就是说,向量在 上均匀分布),并且得到估计

评注:函数 是一个有效的概率密度函数,因为它在整个 上非负,并且 .

根据强大数定律(strong law of large numbers),收敛

的概率为 . 强大数定律不提供关于误差 的信息,因此使用一种概率估计.

定义

公式 (3) 是 的无偏估计(即,对于 的各种集合的期望值 )并且其方差是

其中 代表 的方差,这样 的标准误差是 .

实际操作时, 是未知的,因此它通过下式进行估计:

于是 的标准误差是

蒙特卡洛估计的结果可以写作 .

从方程 (4) 中可以看出,原始蒙特卡洛估计的收敛率不依赖于积分的维数 ,并且如果 个采样点被使用,则收敛率为 .

NIntegrate 的积分规则 "MonteCarloRule" 计算 的估计值.

估计可以逐步得到改善. 也就是说,如果已知 的估计和样本函数值 的一个新的附加集合,则根据 (5) 和 (6),我们有

计算 的估计时,没有必要知道计算 的估计所用的随机点.

"AxisSelector"

当用于多维全局自适应积分时,"MonteCarloRule" 采用两种方式选取一个积分子区域的分裂轴:(i) 随机选择,或 (ii) 如果子区域沿该轴被划分,则通过使每一半子区域的积分估计的方差和最小来选取. 分裂轴选取在积分估计之后进行.

通过下述方式进行轴的随机选择. "MonteCarloRule" 保留了供选择的轴的一个集合 . 初始情况下, 包括所有的轴. 的一个元素被随机选出. 选定的轴从 中被排除. 在下一个积分估计之后,由 中选择一个轴,并被排除在外,如此进行下去. 如果 为空集,则被所有轴填充.

轴的最小方差选择通过下述方式完成. 在对该区域进行积分的过程中,采样点的一个子集及其被积函数值被存储. 然后对于每个轴,使用存储的采样点和相应的积分值来生成沿着这条轴线分开的两个子区域的方差. 使方差和最小的轴被选作分裂轴,因为这将意味着如果该区域被这个轴分开,则新的积分误差估计将最小. 如果对于某些轴,所有存储点集中于一个半区,则选用该轴进行区域的分割.

选项值
Random随机分裂轴选择
MinVariance|{MinVariance, "SubsampleFraction"->frac}分裂轴选择使新的区域的方差和最小的分裂轴选择

"AxisSelector" 选项.

选项名称
默认值
"SubsampleFraction"1/10用来确定分裂轴的采样点比例

MinVariance 选项.

这是一个使用 "MonteCarloRule" 的选项 "AxisSelector" 的例子:

在下面的的例子中,将两种选择轴的算法进行比较. 一般来说,使方差最小的选择使用的采样点较少,然而这种选法会使"MonteCarloRule" 的应用减慢. 因此,如果两种轴选择方法的使用导致相同的采样点数目,则使用随机轴选择方法较快. 此外,使用较大比例采样点,以选择使方差最小的分裂轴会导致积分速度的减慢.

考虑下面的函数:
这是对上述函数进行分裂轴随机选择的自适应蒙特卡洛积分采样点:
这是使方差最小化的分裂轴选法的采样点. 与前面的蒙特卡洛积分相比较,这里的采样点更集中于圆 的周围,并且其数目几乎是先前的1/2:
这是使用随机轴选择方法的自适应蒙特卡洛积分积分:
这是对前面的积分使用方差最小化轴选择方法的自适应蒙特卡洛积分,其速度比使用随机轴选择方法慢:
使用存储点的较大部分进行方差最小化的选择使积分速度减慢:

各种规则的比较

除了 "MonteCarloRule" 外,所有的规则都用于 NIntegrate 中的自适应策略. 一个积分策略的积分规则组件的类型及采样点的数目的更改是导致不同的积分算法. 一般来说,这些不同的积分算法在不同的积分中表现不同. 自然地会出现下述问题.

1.  对于某种类型的任意一个或几个积分,有没有一种类型的规则总是优于其它类型的规则?

2.  给定一种积分策略,哪一种规则的表现较好?是针对何种积分?

3.  已知一个积分,一种积分策略,以及一种积分规则,在这种规则下使用多少点能够在达到积分估计的目标且满足精确度的要求的前提下使采样点的总数最小?

对于一个给定的积分和积分策略,如果一种积分规则满足精确度要求,且采样点的数目最少,则称这种规则为最佳积分规则. 下面几个因素决定了最佳积分规则.

4.  一般来说,规则的次数越高,对于平滑的被积函数及较高的精确度目标,积分速度越快. 另一方面,规则次数有可能相对于被积函数而言过高而导致自适应策略在被积函数的不连续点处工作时的采样点过多.

5.  一个规则的误差估计函数显著影响积分策略的总工作量. 使用较少点的规则可能引起 (1) 由于对积分过低估计的错误结果,或者(2) 由于对积分过高估计而使用过多的采样点. (见 "病态行为范例".) 此外,误差估计函数可能采用一个或几个嵌入的空规则进行计算. 一般来说,空规则的数量越多,误差估计的效果越好因为相位误差越少. 空规则的数量和及其在计算误差估计的加权和中被指定的权重决定了病态积分和难以用该规则进行计算的积分的集合. (NIntegrate 中的一些多维规则使用嵌入的空规则进行误差估计的计算. NIntegrate 中的所有一维积分规则都只使用一种空规则.)

6.  局部自适应策略对于采样点趋于均匀分布的的封闭规则(例如 "ClenshawCurtisRule") 比开放规则(例如GaussKronrodRule)和采样点非均匀分布的的封闭规则(例如 "LobattoKronrodRule")更为有效.

7.  被策略重复使用的点的百分比可在很大程度上决定哪一种规则最好. 对于一维积分,"LocalAdaptive" 重复使用封闭规则的所有点. "GlobalAdaptive" 则放弃误差估计有待改善的区域的几乎所有点.

一个规则中点的数目

这一小节用实例说明规则的次数越高,对于平滑的被积函数及较高的精确度目标,积分速度越快. 同时也用例子说明规则次数相对于被积函数过高而导致自适应策略在被积函数的不连续点处工作时使用的采样点过多. 这些例子均使用高斯规则及 BerntsenEspelid误差估计.

这是高斯规则在区间 上的误差.

(见 [DavRab84].)

这个函数计算规则对于积分 的误差,为了便于比较,使用由 Integrate 计算出来的准确值:
定义一个函数列表:
这是在区间 上函数的图形:
计算对于一定范围的点,对 , , 使用 "GaussBerntsenEspelidRule" 的误差:
这是对于每个函数,误差的对数值如何减小的图形. 可以看出,对于不连续函数和具有不联系导数的函数的积分估计随点数的增加改善很慢:

采样点的最小数目

这是一个求积分中使用的采样点数目的函数:
这里找到在一定范围精确度目标和一定范围积分规则粗略点中使用的采样点的数目:
这里找到对于每一个精确度采样点的最小总数. 通过这种方式,也就找到了粗略积分规则所用的点的个数:
这是使用的采样点的总数最少时的精确度目标和积分规则点数的图形:

规则比较

这个函数计算一个规则对于积分 的误差,使用由 Integrate 计算得到的准确值以作比较:
这里定义一个函数列表:
这是在区间 上的函数图形:
计算对于一定范围的点,对 , , 中的每个积分使用"GaussKronrodRule", "LobattoKronrodRule", "TrapezoidalRule""ClenshawCurtisRule" 的误差:
这是对于每种规则和每个函数,误差的对数值如何减小的图形.

病态行为范例

针对误差估计的错误行为

这一小节将讨论这样一种积分:NIntegrate 在默认设置下由于不能检测被积函数的一部分而对积分估计不足. 如果精确度目标提高则这一部分将被检测到.

错误估计

考虑下述函数:
这是在区间 上的准确积分:
NIntegrate 给出估计:
与准确值相比,这一结果太不准确:
这是函数的图形,也是错的:

较好的结果

使用 NIntegrate 选项 PrecisionGoal 和提高递归深度(recursion depth)可以达到较好的结果:
这个表中找到了没有好结果被计算出来的精确度目标:
如果增加绘图点数,函数的图形看上去不同:
这是放大后尖峰的图形,它在 Plot 默认选项中被丢掉了:
如果对函数的这部分进行积分,结果与默认设置下 NIntegrate 丢失的量相符合:

为什么估计受到误导

这是 NIntegrate 默认使用的 GaussKronrod 规则的横坐标和权重:
定义一个函数应用该规则:
得到自适应策略对被积函数取样的点:
这是采样点的图形. 垂直轴代表点被用来进行被积函数计算的阶数:

可以看出在前面的图形中 NIntegrate 在靠近 的第二个尖峰处进行大量计算. NIntegrate 不在靠近 的不积分的峰处进行同样量的计算.

这是在上一采样点的集合中的高斯-Kronrod 和高斯横坐标,它们在 上的区域:
这里被积函数应用于横坐标:
这是被积函数在横坐标上的多项式近似:
这些图形表明两个多项式近似在 上几乎重合:
如果多项式在放置 的区域进行积分,被 NIntegrate 用作误差估计的两者的差别的确很小.

由于这个差值就是区域 所指定的误差估计,在默认的精确度目标下,NIntegrate 不会用它进行进一步的积分完善.

相位误差

这一小节将讨论导致积分规则对积分估计的实际误差严重过高或过低估计的原因. 在 [LynKag76] 中给出了类似的讨论.

定义一个函数:
考虑 f[x,0.415,1.25] 在区域 上积分的数值和符号计算:
两者相差甚远. 要求达到的精确度目标是2,但相对误差比 要高得多:

(注意对于较高的精确度目标,NIntegrate 给出正确结果.)

下面对这一现象作出解释.

令积分规则 嵌入到规则 中. 偶尔,积分估计 (其中 )的误差估计 与准确误差 相比可能过小.

为了说明这一点,考虑具有11个采样点的高斯Kronrod 规则 ,该规则 具有一5个采样点的嵌入高斯规则. (这是上面两个积分所用的规则.)
定义一个函数,应用该规则:
这是先前定义的 f[x,λ,μ] 的积分
可以将 的误差估计和在 上取不同 值的实际误差绘制到一个图形中. 也就是说,绘制

在上面的图形中,蓝色图线代表估计的误差 . 实际误差 的图线为红色.

可以看到,参数 的值 0.415 非常接近于 的一个局部最小值.

一个一维求积法规则可以看作一个根据规则的横坐标及其对应的被积函数值拟和的多项式的积分结果. 我们可以进一步尝试看到 f[x,λ,μ] 的积分的拟合多项式.

下面的图形中,函数 f[x,λ,μ] 为红色,高斯多项式为蓝色,高斯Kronrod 多项式为紫色,高斯采样点为黑色,高斯Kronrod 采样点为红色.

可以看到由于 f[x,0.415,1.25] 的峰值大致落到两个横坐标的中间,这一近似为过低估计.

相反地,由于 f[x,0.53,1.25] 的峰值大致落到一个横坐标上,该近似为过高估计.

技术术语索引

横坐标

一维积分规则的次数

多维积分规则的次数

准确规则

嵌入规则

空规则

乘积规则

渐进性规则

采样点