Wolfram 语言中的数值积分简介
概述
Wolfram 语言函数 NIntegrate 处理的是一般的数值积分. 它可以处理大量的一维和多维积分.
| NIntegrate[f[x1,x2,…,xn],{x1,a1,b1},{x2,a2,b2},…,{xn,an,bn}] | |
| 在区域 [a1,b1]×[a2,b2]×…×[an,bn] 上找到函数 f 的数值积分 | |
一般来说，NIntegrate 通过在积分区域上对被积函数值进行采样进行积分的估计. 各种数值积分方法对初步采样步骤，以及采样如何演变做出规定.
NIntegrate[(1/Sqrt[x]), {x, 0, 1}]
ListPlot[Last[Reap[NIntegrate[(1/Sqrt[x]), {x, 0, 1}, EvaluationMonitor :> Sow[x]]]]]NIntegrate 使用的算法被称作"积分策略"，尝试计算积分估计，使之满足用户指定的精度和准确度的目标. 积分策略运用 "积分规则" 使用加权和计算被积函数值的积分估计.
NIntegrate[(1/2 + Sin[x]), {x, 0, 5}, Method -> {"GlobalAdaptive", Method -> "ClenshawCurtisRule"}]符号预处理
NIntegrate 使用符号预处理简化特殊结构的积分且自动选择积分方法. 被积函数是偶函数或奇函数或包含分段函数可能导致积分区域被转换成或分成多个不同的积分区域.
NIntegrate[(1/Sqrt[Abs[x - 1]]), {x, 0, 2}]
NIntegrate[(x - 2)^2 Sin[4000 x], {x, 2, 3}]
Plot[(x - 2)^2 Sin[4000 x], {x, 2 + (2/50), 2 + (3/50)}]符号预处理允许自动计算包括不连续性和快速变化区域的各种积分.
NIntegrate[Piecewise[{{(1/Sqrt[Abs[x - 1]]), x < 2}, {(x - 2)^2 Sin[4000 x], 2 < x < 3}}], {x, 0, 3}]求积规则
NIntegrate 包括大部分经典的一维求积规则.
Table[NIntegrate[(1/1 + Cos[x + x^2]), {x, 0, 1}, Method -> rule], {rule, {"NewtonCotesRule", "GaussKronrodRule", "ClenshawCurtisRule"}}]经典的求积规则是对被取样的被积函数值形成线性组合. 对于每个求积规则，线性组合中的加权向量是固定的.
{nodes, weights} = NIntegrate`NewtonCotesRuleData[2, MachinePrecision][[{1, 2}]]对于多维积分，NIntegrate 包含一类基于稀疏网格的规则，并且允许从一维规则的笛卡尔积形成规则.
NIntegrate[(1/1 + x + Cos[y]), {x, 0, 1}, {y, 0, 2}, Method -> "MultidimensionalRule"]NIntegrate[(1/1 + x + Cos[y]), {x, 0, 1}, {y, 0, 2}, Method -> {"GaussKronrodRule", "LobattoKronrodRule"}]Table[ListPlot[Last[Reap[NIntegrate[(1/1 + x + Cos[y]), {x, 0, 1}, {y, 0, 2}, Method -> rule, EvaluationMonitor :> Sow[{x, y}], MaxRecursion -> 0]]], PlotLabel -> rule, ImageSize -> 250], {rule, {"MultidimensionalRule", {"GaussKronrodRule", "LobattoKronrodRule"}}}]//QuietBoole 可用于指定更多复杂的多维区域. 这种方式指定的区域在符号预处理时可进一步简化.
Plot3D[Boole[x^2 + y^2 < 1] * (1 - Sqrt[x^2 + y^2]), {x, -1, 1}, {y, -1, 1}, PlotRange -> All]NIntegrate[Boole[x^2 + y^2 < 1] * (1 - Sqrt[x^2 + y^2]), {x, -1, 1}, {y, -1, 1}]
points = First[Last[Reap[NIntegrate[Boole[x^2 + y^2 < 1] * (1 - Sqrt[x^2 + y^2]), {x, -1, 1}, {y, -1, 1}, EvaluationMonitor :> Sow[{x, y}]]]]];Graphics[{PointSize[0.01], Point[points]}, Axes -> True, PlotRange -> {{-1, 1}, {-1, 1}}]
的边界进行密集采样. 由于没有应用符号预处理，NIntegrate 必须高效地数值定位此边界：points = First[Last[Reap[NIntegrate[Boole[x^2 + y^2 < 1] * (1 - Sqrt[x^2 + y^2]), {x, -1, 1}, {y, -1, 1}, Method -> {Automatic, "SymbolicProcessing" -> None}, EvaluationMonitor :> Sow[{x, y}]]]]];
Graphics[{PointSize[0.005], Point[points]}, Axes -> True, AxesOrigin -> {-1, -1}]振荡积分
NIntegrate 包含适用于各种被积函数的普通的振荡积分方法，可用于一维或多维的有限或无限区域. 另外， NIntegrate 包含许多方法，特别适用于含有 Exp，三角函数诸如 Sin 和 Cos，以及某种其它特殊函数诸如 BesselJ 等特殊形式的一维积分函数.
NIntegrate[(Sin[x + Sqrt[x]] + x BesselJ[0, x^2]/1 + x), {x, 0, ∞}]Plot[(Sin[x + Sqrt[x]] + x BesselJ[0, x^2]/1 + x), {x, 0, 12}]自动奇点处理
NIntegrate 处理奇异被积函数有多种方式. 确定性自适应策略 "GlobalAdaptive" 和 "LocalAdaptive" 采用奇点处理技术（基于变量变换）来加速积分过程的收敛. "DoubleExponential" 策略运用梯形求积法，并对被积函数进行一种特殊的变量变换. 在包含积分区间的复平面上的开放集合中，这种规则变换组合实现了被积函数的最佳收敛.
NIntegrate[(1/Sqrt[x]), {x, 0, 1}, PrecisionGoal -> 10]Do[NIntegrate[(1/Sqrt[x]), {x, 0, 1}, PrecisionGoal -> 10], {100}]//TimingDo[NIntegrate[(1/Sqrt[x]), {x, 0, 1}, Method -> {"GlobalAdaptive", "SingularityHandler" -> None}, MaxRecursion -> 100, PrecisionGoal -> 10], {100}]//Timing特殊策略
"DuffyCoordinates" 策略简化或消去了多维积分中的某些类型的奇点. 某种球对称的积分收敛很快.
NIntegrate[(1/Sqrt[x^2 + y^2 + z^2]), {x, 0, 1}, {y, 0, 1}, {z, 0, 1}, Method -> "DuffyCoordinates"]//TimingNIntegrate[(1/Sqrt[x^2 + y^2 + z^2]), {x, 0, 1}, {y, 0, 1}, {z, 0, 1}]//Timing当解析周期被积函数的积分区间正好是一个周期时，"Trapezoidal" 策略给出最佳收敛.
exact = Integrate[Cos[20x]^4, {x, 0, 2(π/20)}]resTrap = NIntegrate[Cos[20x]^4, {x, 0, 2(π/20)}, PrecisionGoal -> 150, WorkingPrecision -> 200, Method -> "Trapezoidal"];//TimingAbs[exact - resTrap]resDef = NIntegrate[Cos[20x]^4, {x, 0, 2(π/20)}, PrecisionGoal -> 150, WorkingPrecision -> 200];//TimingAbs[exact - resDef]对于高维积分，或者只需要一个粗略的积分估计的情况，蒙特卡洛方法是有用的. NIntegrate 有原始蒙特卡洛、自适应蒙特卡洛以及准蒙特卡洛策略可选.
X = Array[x, 30];
NIntegrate[(1/Total@X), Evaluate[Sequence@@Map[{#, 0, 1}&, X]], Method -> "AdaptiveMonteCarlo", PrecisionGoal -> 3]设计
NIntegrate 框架的主要特征是：
- 面向对象（方法属性规范和交流）
- 方法初始化阶段和运行时间计算的分离
- 分层和折返数值方法
- 类型及精密动态方法
- 通过插件功能的用户可扩展性和原型设计
- 专业化数据结构
策略、规则和预处理器
NIntegrate 积分策略是根据它们如何在积分区域采样，如何应用被积函数的类，它们是否是“基于规则”的策略进行分类.
| "GlobalAdaptive" | 任何被积函数，自适应采样，基于规则 |
| "LocalAdaptive" | 任何被积函数，自适应采样，基于规则 |
| "DoubleExponential" | 任何被积函数，均匀采样 |
| "Trapezoidal" | 任何被积函数，均匀采样 |
| "MultiPeriodic" | 多维被积函数，均匀采样 |
| "MonteCarlo" | 任何被积函数，均匀随机采样 |
| "QuasiMonteCarlo" | 任何被积函数，均匀准随机采样 |
| "AdaptiveMonteCarlo" | 任何被积函数，自适应随机采样 |
| "AdaptiveQuasiMonteCarlo" | 任何被积函数，自适应准随机采样 |
| "DoubleExponentialOscillatory" | 一维无限范围振荡被积函数 |
| "ExtrapolatingOscillatory" | 一维无限范围振荡被积函数 |
NIntegrate 积分策略
自适应采样策略是通过在子区域使用更大的错误估计进行更频繁的采样来尝试改善积分估计，一般是细分那些子区域. 均匀采样策略是通过在整个积分区域均匀增加采样密度来尝试改善积分估计.
基于规则策略将一个给定的积分规则应用于每个子区域，以获得该区域的积分和错误估计. 积分规则可以使用设置 Method->{"strategy",Method->"rule"} 来指定.
NIntegrate[(1/Sqrt[x])Log[(1/x)], {x, 0, 1}, Method -> {"GlobalAdaptive", Method -> "ClenshawCurtisRule"}]NIntegrate 积分规则可以根据它们是应用于一维或多维区域，以及积分规则的类型进行分类.
| "BooleRule" | 一维，加权和 |
| "ClenshawCurtisRule" | 一维，加权和 |
| "GaussBerntsenEspelidRule" | 一维，加权和 |
| "GaussKronrodRule" | 一维，加权和 |
| "LobattoKronrodRule" | 一维，加权和 |
| "LobattoPeanoRule" | 一维，加权和 |
| "NewtonCotesRule" | 一维，加权和 |
| "TrapezoidalRule" | 一维，加权和 |
| "ClenshawCurtisOscillatoryRule" | 一维，专用振荡规则 |
| "LevinRule" | 一维或多维，一般振荡规则 |
| "MultipanelRule" | 一维，加权和，一维规则的组合 |
| "CartesianRule" | 多维，加权和，一维规则的乘积 |
| "MultidimensionalRule" | 多维，加权和 |
可以与基于规则策略 "GlobalAdaptive" 和"LocalAdaptive" 一起使用的积分规则.
经典的“加权和”类型规则估计点集的函数值的预设线性组合的积分. 振荡规则使用求积加权估计积分，其依赖于被积函数的特殊振荡“内核”.
组合规则从一个或多个子规则构建求积规则. 它们由设置 Method->{"rule",Method->{"subrule1",…}} 指定.
NIntegrate[Log[1 + x ^ 2 + y], {x, 0, 10}, {y, 0, 10}, Method -> {"CartesianRule", Method -> {"GaussKronrodRule", "ClenshawCurtisRule"}}]所有策略的能力通过被积函数的符号预处理得到延伸. 预处理是由预处理策略控制的，首先变换或分析积分，然后把积分托付给另一个策略（经常是另一个预处理策略）.
| "SymbolicPreprocessing" | 全局预处理控制器 |
| "EvenOddSubdivision" | 简化奇偶被积函数 |
| "InterpolationPointsSubdivision" | 细分包含内插函数的被积函数 |
| "OscillatorySelection" | 检测振荡被积函数并选择合适的方法 |
| "SymbolicPiecewiseSubdivision" | 细分包含分段函数的被积函数 |
| "UnitCubeRescaling" | 重新调整多维被积函数到单位立方体 |
| "DuffyCoordinates" | 多维奇点消除变换 |
| "PrincipalValue" | 等价于柯西主值的数值积分 |
NIntegrate 预处理器策略
预处理策略是由设置 Method->{"preprocessor",Method->m} 来指定的，其中 m 是预处理完成后积分被托付的策略或规则.
NIntegrate[Boole[x^2 + y^2 < 1] * (1 - Sqrt[x^2 + y^2]), {x, -1, 1}, {y, -1, 1}, Method -> {"EvenOddSubdivision", Method -> "LocalAdaptive"}]Graphics[{PointSize[0.005], Point[Reap[NIntegrate[Boole[x^2 + y^2 < 1] * (1 - Sqrt[x^2 + y^2]), {x, -1, 1}, {y, -1, 1}, Method -> {"EvenOddSubdivision", Method -> "LocalAdaptive"}, EvaluationMonitor :> Sow[{x, y}]]][[2, 1]]]}, Axes -> True, PlotRange -> {{-1, 1}, {-1, 1}}]