NIntegrate 积分策略

介绍

积分策略是一个算法,试图计算满足用户指定的精度和准确度目标的积分估计.

通常,积分策略规定如何管理和建立一个初始积分区域的不相交的子区域的新元素. 每个子区域都可以有自己的被积函数和与其相关的积分规则. 积分估计是所有子区域的积分估计的总和. 积分策略运用积分规则来计算子区域的积分估计. 一个积分规则在一个点集上对被积函数上进行采样,该点集称作采样点(或横坐标).

为了改善积分估计,应该对被积函数的更多的点取样,这有两种主要方法:(1)自适应策略,试图找出存在问题的积分区域并在这个区域集中计算(如采样点);(2)非自适应策略,提高抽样在整个区域上增加采样点,重用前积分估计的被积函数的计算结果,从而计算一个积分次数更高的积分规则估计高的规则估计.

两种方法均能采用符号预处理和变量变换或者是序列求和加速以实现更快的收敛速度.

下面的积分中,NIntegrate 中的符号分段预处理器将被积函数识别为一个分段函数,当 时,被积函数 ,当 时,被积函数为
这是积分中所用采样点的图形. 被积函数按图形中 坐标的次序在 轴上取样. 可以看出采样点集中在奇点 附近. 由于应用了一种奇点处理器,图形上部与下部采样点的格局不同:

"自适应策略" 一节给出自适应策略的一般说明. NIntegrate 的默认(主要)策略是 "GlobalAdaptive",这在"全局自适应策略"中做出解释. 与之相辅相成的是局部自适应策略,这在"局部自适应策略"中做出解释. 两种自适应策略都使用奇点处理机制,这在"奇点处理一节中有所说明.

蒙特卡洛策略的解释见 "原始蒙特卡洛与准蒙特卡洛策略""全局自适应蒙特卡洛与准蒙特卡洛策略".

"GlobalAdaptive"任何被积函数,自适应采样,基于规则
"LocalAdaptive"任何被积函数,自适应采样,基于规则
"DoubleExponential"任何被积函数,均匀采样
"Trapezoidal"任何被积函数,均匀采样
"MultiPeriodic"多维被积函数,均匀采样
"MonteCarlo"任何被积函数,均匀随机采样
"QuasiMonteCarlo"任何被积函数,均匀准随机采样
"AdaptiveMonteCarlo"任何被积函数,自适应随机采样
"AdaptiveQuasiMonteCarlo"任何被积函数,自适应准随机采样
"DoubleExponentialOscillatory"一维无限范围振荡被积函数
"ExtrapolatingOscillatory"一维无限范围振荡被积函数

NIntegrate 积分策略.

NIntegrate 对积分(或被积函数)的特殊类型使用一定的"预处理器"策略. 这在"Duffy 坐标策略"、"振荡策略"以及"柯西主值积分"等节有相应解释. 预处理器策略也处理被积函数的符号式预处理.

"SymbolicPreprocessing"全局预处理控制器
"EvenOddSubdivision"简化奇偶被积函数
"InterpolationPointsSubdivision"细分包含内插函数的被积函数
"OscillatorySelection"检测振荡被积函数并选择合适的方法
"SymbolicPiecewiseSubdivision"细分包含分段函数的被积函数
"UnitCubeRescaling"重新调整多维被积函数到单位立方体
"DuffyCoordinates"多维奇点消除变换
"PrincipalValue"等价于柯西主值的数值积分

NIntegrate 预处理器策略

自适应策略

自适应策略尽力集中处理被积函数不连续区域或存在某种形式的奇点区域改善. 自适应策略的不同在于将积分区域分成不相交子区域的方式. 每个子区域的积分估计共同构成了总积分估计.

自适应策略的基本假设是,给定积分规则 和被积函数 ,若一个积分区域 被分成了两个不连接的子区域 ,则 上的积分估计的和接近于实际积分 . 也就是说

且 (1) 将意味着误差估计 的和小于 的误差估计.

因此,一个自适应策略具有这些组件些组件[MalcSimp75]:

(1) 在一个区域上计算积分和误差估计的积分规则;

(2) 决定一个区域集合 的元素如何被分离的方法;

(3) 决定何时终止自适应策略算法的停止条件

全局自适应策略

全局自适应策略通过子区域的递归二分将最大误差估计分成两半,并计算每一半的积分和误差估计,从而达到积分估计所需要的精度和准确度目标.

NIntegrate 的全局自适应算法通过设定 Method 的选项值为 "GlobalAdaptive" 实现:
选项名称
默认值
MethodAutomatic在每一子区域上计算积分和误差估计所用的积分规则
"SingularityDepth"Automatic应用奇点处理器之前的递归对分的数目
"SingularityHandler"Automatic奇点处理器
"SymbolicProcessing"Automatic进行符号预处理的秒数

"GlobalAdaptive" 选项.

"GlobalAdaptive"NIntegrate 的默认积分策略. 它用于一维和多维积分. "GlobalAdaptive" 适用于迪卡尔积规则完全对称多维规则.

"GlobalAdaptive" 使用一种称为的数据结构来使区域的集合部分排序,最大的误差区域位于堆的顶部. 在算法的主循环中, 最大的误差区域在所估计的对误差负主要责任的层面上被一分为二.

可以说,该算法产生一个二叉树的叶子,节点是这些区域,节点/区域的孩子是二等分后得到的子区域.

经过对一个区域的二等分和在新的(子)区域进行随后的积分后,新的全局积分和全局误差估计被计算出来,它们是对所有作为二叉树的叶子的区域的积分和误差估计的和.

每个区域对于在每个层面上进行二分以便生成这个区域的次数都有记录. 如果一个区域的生成经过的二分次数过多,将应用一种奇平算法,见奇点处理.

"GlobalAdaptive" 在下述表达式为真时停止:

其中 pgag 精度和准确度的目标.

当一个区域的递归二分的次数超过了一定数值时(见 MinRecursion 和 MaxRecursion),或者是全局积分误差过于振荡时(见 "MaxErrorIncreases"),该策略也会终止.

理论和实践证据显示,全局自适应策略通常比局部自适应策略表现要优良[MalcSimp75][KrUeb98].

MinRecursion 和 MaxRecursion

递归二分的最小和最大深度由选项值 MinRecursionMaxRecursion 给出.

如果对于任何子区域,二分的次数在任何维度上都大于 MaxRecursion,则 "GlobalAdaptive" 的积分即终止.

MinRecursion 设置为一个正整数会强制在积分区域上的递归二分在被积函数被计算之前进行. 这样做的目的是为了确保被积函数的狭窄尖峰不被疏漏. (见诱骗误差估算器.)

对于多维积分,在 MinRecursion 的每一层递归的每一个维度上进行二分.

"MaxErrorIncreases"

由于在 "GlobalAdaptive" 中 (2) 预期是成立的,因此在对最大的误差区域进行二分并在新的部分上进行积分后,预期全局误差会降低. 换言之,全局误差预期会相对于积分步骤或多或少呈单调递减.

全局误差可能会由于积分规则的相位误差而振荡. 尽管如此,全局误差还是被假定为会在某些点处开始单调递减.

下面列出的是该假设可能错误的情形.

(i) 实际积分是零.

零积分:

(ii) 指定的工作精度对于指定精度目标来说密度不够.

工作精度的密度不够:

(iii) 积分条件太差 [KrUeb98]. 例如,原因可能是被积函数由复杂的表达式定义,或者是以数学问题近似解的形式给出(如微分方程或非线性代数方程等).

"GlobalAdaptive" 时刻追踪着策略对误差估计最大的区域进行二分后,总误差估计没有降低的次数. 当这个数目开始大于 "GlobalAdaptive" 选项的 "MaxErrorIncreases" 的值时,积分终止,并给出信息 (NIntegrate::eincr).

"MaxErrorIncreases" 默认值对于一维积分是400,对于多维积分是2000.

"MaxErrorIncreases" 取默认值时,下面的积分会导致出现这个信息 NIntegrate::eincr
增大 "MaxErrorIncreases" 的值后信息 NIntegrate::eincr 不再出现:
这个结果可以与准确值媲美:

全局自适应策略的示例

计算高斯Kronrod 横坐标、权重和误差权重:
这是一个函数的定义,应用的积分规则是在区间 {a,b} 上计算函数 f 的横坐标和权重:
定义一个简单的全局自适应算法,找到函数 f 在区间 {aArg,bArg} 上的积分,相对误差为 tol
定义一个被积函数:
前面定义的全局自适应策略给出下述结果:
这是准确结果:
相对误差在规定的公差范围内:

局部自适应策略

为了达到积分估计所需的精度和准确度目标,局部自适应策略将子区域递归分区为更小的不相交的区域,并对它们中的每一个进行积分和误差估计的计算.

NIntegrate 的局部自适应算法通过设定选项 Method 的值为 "LocalAdaptive" 实现 :
选项名称
默认值
MethodAutomatic在子区域上用于计算积分和误差估计的积分规则
"SingularityDepth"Automatic应用奇点处理器前进行递归二分的次数
"SingularityHandler"Automatic奇点处理器
"Partitioning"Automatic为了改善积分估计如何进行区域的划分
"InitialEstimateRelaxation"True为避免不必要的计算,试图调整初始积分估计的幅度
"SymbolicProcessing"Automatic进行符号预处理的秒数

"LocalAdaptive" 选项.

"GlobalAdaptive" 相同,"LocalAdaptive" 能够用于一维和多维积分. "LocalAdaptive" 适用于迪卡尔乘积规则完全对称多维规则.

"LocalAdaptive" 策略有一个初始化例程和递归例程(RR). RR 生成树叶,叶的节点是区域. 节点/区域的孩子是通过分区得到的子区域. RR 去一个区域作为自变量,并返回对这个区域的积分估计.

RR 使用一种积分规则计算区域自变量的积分和误差估计. 如果误差估计过大,RR 在分区得到的不相交的区域上调用自身.由这些递归调用返回的积分估计的和称为区域的积分估计.

RR 对继续递归与否做出决定,只知道在执行递归的区域上的积分和误差估计. (这就是为什么该策略被称为局部自适应性")

初始化例程计算计算在一个初始区域上的初始估计. 初始积分估计在 RR 的终止条件中被使用:如果一个区域的误差相对于初始积分估计显著,则该区域被分成不相交的区域,并在这些区域上调用 RR;如果误差不够显著则递归终止.

一个区域的误差估计,regionError 在满足下述条件时被认为是不显著的:

终止条件 (3) 将计算工作精度下的积分. 如果想要计算用户自行指定的精度和准确度目标下的积分估计,则使用下面的终止条件:

其中 eps 是使得 1+eps1 在工作精度下成立的最小数,pgag 代表用户自行指定的精度和准确度目标.

"LocalAdaptive" 的递归例程在下述情况时递归终止:

1.  区域的边界之间没有指定工作精度的数值;

2.  达到了最大的递归层;

3.  区域的误差不显著,例如条件 (4) 为真.

MinRecursion 和 MaxRecursion

选项 MinRecursionMaxRecursion"LocalAdaptive" 中的意义和功能与在 "GlobalAdaptive" 中相同. 见 MinRecursion 和 MaxRecursion.

"InitialEstimateRelaxation"

经过第一次递归后,一种较好的积分估计 将可用. 这个较好的估计与两个积分估计 作比较,其中 曾在初始过程中被积分规则用于给出积分估计 () 和误差估计 (). 如果

则在 (5) 中的积分估计 integralEst 能够通过公式

被增大也就是说,条件 (6) 的限制可以放松,因为 意味着该规则的积分估计与该规则的误差估计所预测的结果相比更加准确.

"Partitioning"

"LocalAdaptive" 通过选项 "Partitioning" 来指定如何对不满足 (7) 的区域进行分区. 对于一维积分来说,如果 "Partitioning" 被设置为 Automatic"LocalAdaptive" 将在重新调整的积分规则的采样点之间对一个区域进行分区. 通过这种方式,如果积分规则是封闭型的,每一个积分规则都可被重新使用. 如果 "Partitioning" 被给出一个整数列表 ,其中列表长度 等于积分变量的个数,则积分区域的每一个方向 被平均分为 等份. 如果 "Partitioning" 被赋予一个整数值 ,则所有的方向都被平均分为 等份.

考虑下述函数:
这是 "LocalAdaptive" 在默认分区设置时所使用的采样点. 可以看出每一个递归层的采样点都位于上一个递归层的采样点之间:
这是 "LocalAdaptive" 积分所用的采样点,其中将误差大的区域被分成三个子区域. 形成的这些模式清楚地表明,第一和第二层递归的每个区域的未来的下个递归层的三个子区域:
使用 "Partitioning" 选项的多维范例. 为了画出图形,被积分的第一个区域的采样点 被去除:

被积函数值的重复利用

在一维积分的默认分区设置下,"LocalAdaptive" 重复利用子区间终点处积分和误差估计不符合条件 (8) 的被积函数值.

通过 "LocalAdaptive" 进行积分的采样点. 变量 rulePoints 决定由"LocalAdaptive" 所使用的积分规则中点的个数:
积分中重复利用的点的百分比:

局部自适应策略的示例

计算 ClenshawCurtis 横坐标、权重和误差权重:
这是将积分规则应用到函数 f 的区间 {a,b} 上对一个函数的定义,其中横坐标和权重在上例中已经计算出来:
定义一个简单局部自适应算法,得到函数 f 在区间 {aArg,bArg} 上的积分,其相对误差为 tol
定义一个函数:
局部自适应策略给出结果:
这是准确结果:
相对误差在规定的公差范围内:

"GlobalAdaptive" 与 "LocalAdaptive"

通常全局自适应策略比局部自适应策略的表现要好,但在某些情况下,局部自适应策略更可靠和/或提供更好的性能.

"GlobalAdaptive""LocalAdaptive" 存在两个主要区别

1. "GlobalAdaptive" 当所有区域的误差之和达到精确度目标时即终止,而 "LocalAdaptive" 则当每一个区域的误差与某个积分估计相比足够小时终止.

2. 为了改善积分估计,"GlobalAdaptive" 对误差最大的区域进行二分,而 "LocalAdaptive" 则对所有误差不够小的区域进行分区.

对于多维积分,"GlobalAdaptive" 要快得多,因为 "LocalAdaptive" 沿每一个轴进行分区,因此区域的数量可能会发生组合爆炸.

关于全局自适应策略对于一维平滑被积函数较快的原因和方式在 [MalcSimp75] 中得到证明(和解释).

"LocalAdaptive""GlobalAdaptive" 速度较快,性能较好时,原因是 "LocalAdaptive" 的精确度目标终止条件和分区策略更适于该被积函数的特性. 另一个原因是 "LocalAdaptive" 重复利用已有采样点处被积函数值的能力. "GlobalAdaptive" 能够重复利用的积分值很少(每一个规则至多三个,对于默认一维规则高斯Kronrod 规则为零).

下一小节展示了 "GlobalAdaptive""LocalAdaptive" 性能的不同.

"GlobalAdaptive" 一般要优于 "LocalAdaptive"

在接下来的表中,展示了 "GlobalAdaptive" 要优于 "LocalAdaptive" 的大多数情形,并给出了两者的时间比率和被积函数被计算的次数. 所有被考虑的积分是在区间 上的一维积分,原因是 (1)对于多维积分,所期望的性能差异会更大,因为 "LocalAdaptive" 沿每一个轴进行分区,以及 (2) 在不同范围上的一维积分可以被重新调整成在区间 上的积分.

这是对一些函数、精确度目标、积分次数和积分规则的定义. 可以更改变量 integrationRule 以便比较同一积分规则时的运行情况. 最后一个函数是由 通过变量变换 ,将 变换为 而得到:
准确的积分值. 所有积分都在区间 上:
"GlobalAdaptive" 计时:
"LocalAdaptive" 计时:
"GlobalAdaptive" 函数计算:
"LocalAdaptive" 函数计算:
关于时间比率和被积函数计算的表:
带有积分误差的表. "GlobalAdaptive""LocalAdaptive" 都达到了要求的精确度目标:

奇点处理

NIntegrate 的自适应策略通过在积分区域边界的变量变换和用户指定的奇点或流形来使收敛数度加快. 自适应策略忽略在奇点处的被积函数值.

关于奇点的说明在"用户指定的奇点"中讨论.

带有变量变换的多维奇点处理要小心慎用,见 "IMT 多维奇点处理". 多维积分的坐标变换能够简化或消去奇点;见"Duffy 坐标变换法应对多维奇点处理".

关于 NIntegrate 如何忽略奇点的更多细节,请参见"忽略奇点".

关于柯西主值积分的计算,在"柯西主值积分" 中会介绍到.

用户指定的奇点

点类奇点(Point Singularities)

如果已知奇点的位置,则这些奇点可以在积分范围中指定,也可以通过选项 Exclusions 指定.

这个例子中的积分有两个奇点,分别在 处:
这是一个二维积分的例子,在 处有一个奇点:
这个例子中,积分的两个分别在 处的奇点通过 Exclusions 选项指定:
Exclusions 选项指定一个二维积分在 处的奇点:

曲线、曲面和超曲面的奇点

在曲线、曲面或超曲面上存在的奇点可以通过 Exclusions 选项用它们的方程进行指定. 这种奇点通常不能通过变量的范围来指定.

这个二维函数的奇点构成曲线
使用 Exclusions,积分得以迅速计算:
如果不给出奇点,NIntegrate 的收敛要慢得多:
这个例子中,奇点可以通过变量的范围指出. 如果 ,做到这一点是不可能的见下面的例子:

曲线、曲面和超曲面的奇点处理的示例

如果存在奇点的曲线、曲面或超曲面是隐式形式(例如能被描述成一个方程),可用函数 Boole来形成积分区域,而用奇点组成的曲线、曲面或超曲面来作为区域的边界.

二维函数的奇点构成曲线
使用 Boole,积分被迅速地计算出来:
二维函数的奇点构成曲线
使用 Boole,积分在一次被迅速地计算出来:
这是积分过程中的采样点:
这是一个函数,取出由奇点构成的曲线、曲面或超曲面,使用函数 Boole 来划分积分区域,使这些奇点出现在边界上:
定义一个三维函数:
这是三维函数的激愤,奇点是曲面
这是积分的采样点:

"SingularityHandler" 和 "SingularityDepth"

自适应策略通过区域二分法来改善积分估计. 如果一个自适应策略的子区域由选项 "SingularityDepth" 指定的二分次数得到,则可以断定子区域存在一个奇点. 这样,在子区域上的积分可由 "SingularityHandler" 指定的奇点处理器来完成.

选项名称
默认值
"SingularityDepth"Automatic应用奇点处理器之前进行递归二分的次数
"SingularityHandler"Automatic奇点处理器

"GlobalAdaptive""LocalAdaptive" 的奇点处理选项.

如果在一个给定积分区域的边界上有一个可积的奇点,二分很容易在收敛前重复进行至 MaxRecursion. 为了应对这种情况,NIntegrate 的自适应策略使用变量变换 (IMT, "DoubleExponential", SidiSin) 来使积分收敛加速,或者使用区域变换 (Duffy 坐标) 来使奇点的秩序放宽. 变量变换奇点处理器的理论背景由 EulerMaclaurin 公式[DavRab84] 给出.

IMT 变量变换的使用

IMT 变量变换是一种由 Iri、Moriguti 和 Takasawa 提出的在正交变换方法中进行的变量变换,在文献中被称作 IMT 规则 [DavRab84][IriMorTak70]. IMT 规则基于一种独立变量变换的思想,这种变换方式使得新的被积函数的所有导数在积分区间的终点处消失. 然后在新的被积函数上应用梯形规则,在适当的条件下,可能得到高准确度的结果 [IriMorTak70][Mori74].

这是一个使用 IMT 变量变换进行奇点处理的数值积分:
选项名称
默认值
"TuningParameters"10在变换公式 中作为调节参数的一对数值 {a,p};如果只给出 a 的值,则被解释为 {a,1}

IMT 奇点处理器选项.

NIntegrate 的自适应策略仅使用 NIntegrate 规则的变换. 当认定一个区域可能存在奇点后,IMT 变换被应用到它的被积函数. 积分继续进行,尽管不是使用梯形规则,而是使用与变换前相同的积分规则. ("DoubleExponential" 时的积分处理转成一种梯形积分规则.)

另外,NIntegrate 的自适应策略使用原始 IMT 变换的一种变体,变换后的被积函数仅在一个端点处消失.

定义 IMT 变换

参数 ap 被称作调节参数[MurIri82].

IMT 变换的导数的极限是0:
这是 IMT 变换的的图形:

从上图得出,变换后的采样点在 0 附近更加密集. 这意味着如果被积函数的奇点在 0 处,它可以被更有效地采样,这是由于积分规则采样点的较大一部分将对积分规则的积分估计贡献大的被积函数值.

由于对于任意给定的工作精度, 周围的数值比在 周围密集得多,经过一次区域二分后,NIntegrate 的自适应策略扭转了二分后区间的右端点所在的子区域的二分变量. 这可以从下面的图形中看出来:

对于应用了 IMT 变量变换的区域的子区域,将不再应用其它的奇点处理器.

IMT 变换实例

考虑在区间 上的函数 ,这个函数在0处有奇点:
假设积分通过 "GlobalAdaptive" 结合使用奇点处理器 IMT 及奇点深度4来完成. 经过四次二分,"GlobalAdaptive" 将有一个边界为 、且包含奇点端点的区域. 对于这个区域,IMT 变量变换的边界改为 ,并将其被积函数变成下面的函数:
这是新的被积函数的图形:

奇点消失了!

然而也有一些采样点在奇点的端点过于密集,由于 IMT 变换的原因,对与奇点重合的采样点要特别注意. NIntegrate 跳过在奇点处的计算,见"忽略奇点".

考虑高斯-Kronrod 规则的采样点和权重:
随后是在区域 上的高斯Kronrod 采样点和比例调整后的导数:
这是积分估计:
这是经过 IMT 变换后的采样点和导数:
这是采用 IMT 变换的积分估计:
采用 IMT 变量变换的估计非常接近于准确值:

双指数正交法的使用

当自适应策略使用 IMT 变量变换时,在IMT 变换后的区域上的积分规则并没有变化. 与此形成对照的做法是,使用变量变换的同时,在可能存在奇点的区域上运用不同的积分规则. (这是 IMT 规则的精髓[DavRab84].) 这正是使用双指数正交法所采用的方法双指数正交法使用梯形规则.

NIntegrate 仅对一维积分能够使用双指数正交法来处理奇点.

这是一个数值积分,使用双指数正交法处理奇点:

IMT、"DoubleExponential" 以及不对一维积分进行奇点处理等方法的比较

两种奇点处理器 (IMT"DoubleExponential") 适用于通过过多次二分后得到的区域 (在"SingularityDepth" 中指出).

两者的主要区别是 IMT 不改变对应用 IMT 的区域上的积分估计的积分规则IMT 只是一种变量变换. 另一方面,"DoubleExponential" 既使用变量变换又使用不同的积分规则梯形规则来计算所应用的区域上的积分估计. 换言之,奇点处理器 "DoubleExponential" 将积分任命给双指数策略中介绍的双指数正交法.

因此,应用 IMT 奇点处理器的区域仍然受到自适应积分策略的二分法的控制. 因此,在达到精确度目标之前,最后一次二分之前的被积函数值的计算都将被废弃不用. 另一方面,应用 "DoubleExponential" 奇点处理器的区域将不再进行二分. 由 "DoubleExponential" 所使用的梯形规则求积法所进行积分估计的计算过程中,随着每一步地进行,采样点的数目逐渐增加,并对上一步中所有采样点的被积函数值进行重复利用.

因此,如果被积函数在含有奇点的区域上具有明显的解析性质(例如,被积函数及其导数相对于积分变量没有快速或者突然的变化),则 "DoubleExponential" 奇点处理器将比 IMT 奇点处理器要快得多. 对于被积函数在应用 "DoubleExponential" 奇点处理器的区域不具有解析性质的情形,或者是被积函数的双变换收敛过慢的情形,最好是转为 IMT 奇点处理器. 这可以通过将选项 "SingularityHandler" 设置为 Automatic 来实现.

下面的表中将不同深度二分时的奇点处理器 IMT"DoubleExponential" 以及 Automatic 进行比较.

加载一个软件包,定义分析函数 NIntegrateProfile 以给出采样点的数量和执行数值积分命令所需的时间:
关于一个明显具有解析性质的被积函数 的表格,使用 "DoubleExponential" 奇点处理器可以轻松计算出来:
一个关于被积函数 的表格,该函数不存在奇点,并有一个几乎不连续的导数(也就是说,不具有明显的解析性质). Automatic 奇点处理器开始使用 "DoubleExponential",然后转成 IMT
一个关于被积函数 的表格,其中 Automatic 奇点处理器开始使用"DoubleExponential",然后转成 IMT

IMT 多维奇点处理

当用于多维积分时,IMT 奇点处理器仅当奇点沿着一个坐标轴时才能够加速积分进程. 当奇点落在积分区域的一角时,使用 IMT 时会适得其反. 先前定义的 NIntegrateProfile 被用在下面的例子中.

一个仅沿 轴有一个奇点的被积函数的计算次数与计时. 默认(自动)奇点处理器选择对经过默认的(4次)二分后得到的区域应用 IMT:
一个仅沿 轴有一个奇点的被积函数,在不应用任何奇点处理器时的计算次数与计时:
被积函数的奇点落在积分区域的一角时,其计算次数和计时. 默认(自动)奇点处理器选择对经过默认的(4次)二分后得到的区域应用奇点处理器 DuffyCoordinates
被积函数的奇点落在积分区域的一角时,其计算次数和计时. 对经过默认的(4次)二分后得到的区域应用 IMT:
被积函数的奇点落在积分区域的一角时,在不应用任何奇点处理器时的计算次数和计时:

处理多维奇点的Duffy 坐标变换法

Duffy 坐标变换法是一种将被积函数在一个正方形、立方体或超立方体上进行变换的技术,通过这种变换使得在角落处的奇点落在新的被积函数的一条线上,从而可能使积分变容易.

下面的积分使用 Duffy 坐标变换法:
下面的积分不使用 Duffy 坐标变换法:

NIntegrate"GlobalAdaptive""LocalAdaptive" 仅在积分区域的角落处应用 Duffy 坐标变换技术.

当一个多维积分的奇点在某点出现时,变量的耦合会使得用于一维积分的奇点变量的变换起到相反的效果. 由 Duffy 在 [Duffy82] 中提出的具有几何性质的变量变换,对变量做出改变,将在积分区域的角落上的奇点用平面上的一个叫的点来替换.

如果 是积分的维数,且有 ,则对于下述类型的奇点, Duffy 坐标变换法是一种合适的技术(仍参见 [Duffy82]):

1.  , , ;

2.  , , ;

3.  , , , .

例如下述积分

如果积分区域 根据规则 变成 , 其雅可比式(Jacobian )为 ,积分变成

最后一个积分根本不存在奇点!

现在考虑积分

它等价于和式

该和式的第一个积分变形如 (9);而对于第二个积分,通过 的变换存在雅可比量 ,这不会理想地进行项的消去. 幸运的是,积分次序的改变:

是第二个积分改造成适合 (10) 的变形:

(在方程 (3) 的第二个积分中,变量被置换,这并不需要证明数学上的等价性,但进行积分的计算时速度更快.)

于是积分 (11) 可以被改写为一个没有奇点的积分:

如果积分变量在 (12) 中没有被置换,则积分 (13) 将被改写成下述形式

这是一个更加复杂的积分,被积函数在两个轴上都不是简单的形式,因此其积分的计算比前一个更困难.

这是这个较简单的积分的采样点个数:
这是较复杂的积分的采样点个数:

NIntegrate 使用的是上例中所用技术在任意维度上的一种归纳. (在 [Duffy82] 中仅介绍了三维的情形.) 一个与归纳描述结合使用的实例在下面给出.

这个表对 所用的不同奇点处理方法进行比较.(这里采用的是先前定义的分析函数 NIntegrateProfile.)

Duffy 坐标变换策略

当 Duffy 坐标变换适用时,如果 Duffy 坐标变换发生在实际积分开始之前,则可以较快地得到一个数值积分结果. 然而提前进行这种变换需要知道奇点出现在积分区域的哪个(些)角上. NIntegrate 中的 "DuffyCoordinates" 策略为这种积分前进行的变换提供了方便.

这个例中的被积函数在积分区域的两个角上存在奇点:
选项名称
默认值
Method{"GlobalAdaptive","SingularityDepth"->}使积分在 Duffy 坐标变换之后进行的策略
"Corners"All指定在哪个(些)角上应用 Duffy 坐标变换的向量或向量列表;向量中的元素是 0 或 1;每个向量的长度与积分的维数相同

"DuffyCoordinates" 选项.

"DuffyCoordinates" 所作的第一件事是进行积分坐标的调整,使之出现在一个单位超立方体(也可以是正方形或者立方体). 如果只指定了一个角,"DuffyCoordinates"前面介绍的一样应用 Duffy 坐标变换. 如果指定的角多于一个,则前一步中的单位超立方体被分成不相交的立方体,边长为原来的一半. 考虑在这些不相交的立方体上的积分. Duffy 坐标变换应用于顶点被指定为奇点的立方体. 其余的立方体被变形成单位立方体上的积分. 由于到了这一步,所有的积分的积分区域都是单位立方体,将这些积分相加,并将结果赋与 NIntegrate,它的 Method 选项与与赋给 "DuffyCoordinates"Method 选项相同.

"DuffyCoordinates" 使用的实际积分可以通过 NIntegrate`DuffyCoordinatesIntegrand 得到,它具有与 NIntegrate 相同的参数.

这个例子中,"DuffyCoordinates" 所用的被积函数是一个三维的函数,并且这个函数在积分区域的一个角上存在奇点:
这个例子中,"DuffyCoordinates" 所用的被积函数是一个二维的函数,并且这个函数在积分区域的两个角上存在奇点:

对于"Duffy 坐标变换应对多维奇点处理"中介绍的积分函数类型"DuffyCoordinates" 可能可以大大提高速度

使用 "DuffyCoordinates" 的积分.
如果使用默认的 NIntegrate 选项设置,积分速度将大大慢于前面的例子.
这是又一个通过 "DuffyCoordinates" 使速度加快的例子:
如果使用默认的 NIntegrate 选项设置,积分速度将大大慢于前面的例子.

Duffy 坐标变换的归纳和示例

关于 Duffy 坐标变换的理论,见"Duffy 坐标变换应对多维奇点处理".

执行基于以下两个定理.

定理1: 一个 维立方体能够被分成 个不相交的几何上等价的 维金字塔 (其基底为 维立方体),并且具有相同的顶点.

证明: 假设立方体的边长为 ,一个顶点在原点,所有其它顶点的坐标为 . 考虑该 维立方体的面,其中 . 它们的个数正好是 ,并且原点不在这些面上. 每一个 面可以形成一个顶点在原点的金字塔. 证明完毕.

这是描述该定理的三维图形:

如果 个轴分别被记为 ,构成这个金字塔的面 可以被描述成. 令 表示 向左旋转 次后的排列(即对 应用 RotateLeft). 那么有下面的定理成立:

定理2: 对于任何在单位立方体上的任何积分,有下述性质成立:

证明: 第一个性质可由定理1得到. 第二个性质则仅是由金字塔到立方体的变量变换.

这个函数给出对一个超立方体进行变形的规则及雅可比量,该立方体的一条边在一个区域上:
下面是一个单位正方形在无限区域缩放的例子:
这个函数计算由 Duffy 坐标变换技术得到的积分,其奇点在原点:
这个函数计算由 Duffy 坐标变换技术得到的积分,其中奇点出现在超立方体的一个指定的角上:
这是一个符号化的例子:
又一个一个符号化的例子:
下面是一个算例:
使用 Duffy 坐标变换比不进行奇点处理要快得多 (见下例):
不进行奇点处理的积分.
当然,NIntegrate 的内部执行给出表现相似的结果:

忽略奇点

另一个处理奇点的方法就是无视它的存在. NIntegrate 忽略与奇点重合的采样点.

考虑下面的积分,该积分在 1 处有一个奇点.

当积分变量向 1 靠近时,被积函数趋向于 .

这是被积函数的图形:
NIntegrate 给出接近准确值的结果:
当增加 MaxRecursion 时实现了收敛:

使用默认选项时,NIntegrate 在1处有一个采样点.

验证得知 NIntegrate 确实有一个采样点为1:

但对于 NIntegrate[Log[(1-x)2],{x,0,2}],评估监视器没有选择采样点1.

属于区间 的采样点:

换言之,在1处的奇点被忽略了. 忽略奇点等价与被积函数在奇点采样点的值为零.

注意如果奇点在变量范围中被指名,积分则容易积出. 下面是 NIntegrate 在给或不给出奇点说明时采样点的个数和计时.

具有奇点说明的积分:
忽略奇点的积分:

忽略奇点的一个更加有趣的例子是在被积函数的分母上使用贝塞尔函数.

具有几个(5个)可积奇点的积分:

结果可以使用 NIntegrate 检查,其奇点范围说明用 BesselJ[2,x] 的零值给出(见 BesselJZero).

贝塞尔零值作为奇点的积分:

不用说,最后一个积分要求计算 BesselJ 的零值. 前一积分是在不进行任何被积函数的分析的情况下进行积分.

对于振荡型的被积函数,无视奇点可能行不通.

例如,这两个积分是等价的:
NIntegrate 可以计算第一个:
NIntegrate 不能计算第二个:

然而,如果被积函数在它的奇点附近是单调的,或更确切地说,如果它能够被一个单调可积函数主导,可以证明通过忽略奇点实现了收敛.

关于忽略奇点的理论证明和实践推荐,参见 [DavRab65IS] 和 [DavRab84].

自动奇点处理

一维积分

对于一维积分,当选项 "SingularityHandler" 被设置成 Automatic 时,"DoubleExponential" 用于多次分区后得到的区域上,分区次数由 "SingularityDepth" 指定. 正如前面所解释的,只要 "DoubleExponential" 奇点处理器在这个区域上工作,这个区域将不再继续进行分区. 如果由 "DoubleExponential" 所计算的误差估计不按照双指数求积法的理论所预测的方式进行,对该区域的奇点处理则转为 IMT.

正如"收敛速率"所解释的,可以预见该误差与双指数采样点个数存在下面的依赖关系:

其中 为正常数. 考虑两个连续双指数求积法计算的相对误差 ,其中采样点的个数分别为 ,且 . 假设 , , 及 ,能够预见到

"DoubleExponential"IMT 的转换发生在这些情况下:

(i) 区域的误差估计大于区域积分估计的绝对值 (因此相对误差不小于1);

(ii) 不等式 (2) 在两种不同情形中不成立;

(iii)通过双指数变换计算得到的被积函数值的衰变速度不够快.

这是一个由 "DoubleExponential"IMT 奇点处理进行转换的例子. 在图形中,被积函数按照 坐标的顺序在 坐标上采样. 采样点在 上的布局表明由高斯求积法 () 到双指数求积法 () ,再到后来通过 IMT 变量变换 () 再次被高斯求积法替代的变化:

多维积分

对于多维积分,当选项 "SingularityHandler" 被设置成 Automatic"DuffyCoordinates"IMT 都被使用.

一个区域需要满足下述条件才能应用 "DuffyCoordinates"

该区域通过沿每一个轴线进行 "SingularityDepth" 次的二分 (或分区) 得到;

该区域是一个原始积分区域的一角(指定的积分区域可以通过分段处理或用户指定奇点进一步划分积分区域).

一个区域需要满足下述条件才能应用 IMT

该区域通过沿一个主要轴线进行 "SingularityDepth" 次的二分 (或分区) 得到;

该区域不是一个角区域,并且在一个原始积分区域的一条边上.

换言之,IMT 应用于经过次数为 "SingularityDepth"的分区后得到的、但不满足 "DuffyCoordinates" 的自动应用条件的区域.

如果奇点沿着一个轴线,则 IMT 有效. 使用 IMT 处理点式奇点可能会适得其反.

二维积分 的采样点,采用的是 Automatic (左) 和 "DuffyCoordinates" (右) 的奇点处理. 可以看出自动奇点处理使用的点数几乎是 "DuffyCoordinates" 所用点数的两倍. 为了说明奇点处理器的效果,它们应用于两次二分后:
积分 的计时,其中奇点处理器是 Automatic"DuffyCoordinates"IMT,或者不用奇点处理. 积分在 处有奇点:
积分 的计时,奇点沿 轴,奇点处理器是 Automatic"DuffyCoordinates"IMT,或者不用奇点处理:

柯西主值积分

为计算一个积分的柯西主值,NIntegrate 使用策略 PrincipalValue.

奇点在2的柯西主值积分:

NIntegrate 中,PrincipalValue 使用由 Method 选项指定的策略,直接作用于没有困难的区域,并且通过对指定的奇点进行对称配对,从而利用正负值能够相互抵消的优势.

选项名称
默认值
MethodAutomatic用于计算子区域上估计的方法指定
SingularPointIntegrationRadiusAutomatic在范围指定时对应于奇点 的一个数值 ϵ1 或一个数值列表 {ϵ1,ϵ2,,ϵn};对于每一对 ,构成一个形如 的积分

"PrincipalValue" 选项.

于是指定

的形式进行计算,其中使用 Method->methodspecNIntegrate 计算每一个积分. 如果 ϵ 没有显式给出,则根据 b-ac-b 的差选出一个值. 选项 SingularPointIntegrationRadius 能够取与奇点个数相等的一列数值. 关于攻势的推导,见 [DavRab84].

找到 的柯西主值:
这是 的柯西主值. 注意有两个奇点需要指定:
奇点可以使用 Exclusions 选项来指定:
检查该值. 如果一切正确的话,结果应该是0:

应该注意的是,奇点的定位必须准确. 由于该算法将奇点两边的点配对,如果奇点位置稍微有所偏差,则抵消将在极点附近不够理想,从而导致即使 NIntegrate 收敛,其结果也有可能出现明显的差错.

采样点的可视化

考虑下式的主值的计算:

下面的例子展示了两种将采样点可视化的方法. 第一个显示了所用的采样点. 为了进行柯西主值积分,被积函数被改变,有可能更希望看到的是在原始被积函数的点处的计算. 这在第二个例子中说明.

NIntegrate 使用的采样点. 在区间 上没有点,原因是 PrincipalValue 积分 ,在 上有采样点:
定义一个函数,将赋给被积函数的自变量的值累计:
在这些点上进行的被积函数的计算. 注意在区间 上的对称布局:

双指数策略

双指数求积法在变量变换后应用梯形规则. 双指数求积法由 Mori 和 Takahasi 在 1974 年提出,它的灵感来自于所谓的IMT 规则和 TANH 规则. 这种变换被冠以"双指数"之名是由于在被积函数的变量到达积分区域的尾部时,它的导数将呈双指数递减.

NIntegrate 的双指数算法通过令 Method 选项值为 "DoubleExponential" 实现:
选项名称
默认值
"ExtraPrecision"50内部使用的最大额外精度
"SymbolicProcessing"Automatic符号预处理所用的秒数

"DoubleExponential" 选项.

双指数策略可用于一维和多维积分. 当应用于多维积分时,它使用梯形规则的迪卡尔积.

一个双指数变换 将积分

变形为

其中 可以为有限,半无穷 (),或者是无穷(). 被积函数 上必须解析,并且可能在一个或两个端点上有奇点.

变换后的被积函数双指数递减,也就是,当 时, .

函数 具有解析性. 像 (14) 中这样的解析型的被积函数的积分,梯形规则是已知的最佳规则 [Mori74].

对于不同类型的积分区域,变形如下:

其中 取有限值.

梯形规则应用于 (15):

对于足够大的 ,(16) 中的各项呈双指数衰减. 因此,在 (17) 的求和中,对总和贡献过小的项被删去. (对于局部自适应策略,一种与 (18) 相似的判据被使用. 另请参见下面的双指数示例.)

"DoubleExponential" 策略使用双指数求积法.

"DoubleExponential" 策略对于解析型被积函数效果最好;见"双指数和高斯求积法的比较".

"DoubleExponential" 使用多维积分中的双指数求积法的笛卡尔积(Cartesian product).

笛卡尔双指数求积法:

与其它的笛卡尔积的规则一样,如果 "DoubleExponential" 用于大于三的维数,由于组合爆炸的缘故,速度可能会非常慢.

下面的图形描述了 "DoubleExponential" 多维积分的笛卡尔积特征:

双指数求积法可用于自适应策略的奇点处理;见"奇点处理".

MinRecursion 和 MaxRecursion

正如"MinRecursion 和 MaxRecursion"中所介绍的,选项 MinRecursion"GlobalAdaptive" 中与在 "LocalAdaptive" 中具有相同的意义和功能. "DoubleExponential"MaxRecursion 限制了体形求积法估计被改善的次数;见"双指数求积法的示例".

双指数与高斯求积法的比较

"DoubleExponential" 策略对于解析型被积函数效果最好. 例如,下面的积分用 "DoubleExponential" 完成,其速度是高斯求积法的三倍 (使用一种全局自适应算法).

使用 "DoubleExponential" 积分.
使用高斯求积法积分. ("GlobalAdaptive" 作为 NIntegrate 的默认策略,默认采用一种高斯Kronrod积分规则,其中高斯点为 个,Kronrod 点为 个points.)

由于 "DoubleExponential" 相对于所计算点的数目呈双指数收敛,精确度目标的提高使由 "DoubleExponential" 所完成的工作稍微有所增加. 这通过两个积分 来说明. 各表项表明了计算的误差和次数.

采用的双指数求积法与高斯求积法. 精确度目标的提高不能改变 "DoubleExponential" 所用的采样点的个数:
采用的双指数求积法与高斯求积法. 精确度目标的提高不能改变 "DoubleExponential" 所用的采样点的个数. (积分在不进行符号预处理的情况下完成.)

另一方面,对于非解析被积函数,"DoubleExponential" 速度很慢,使用一种利用高斯求积法的全局自适应策略可以轻松求解起点.

"DoubleExponential" 需要对被积函数进行大于1000次的计算才能计算出该非解析被积函数的积分:
对于该积分高斯求积法要快得多:

此外,"DoubleExponential" 可能会被具有几乎不连续导数(也就是不非常具有解析性质)的被积函数拖慢.

这是一个不非常具有解析性质的被积函数:
同样,高斯求积法要快得多:
这是被积函数 及其导数的图形:

收敛速度

这一节说明的是,双指数求积法的渐近误差以所用计算点的个数 的形式表示是

其中 为正常数.

定义一个双指数积分函数,返回积分估计和所用点的个数:
定义一个函数:
这是准确的积分值:
找到梯形规则在一定范围的步骤内误差和计算点的个数:
通过误差的对数进行 的拟和;见 (19):
这是拟和后的函数. 和项 30.48 仅是一个转换参数:
可以看到误差与模型 (20) 拟和:

双指数求积法的示例

下面是一个进行了有限区域变量变换(先前的变换 (21))的双指数求积法的示例.

定义一个函数,对一个变换后的被积函数应用梯形规则. 该函数使用 (22),并且它重复使用由一个两倍大的步骤计算得到的积分估计:
一个简单的双指数策略的定义,找到函数 f 在有限区间 {a,b} 上的积分,相对误差为 tol
定义一个奇点在0的函数:
这是由双指数策略得到的积分估计:
这是准确结果:

两个结果相同.

定义一个振荡函数:
这是由双指数策略给出的积分估计:
这是准确结果:
这是机器精度下的准确结果:
相对误差在规定的公差范围内:

"Trapezoidal" 策略

当积分区间正好是一个周期时,"Trapezoidal" 策略给出解析型周期被积函数的最佳收敛.

选项名称
默认值
"ExtraPrecision"50内部使用的最大额外精度
"SymbolicProcessing"Automatic符号预处理所用的秒数

"Trapezoidal" 选项.

"Trapezoidal" 的选项与 "DoubleExponential" 相同. 如果积分范围是无穷或半无穷,"Trapezoidal" 就变成 "DoubleExponential".

关于周期性函数积分(使用梯形求积法)的理论背景、实例以及注解见 [Weideman2002].

这个表显示了积分 ,当 "GlobalAdaptive""Trapezoidal" 的参数 分别取不同值时采样点的个数:

示例

这个函数用指定的点给出了梯形求积法的积分估计:
这个函数通过将原来的点之间的点作为采样点,改善了梯形求积法的积分估计:
这个函数是前面一个的接口:
定义一个(贝塞尔)函数:
这是梯形求积法估计:
这是准确值:
相对误差在规定的公差范围内:

振荡策略

NIntegrate 包括振荡积分的多种方法. 积分规则 "LevinRule" 适用于一大类一维和多维振荡积分. 它与"GlobalAdaptive" 等基于自适应规则的策略结合使用.

NIntegrate 也包括多个专用于特定类型一维振荡积分的积分策略. 一般地,这些算法或者仅用于有限区域上的积分,或者仅用于无限区域上的积分. 对于有限区域上的振荡积分,NIntegrate 使用被积函数的切比雪夫展开和全局自适应积分策略. 对于无限区域上的振荡积分,NIntegrate 使用双指数算法的修正,或在被积函数零点之间的区域上进行积分序列的序列加和加速.

下面的示例在不同子区域上使用各种不同的专用振荡策略:

NIntegrate 自动检测到振荡型(一维)被积函数,并根据被积函数的范围和特定振荡形式选择要用的策略或规则.

有此处所介绍的专用策略所处理的积分具有下述形式:

其中振荡内核 的形式如下:

1.   为有限时,, ,

2.   为无穷或半无穷时,, .

在这些振荡内核形式中, 为实数常数, 为正整数.

关于该积分规则所处理的振荡被积函数的更一般形式描述,请参见 "LevinRule".

有限区域振荡积分

修正 ClenshawCurtis 求积法 ([PiesBrand75][PiesBrand84]) 针对的是下述形式的有限区域一维积分:

其中 , , , , 为有限实数.

修正 ClenshawCurtis 求积规则使用一个单一多项式通过切比雪夫多项式展开来近似. 由于正弦和余弦函数的切比雪夫多项式的正交性,从而导致计算的简化. 修正 ClenshawCurtis 求积规则与 "GlobalAdaptive" 策略联合使用. 对于平滑的 ,通常,修正 ClenshawCurtis 求积法具有其它处理振荡积分的方法 (例如 Filon 求积法和被积函数零值间的多面板积分) 所无法比拟的优越性[KrUeb98] .

修正 ClenshawCurtis 求积法非常适用于形如(23)的高度振荡的积分. 例如,对于 ,修正 ClenshawCurtis 求积法对被积函数的计算次数都小于100.

对于振荡缓慢的内核,修正 ClenshawCurtis 求积法进行被积函数计算的次数:
对于振荡缓慢的内核,修正 ClenshawCurtis 求积法的计时和积分估计:
对于高度振荡的内核,修正 ClenshawCurtis 求积法进行被积函数计算的次数:
对于高度振荡的内核,修正 ClenshawCurtis 求积法的计时和积分估计:

另一方面,如果不进行符号预处理,默认的 NIntegrate 方法使用高斯Kronrod 规则的"GlobalAdaptive" 策略对 进行几千次的计算,对 则无法积分.

对于振荡缓慢的内核,高斯求积法进行被积函数计算的次数:
对于振荡缓慢的内核,高斯求积法的计时和积分估计:
对于高度振荡的内核,高斯求积法进行被积函数计算的次数:
对于高度振荡的内核,高斯求积法的计时和积分估计:

外推振荡策略

NIntegrate 策略 "ExtrapolatingOscillatory" 针对的是一维无穷区域上的振荡积分. 该策略运用序列和的序列收敛加速,该序列和由被积函数的两个连续零点之间区域上的积分组成.

一个运用 "ExtrapolatingOscillatory" 进行积分的例子:
选项名称
默认值
MethodGlobalAdaptive如果 ExtrapolatingOscillatory 失败,在零值之间积分所用的积分策略
"SymbolicProcessing"Automatic符号预处理所用的秒数

考虑下面的积分

其中函数 是振荡内核,函数 是平滑的. 令 的零值,由(有限)积分下限开始依次枚举,也就是说,不等式 成立. 如果积分 (24) 收敛,则序列

亦收敛. 序列 (25) 中的元素是序列的部分和

往往可以通过某种收敛加速技术,使用较少的元素就可将序列 (26) 的极限计算出来.

"ExtrapolatingOscillatory" 策略对 (27) 中的积分使用结合Wynn 外推法的 NSum. (28) 中的每一个积分通过NIntegrate 计算,而不使用振荡法.

"ExtrapolatingOscillatory" 策略将它的算法应用于 (29) 中具有如下形式的的振荡内核: 或者 ,其中 是实数常数.

示例

下面的例子介绍了 "ExtrapolatingOscillatory" 策略是如何工作的.

定义一个将在区间 上被积分的振荡函数. 该振荡函数 的零值是
这是该振荡函数在区间 上的图形:
定义一个函数使积分在两个连续零值之间进行. 该振荡函数 的零值是
通过序列收敛加速(外推)计算的积分估计:
这是准确的积分值:
该积分估计非常接近于准确值:
使用 "ExtrapolatingOscillatory" 策略再次验证:
通过 "ExtrapolatingOscillatory" 所作的积分估计非常接近于准确值:

双指数振荡积分

策略 "DoubleExponentialOscillatory" 针对的是缓慢衰减的振荡积分,积分区域是一维无限区域,被积函数形如 , ,其中 为积分变量, 为常数.

使用 "DoubleExponentialOscillatory" 进行积分:
选项名称
默认值
MethodNone"DoubleExponentialOscillatory" 失败时所用的积分策略
"TuningParameters"Automatic误差估计的调节参数
"SymbolicProcessing"Automatic符号预处理所用的秒数

"DoubleExponentialOscillatory" 选项.

"DoubleExponentialOscillatory" 基于 "DoubleExponential" 策略,但 "DoubleExponentialOscillatory" 使用的变换不是以双指数形式到达积分区间的两端,而是以双指数形式到达 的零值. 关于该算法的理论基础和性质在 [OouraMori91],[OouraMori99],[MoriOoura93] 中有解释. "DoubleExponentialOscillatory" 的运行示例使用的是[OouraMori99] 中的公式和积分器设计.

"DoubleExponentialOscillatory" 的算法将用正弦积分

解释. 考虑下述变换

其中 是常数,满足

选择参数 以满足

(选自 [OouraMori99]).

将变换 (30) 应用于 (31) 得到

注意 在正弦项中消失. 将具有相同网格尺寸 的梯形公式应用于 (32) 得到

这被近似于截断级数的和(truncated series sum)

选择合适的 以满足

正如 (33) 中所看到的,被积函数在 取大的负值时呈双指数形式衰减. 尽管"双指数策略"一节中,(34) 的双指数变换使得被积函数在 取大的正值时呈双指数形式衰减,(35) 的变换不会在 取大的正值时衰减,而是使得采样点在 取大的正值时,以双指数形式靠近 的零点. 并且有

如在 [OouraMori99] 中所解释的,由于 在它的零值附近是线性的,当 接近于 的零值时,被积函数呈双指数形式衰减. 这也是为什么 (36) 被认为是双指数公式.

相对误差假定满足:

在 [OouraMori99] 中, 的建议值是 .

由于 公式不能成为渐进式的,"DoubleExponentialOscillatory" (如在 [OouraMori99] 中所提议的) 对于 之间的积分估计使用不同的 . 如果所需的相对误差为 ,积分步骤如下:

1. 选择 使满足

并且计算 时的 (37). 结果记为 .

2. 然后,令 ,并且计算 时的 (38). 结果记为 . 假定第一步积分中的相对误差 . 由 (39) ,因此,若

成立,其中 为鲁棒性因素 (默认值为10) "DoubleExponentialOscillatory" 退出,结果为 .

3. 如果 (40) 不成立,则计算

并计算 时的 (41). 若

"DoubleExponentialOscillatory" 退出,结果为 .

4. 如果 (42) 不成立,则计算

并计算 时的 (43). 令结果为 . 若

不成立,"DoubleExponentialOscillatory" 发出讯息 NIntegrate::deoncon. 如果"DoubleExponentialOscillatory" 方法选项的值为 None,则返回 . 否则"DoubleExponentialOscillatory" 将与 "DoubleExponentialOscillatory" 方法选项合用,返回 NIntegrate 的结果.

对于余弦积分

对应于 (44) 的变换为

广义积分

这里是正规化发散积分 的符号计算:
"DoubleExponentialOscillatory" 计算广义意义上的上述非正规化积分:

用在发散傅立叶型积分中的关于 "DoubleExponentialOscillatory" 的更多性质,见 [MoriOoura93].

非代数被乘数

一个振荡积分的符号积分:
如果振荡内核被乘以一个非代数函数,"DoubleExponentialOscillatory" 仍能给出一个好的结果:
被积函数的图形和它的振荡内核:

原始蒙特卡洛和准蒙特卡洛策略

原始蒙特卡洛算法对给定的积分通过取在积分区域上均匀分布的随机点处被积函数值的平均值进行估计. 点的数量递增,直到估计标准偏差小到足以满足指定精度或准确度的目标. 阿蒙特卡洛算法叫做准蒙特卡罗算法,如果一种蒙特卡洛算法使用的是均匀分布的、确定性产生的点序列,而不是均匀分布的随机点,它就叫做准蒙特卡洛算法.

这是一个原始蒙特卡洛积分:
这是一个原始准蒙特卡洛积分:
选项名称
默认值
Method"MonteCarloRule"蒙特卡洛积分规则规范
MaxPoints50000采样点的最大数目
"RandomSeed"Automatic来重置随机数发生器的种子
"Partitioning"1沿各轴划分积分区域
"SymbolicProcessing"0符号预处理所用的秒数

"MonteCarlo" 选项.

选项名称
默认值
MaxPoints50000采样点的最大数目
"Partitioning"1沿各轴划分积分区域
"SymbolicProcessing"0符号预处理所用的秒数

"QuasiMonteCarlo" 选项.

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

其中 是在 上的均匀分布的函数 作为随机变量的平均值(期望),也就是说, 的概率密度分布为 vol(V)-1Boole(xV). Boole(xV) 代表区域 上的特征函数,而 代表 的体积.

原始蒙特卡洛估计由 "MonteCarloRule" 实现. 积分和误差估计的公式在教程 "NIntegrate 积分规则" "MonteCarloRule" 一节中给出.

考虑积分

如果原始积分区域 被划分为一组不相交的子区域的集合 ,则积分估计为

积分误差为

在每一个子区域上所用的采样点的个数通常可以不同,但在蒙特卡洛算法中所有的 都相等().

分区 被称作层化,每一个 被称作层. 层化可用来改善原始蒙特卡洛估计. (自适应蒙特卡洛算法使用的是递归层化.)

AccuracyGoal 和 PrecisionGoal

当使用 NIntegrate 的蒙特卡洛算法时,AccuracyGoalPrecisionGoal 的默认值分别是 Infinity 和 2.

MaxPoints

选项 MaxPoints 指定的是计算一个积分的蒙特卡洛估计时,所使用的(准)随机采样点的最大数目.

这个例子中,采样点的数目达到了最大值,NIntegrate 终止并给出一个提示信息:

"RandomSeed"

选项 "RandomSeed" 的值用于给采样点随机生成器提供种子以便积分. 关于这一点,蒙特卡洛方法中 "RandomSeed" 的用法与 SeedRandomRandomReal 的用法相似.

通过使用 "RandomSeed",蒙特卡洛积分的结果可以复制. 下面的两种运行结果相同.

使用 "RandomSeed" 的蒙特卡洛积分:
这个蒙特卡洛积分给出同一结果:
蒙特卡洛积分中所用的前20个点:
这些点与使用 SeedRandomRandomReal 得到的点是重合的.

层化原始蒙特卡洛积分

在层化蒙特卡洛积分中,区域被分成几个子区域,并对每个子区域分别应用原始蒙特卡洛估计.

原始蒙特卡洛和准蒙特卡洛策略开端的期望 (平均) 值公式方程 (45) 得到

令区域 被二分为两个半区,. 上的期望值 上的方差. 由定理 [PrFlTeuk92]

可以看到层化采样的方差决不会大于原始蒙特卡洛采样的方差.

指定 "MonteCarlo" 策略的层(strata)有两种方式. 一种是在变量范围指定中指明奇点,另一种是使用该方法的子选项 "Partitioning".

使用变量范围指定的层化原始蒙特卡洛积分:
使用子选项 "Partitioning" 层化原始蒙特卡洛积分:

如果赋以 "Partitioning" 一个整数列 ,长度 等于积分变量的个数,积分区域的每个维度 被平分为 个部分. 如果赋以 "Partitioning" 一个整数 ,则所有的维度上都被平分为 个部分.

这个图展示了由 "Partitioning" 指定的层化采样. 每个单元包含 个点,这是由 "MonteCarloRule" 选项 "Points" 指定的:

可以通过中间奇点指定变量范围,从而对层化蒙特卡洛采样进行指定.

通过中间奇点指定层化蒙特卡洛采样:

层化采样提高了原始蒙特卡洛积分的效率:如果层数为 ,层化蒙特卡洛估计的标准偏差是原始蒙特卡洛估计的标准偏差的1/. (见下面的例子.)

下面的基准显示,层化提高了的收敛速度:

层化蒙特卡洛积分的收敛加速

下面的例子证实,如果层数为 ,层化蒙特卡洛估计的标准偏差是原始蒙特卡洛估计的标准偏差的1/.

这是一个基于期望(平均)值公式 (46) 的层化积分定义:
考虑这个积分:
层数由1 到40 变化,用1000个点对上述积分近似:
这是标准偏差与非层化、原始蒙特卡洛估计的比值:
注意 ratiosi 是层数为 i 的蒙特卡洛估计的比例. 这允许用户尝试函数 ratios 的最小二乘拟合:
的拟合表示出一个非常接近于1的系数,这印证了一种经验法则:层数为 ,收敛速度为原来的 倍. 这是 ratios 最小二乘拟合的图形:

全局自适应蒙特卡洛和准蒙特卡洛策略

全局自适应蒙特卡洛和准蒙特卡洛策略将方差估计最大的子区域进行递归二分,然后计算每一半的积分和方差估计.

这是一个自适应蒙特卡洛积分的例子:
选项名称
默认值
MethodMonteCarloRuleMonteCarloRule 规定
"Partitioning"Automatic在每个轴上对积分区域的的初始分区
"BisectionDithering"0从区域与二分轴平行的边的中间进行抵消
"MaxPoints"Automatic所用 () 随机采样点的最大数目
"RandomSeed"Automatic用于产生 () 随机采样点的随机种子

自适应 (拟) 蒙特卡洛在每个子区域上使用原始 (拟) 蒙特卡洛估计规则.

子区域二分和后继半积分区域的过程将减少全局方差,它被称为递归分层采样. 这是出于这样一则定理,该定理指出,如果一个区域被分成不相交的子区域,则在该区域上的随机变量方差大于或等于在每一个子区域上的随机变量方差之和. (见"原始蒙特卡洛和准蒙特卡洛策略"一节中的"层化蒙特卡洛积分".)

全局自适应蒙特卡洛 "AdaptiveMonteCarlo""GlobalAdaptive" 相似,但也有一些重要的差异.

1.  "AdaptiveMonteCarlo" 不使用奇点平缓,不具有慢速收敛和噪音积分的探测器.

2.  "AdaptiveMonteCarlo" 随机选择二分维度. 为了避免不同坐标的不规则分离,仅当其它维度都已被选过,才能重复一个维度.

3.  可以调节 "AdaptiveMonteCarlo",使对子区域的二分远离中间. 更多内容参见 "BisectionDithering".

MinRecursion 和 MaxRecursion

自适应蒙特卡洛方法的选项 MinRecursionMaxRecursion 与在 "GlobalAdaptive" 中的意义和功能相同. 参见 "MinRecursion 和 MaxRecursion".

"Partitioning"

"AdaptiveMonteCarlo" 的选项 "Partitioning" 提供积分的初始层化. 它的意义和功能与 "MonteCarlo" 策略的 "Partitioning" 相同.

"BisectionDithering"

当被积函数具有某种特殊对称性,使其重要的部分置于区域中间时,稍微远离中间再进行区域的二分比较好. 选项值 "BisectionDithering"->dith 设定该区域的分裂维度边的分裂分数应该是 而不是 . dith 的符号以随机方式进行变化. 赋给 "BisectionDithering" 的默认值是 . dith 的允许值是在区间 上的实数.

考虑这个函数:
考虑这个积分:
如果不进行二分抖动,即 "BisectionDithering" 为 0,则该积分严重估计不足:
下图表明了积分估计不足的原因. 黑色的点是积分采样点. 可以看到被积函数中的一半峰值的采样不足:
指定10%的二分抖动给出了一个满意的估计:
从这幅图中可以看出,对被积函数的峰值的采样改善了:

二分轴的选择

对于多维积分,自适应蒙特卡洛算法选择积分区域的分裂轴有两种方式:(1) 随机选择,或 (2) 使每一半积分估计的方差最小化. "MonteCarloRule"轴选法负责.

实例:与原始蒙特卡洛比较

一般地,"AdaptiveMonteCarlo" 策略的性能好于 "MonteCarlo". 这可以通过这一小节中的例子证明.

考虑这个函数:
这是该函数的图形:

由下面的分析可以看出,"AdaptiveMonteCarlo" 使用的采样点接近于原始 "MonteCarlo" 策略的1/3.

这是 "MonteCarlo" 的采样点和计时:
这是 "AdaptiveMonteCarlo" 的采样点和计时:
这是准确结果:
这是使用 "MonteCarlo" 进行100 次积分所用的的时间:
"MonteCarlo" 积分结果可以与准确结果媲美. 下面的数字显示了积分估计平均值的误差,积分估计的相对误差的平均值,以及积分估计的方差:
这是使用 "AdaptiveMonteCarlo" 进行100 次积分所用的的时间,其速度是 "MonteCarlo" 积分的好几倍:
"AdaptiveMonteCarlo" 积分结果可以与准确结果媲美. 下面的数字显示了积分估计平均值的误差,积分估计的相对误差的平均值,以及积分估计的方差:

"MultiPeriodic"

策略 "MultiPeriodic" 将所有积分变换到一个单位立方体上进行,并将被积函数的周期变换为相对于每一个积分变量为1. 可以对不同的变量应用不同的周期变换函数(或不用). "MultiPeriodic" 用于维数小于或等于12的积分. 如果 "MultiPeriodic" 用于较高维数的积分,则采用 "MonteCarlo" 策略.

这是一个使用 "MultiPeriodic" 进行积分的例子:
选项名称
默认值
"Transformation"SidiSin应用于被积函数的周期化变换
"MinPoints"0采样点的最小数目
"MaxPoints"105采样点的最大数目
"SymbolicProcessing"Automatic符号预处理所用的秒数

可以将 "MultiPeriodic" 看作 "Trapezoidal" 策略在多维上的推广,也可以将其看作一种准蒙特卡洛方法.

"MultiPeriodic" 使用格积分规则,见[SloanJoe94] [KrUeb98].

这里在 上的积分格被理解为 的一个离散子集,其中 对于加减法均封闭,并且包含 . 一个格积分规则 [SloanJoe94] 是一个形如下式的规则:

其中 一个包含在 中的积分格的所有点.

如果 "MultiPeriodic" 被调用,一个 维积分选项 "Transformation" 取单自变量函数的一个列表 ,该列表用于对相应的变量进行变换. 如果赋给 "Transformation" 一个长度为 ,且小于 的列表,则最后一个函数 被用于最后的 个积分变量. 如果赋给 "Transformation" 一个函数,则该函数将被用于对所有变量进行变换.

为积分的方向(维度). 如果 "MultiPeriodic" 在周期性变换后调用 "Trapezoidal". 维数高于12 时,则不进行周期性变换而直接调用"MonteCarlo". 在 "MultiPeriodic" 使用所谓的 复制规则. 对于每一个 "MultiPeriodic" 有一个复制规则的集合,用于计算积分估计的序列. 点数较少的规则首先被使用. 如果一个规则的误差估计达到精确度目标,或者是在该序列中两个积分估计的差值达到精确度目标,积分终止.

不同维数的规则集合中用于 复制规则的点数:

与 "MultidimensionalRule" 进行比较

一般地,"MultiPeriodic" 慢于使用 "MultidimensionalRule""GlobalAdaptive". 对于低精确度、高维数积分的计算,"MultiPeriodic" 可能较快地给出结果.

定义一个具有8个自变量的函数:
使用 "MultidimensionalRule""MultiPeriodic""GlobalAdaptive" 计算 的时间(以秒计):
使用 "MultidimensionalRule""MultiPeriodic""GlobalAdaptive" 计算 时,被积函数进行计算的次数:

预处理器

所有策略的功能都可以通过积分的符号预处理进行扩展. 预处理器可以被视为将积分任命给其它策略(包括预处理器)的策略.

"SymbolicPiecewiseSubdivision"

"SymbolicPiecewiseSubdivision" 是这样一种预处理器,它将被积函数为分段形式的积分分成不相交积分区域的积分,在每一个不相交积分区域上,被积函数不再是分段函数.

选项名称
默认值
MethodAutomatic积分被传递至的积分策略或预处理器
"ExpandSpecialPiecewise"True应该扩展哪一个分段函数
TimeConstraint5试图进行分段划分所用的最大秒数
"MaxPiecewiseCases"100分段预处理器能够返回的子区域的最大个数
"SymbolicProcessing"Automatic符号预处理所用的秒数

"SymbolicPiecewiseSubdivision" 选项.

正如曾在该教程的一开始所提到的,NIntegrate 能够同时进行不相交的区域上不同被积函数的积分. 因此,在通过"SymbolicPiecewiseSubdivision" 进行预处理之后,积分按相同的方式继续进行,正如给出了NIntegrate 的奇点指定的范围(可以看作不相交区域上的积分指定相同的被积函数). 例如,策略 "GlobalAdaptive" 尽量改善通过一种区域的积分估计,这种区域经二分后误差最大,并且不管它对应于哪个被积函数,都将选择最大的误差区域.

下面是积分 的数值估计的采样点. 在图中,按照 ord 坐标的顺序在 上对被积函数进行采样. 对于函数段 ,可以看到 "GlobalAdaptive" 将函数段 的采样变更为对函数段 的采样:
这是积分 的数值估计的采样点. 在左边画出了被积函数,在右边画出了采样点. 该积分被划分为 + + + ,这是为什么采样点在 上构成了不同的布局:

"ExpandSpecialPiecewise"

在某些情况下,最好是只对某些分段函数进行分段扩展. 这时选项 "ExpandSpecialPiecewise" 可以被赋以一个函数列表以进行分段扩展.

这个蒙特卡洛积分只有在 Boole 时进行分段扩展,速度较快:
既在 Boole 也在 Abs 上进行分段扩展的蒙特卡洛积分:

"EvenOddSubdivision"

"EvenOddSubdivision" 是这样一种预处理器,如果区域关于原点对称,并且被积函数被确定为奇函数或偶函数的一种,则进行积分区域的简化. 奇函数积分的收敛性被默认验证.

选项名称
默认值
MethodAutomatic积分被传递至的积分策略或预处理器
VerifyConvergenceAutomatic如果检测到奇函数,是否验证其收敛性
"SymbolicProcessing"Automatic符号预处理所用的秒数

"EvenOddSubdivision" 的选项.

当被积函数是一个偶函数,且积分区域关于原点对称时,可以仅在某些积分区域上进行积分然后乘以一个相应的因子来计算积分.

这是一个没有经过任何预处理的偶函数及其采样点的图形:
在应用 "EvenOddSubdivision" 之后被 NIntegrate 所使用的采样点. 注意采样点仅在区域 上:

变换定理

预处理器 "EvenOddSubdivision" 基于下述定理.

定理:已知 维积分

假设对于某个 下列不等式成立:

a) ,

b) 对于所有 :

换言之, 的范围关于原点对称,并且变量 的边界关于 为偶函数.

则:

a) 如果被积函数关于 为偶函数,亦即

积分等价于

b) 如果被积函数关于 为积函数,亦即

则积分等于0.

上面的定理可以在一个积分上多次应用.

为了解释该定理,考察积分 . 它关于 轴对称,并且被积函数与 的边界是关于 轴的偶函数.

这是采样点分别在不应用 "EvenOddSubdivision" (黑色)和应用 "EvenOddSubdivision"(红色)的图形:

如果 的边界关于 不是偶函数,就破坏了沿 的对称性. 例如,积分 没有 NIntegrate 可以利用的对称性.

这是采样点应用 "EvenOddSubdivision"(红色)时的图形. 该区域沿 轴无对称性:

"VerifyConvergence"

考察下面的发散积分 . NIntegrate 检测到其在一个对称区域上是奇函数,并试图对 进行积分(也就是说,检查 的收敛性). 由于提示信息 ncvb 表明没有达到收敛,NIntegrate 给出积分可能发散的提示信息 oidiv
如果选项 VerifyConvergence 被设置为 False,则当发现被积函数是奇函数后,不进行收敛性验证,也不进行被积函数的计算:

"OscillatorySelection"

"OscillatorySelection" 是一个预处理器,它选择有效计算振荡积分的专用算法.

在每个积分区域上,"OscillatorySelection" 检测被积函数是否为傅立叶有限范围类型、傅立叶无限范围类型(见 "DoubleExponentialOscillatory")、贝塞尔无限范围类型(见"ExtrapolatingOscillatory")、 Levin 类型(见 "LevinRule"),还是没有特殊振荡类型. 选项控制的是对各类被积函数所选择的方法.

选项名称
默认值
"BesselInfiniteRangeMethod"Automatic用于无限范围贝塞尔积分的方法
"FourierFiniteRangeMethod"Automatic用于有限范围傅立叶积分的方法
"FourierInfiniteRangeMethod"Automatic用于无限范围傅立叶积分的方法
"LevinMethod"Automatic用于 Levin 类型振荡积分的方法
Method"GlobalAdaptive"用于非振荡积分的方法
"TermwiseOscillatory"False是否单独处理和式中的项
"SymbolicProcessing"5符号处理所用的秒数

"OscillatorySelection" 方法选项.

"OscillatorySelection" 预处理器默认用于 NIntegrate 中:
如果没有 "OscillatorySelection" 预处理器,NIntegrate 在默认选项设置时不会达到收敛:

预处理器 "OscillatorySelection" 是被设计用来与预处理器 "SymbolicPiecewiseSubdivision" 的内部输出合作. "OscillatorySelection" 自身对下述振荡积分进行分区:振荡积分包含原点,或者是含有需要扩展的振荡内核,或者是振荡内核需要变换成针对振荡算法所设计的形式

这是一个分段函数的积分,其中一维积分的所有专用积分策略由 "OscillatorySelection" 自动选择. 对于这个积分,预处理器 "SymbolicPiecewiseSubdivision" 首先将积分分成四个不同的积分;针对每一个积分,"OscillatorySelection" 分别选择一个适当的专用算法:
下表表明了用于指定上述积分的每一个子区间算法的 "OscillatorySelection" 选项的名称:
x(-,0]"BesselInfiniteRangeMethod"
x[0,20]"FourierFiniteRangeMethod"
x[30,)"FourierInfiniteRangeMethod"
x[20,30]Method
在这个例子中,"DoubleExponentialOscillatory" 被调用了两次. "DoubleExponentialOscillatory" 是针对傅立叶积分的一种专用算法,公式 将被积函数变成两个傅立叶型被积函数的和:
为了表明 "OscillatorySelection" 使用了公式 ,这里手动将上述积分分开,结果与上面的结果完全相同:

选项 "FourierFiniteRangeMethod" 的值为 Automatic 意味着,如果由选项 Method 指定的积分策略是 "GlobalAdaptive""LocalAdaptive" 的一种,则那个策略将被用于有限范围上的傅立叶积分,否则使用 "GlobalAdaptive".

这是一个分段函数的积分,对于非振荡积分使用 "DoubleExponential" 策略,对于有限范围的振荡积分使用 "LocalAdaptive" 策略:
这是默认选项设置下,前面积分的采样点和积分. 左图中,在 之间的模式是局部自适应求积法的典型形式可以看到积分被递归分成了三个部分(由于选项 "Partitioning"->3 被赋给了 "LocalAdaptive"). 右图中,在 之间的模式T来自于 "GlobalAdaptive". 第一副图形中,在 之间的模式是双指数求积法的典型形式. 在第二副图形中,在 之间可以看到相同的模式,这是由于"GlobalAdaptive" 默认时使用 "DoubleExponential" 奇点处理器:

对于一种振荡积分的特殊类型,如果应用一个特定的振荡方法是理想的,则要么是 "OscillatorySelection" 相应的选项应该改变,要么是在不使用预处理器 "OscillatorySelection" 的情况下使用 NIntegrate 中的 Method 选项.

这是一个分段函数的积分,对于任何一个无穷范围振荡积分使用 "ExtrapolatingOscillatory"
如果 "ExtrapolatingOscillatory" 被作为方法给出,"OscillatorySelection" 则使用它进行无穷范围振荡积分:
使用默认选项 NIntegrate,可以较快地进行上述积分. 对于这个积分,默认使用的 "OscillatorySelection" 使用的是 "DoubleExponentialOscillatory"

"UnitCubeRescaling"

"UnitCubeRescaling" 是一个将积分区域变换到一个单位立方体或单位超立方体上的预处理器. 原始被积函数的变量被替代,结果被乘以变换的雅可比量.

选项名称
默认值
"FunctionalRangesOnly"True哪个范围应该被变换到单位立方体上
Method"GlobalAdaptive"积分被传递至的积分策略或预处理器
"SymbolicProcessing"Automatic符号预处理所用的秒数

"UnitCubeRescaling" 选项.

这里使用单位立方体调整尺度,其计算速度较快:
这个积分不使用单位立方体调整尺度,它的计算速度是上面一个的将近1/3:

"UnitCubeRescaling" 将积分

变换到单位超立方体 上. 假设 为有限值, 为分段连续函数,则由"UnitCubeRescaling" 进行的变换为

该变换的雅可比量是

如果对于第 个轴, 或两个都是有限的,则 (47) 中 的公式是一个将 映射到 的非仿射变换. NIntegrate 使用的是下述变换:

其中 .

如果积分区域的边界是有限或无限常数,"UnitCubeRescaling" 的应用会使被积函数更加复杂. 由于 NIntegrate 有效率的仿射变换和无限内部变量变换,积分过程会变得较慢. 如果积分区域的边界中有一些是函数,应用 "UnitCubeRescaling" 将使积分变快,这是由于牵涉到积分变量的计算在仅当被积函数被计算时才进行. 出于这些性能的考虑,"UnitCubeRescaling" 中存在选项 "FunctionalRangesOnly". 如果 "FunctionalRangesOnly" 被设置为 True,尺寸的重新调整仅应用于多维函数的范围.

这里积分使用单位立方体调整尺度:
这个积分不使用单位立方体调整尺度,它的计算速度是上面一个的将近2倍:

示例

"UnitCubeRescaling" 所使用的变换过程与下面由 FRangesToCube 所执行的操作相同(在"Duffy 的坐标推广与示例"中也有定义).

这个函数提供了一个积分范围的列表和一列平行六面体各面或一个超立方体面的列表的变换 (48) 及其雅可比量 (49):

变换 (50) 的每一个变换可以通过 Rescale 完成.

注意对于给定的轴 ,已经推导得到的针对 的变换规则,要在沿第 个轴进行边界的尺寸调整之前被应用到原始的边界上.

用于 的变换规则及其雅可比量:
将变换应用到一个函数:
用于 的变换规则及其雅可比量:
用于 的变换规则及其雅可比量:

"SymbolicPreprocessing"

"SymbolicPreprocessing" 是一个复合预处理器,用于简化其它预处理器的打开和关闭.

选项名称
默认值
MethodAutomatic积分被传递至的积分策略或预处理器
"SymbolicPiecewiseSubdivision"True分段再分
"EvenOddSubdivision"True偶-奇再分
"OscillatorySelection"True振荡函数乘积的检测
"UnitCubeRescaling"Automatic调整尺寸至单位超立方体上
"SymbolicProcessing"Automatic符号预处理所用的秒数

"SymbolicPreprocessing" 选项.

"UnitCubeRescaling" 被设置为 Automatic 时,它只应用于多维函数范围.

这是应用不同的预处理器组合对 进行积分的例子:

实例和应用

封闭式围道积分

这个函数计算的是一个极坐标参数化的封闭式围道的质量:
这是使用 Integrate 得到的半径为2和3的椭圆的周长:
这是使用同一函数得到的半径为2和3的椭圆的周长近似值:
这个近似值可以与准确值相媲美:

傅立叶级数计算

这是一个 Wolfram 语言中的函数,计算一个函数的截断傅立叶级数近似:
使用 Integrate 得到的 上的傅立叶近似:
这是 的图形和傅立叶级数近似:
这里使用 NIntegrate 计算得到了 Sin[x3+] 上的60项的傅立叶近似. 如果使用的是 Integrate,计算过程会非常慢:
这是 Sin[x3+] 和傅立叶级数近似的图形: