Python SciPy integrate.nquad用法及代码示例

本文简要介绍 python 语言中 scipy.integrate.nquad 的用法。

用法:

scipy.integrate.nquad(func, ranges, args=None, opts=None, full_output=False)#

集成多个变量。

包装 quad 以启用对多个变量的集成。各种选项允许改进不连续函数的积分，以及加权积分的使用，以及对积分过程的一般更精细的控制。

参数 ：：

func： {callable, scipy.LowLevelCallable}

要集成的函数。具有 x0, ... xn 、 t0, ... tm 的参数，其中积分是在 x0, ... xn 上进行的，它必须是浮点数。其中 t0, ... tm 是在 args 中传递的额外参数。函数签名应该是 func(x0, x1, ..., xn, t0, t1, ..., tm) 。整合有序进行。也就是说，x0 上的积分是最内层的积分，而xn 是最外层的积分。

如果用户希望提高集成性能，那么f可能是一个scipy.LowLevelCallable带有其中一个签名：

double func(int n, double *xx)
double func(int n, double *xx, void *user_data)

其中n 是变量和参数的数量。 xx 数组包含坐标和额外参数。 user_data 是 scipy.LowLevelCallable 中包含的数据。

ranges： 可迭代对象

范围的每个元素可以是 2 个数字的序列，也可以是返回此类序列的可调用对象。 ranges[0] 对应于 x0 上的积分，依此类推。如果范围的元素是可调用的，那么将使用所有可用的集成参数以及任何参数参数来调用它。例如，如果 func = f(x0, x1, x2, t0, t1) ，则 ranges[0] 可以定义为 (a, b) 或 (a, b) = range0(x1, x2, t0, t1) 。

args： 可迭代对象，可选

func 、 ranges 和 opts 所需的附加参数 t0, ... tn 。

opts： 可迭代对象或字典，可选

要传递给 quad 的选项。可能是空的、一个字典、或者一系列字典或返回字典的函数。如果为空，则使用 scipy.integrate.quad 中的默认选项。如果是字典，则相同的选项用于所有级别的集成。如果是序列，则序列的每个元素对应于特定的积分。例如， opts[0] 对应于 x0 上的积分，依此类推。如果是可调用的，则签名必须与 ranges 相同。可用选项及其默认值是：

• epsabs = 1.49e-08

• epsrel = 1.49e-08

• limit = 50

• points = None

• weight = None

• wvar = None

• wopts = None

有关这些选项的更多信息，请参阅 quad 。

full_output： 布尔型，可选

来自 scipy.integrate.quad 的 full_output 的部分实现。被积函数求值次数neval可以通过调用nquad时设置full_output=True来获得。

返回 ：：

result： 浮点数

整合的结果。

abserr： 浮点数

各种积分结果中绝对误差估计值的最大值。

out_dict： 字典，可选

包含有关集成的附加信息的字典。

注意：

为了获得有效的结果，积分必须收敛；不保证发散积分的行为。

QUADPACK关卡例程的详细信息

nquad调用 FORTRAN 库 QUADPACK 中的例程。本节提供有关调用每个例程的条件的详细信息以及每个例程的简短说明。调用的例程取决于重量,积分和积分限制a和b.

四包例程	重量	积分	无限界限
qagse	None	No	No
qagie	None	No	Yes
qagpe	None	Yes	No
qawoe	‘sin’, ‘cos’	No	No
qawfe	‘sin’, ‘cos’	No	A或B
qawse	‘阿尔格*’	No	No
qawce	‘cauchy’	No	No

下面提供了 [1] 中每个例程的简短说明。

察格色：

是一个基于全局自适应区间细分与外推法相结合的积分器，它将消除几种类型的被积函数奇异性的影响。

卡吉：

处理无限间隔内的积分。无限范围映射到有限区间，随后应用与 QAGS 中相同的策略。

察格佩：

与 QAGS 的用途相同，但还允许用户提供有关 trouble-spots 的位置和类型的明确信息，即内部奇点、不连续性和被积函数的其他困难的横坐标。

卡沃伊：

是用于在有限区间 [a,b] 上评估 \(\int^b_a \cos(\omega x)f(x)dx\) 或 \(\int^b_a \sin(\omega x)f(x)dx\) 的积分器，其中 \(\omega\) 和 \(f\) 由用户指定。规则评估组件基于修改后的Clenshaw-Curtis技术

自适应细分方案与外推程序结合使用，该方案是对 QAGS 中的方案的修改，并允许算法处理 \(f(x)\) 中的奇点。

夸夫：

计算用户提供的 \(\omega\) 和 \(f\) 的傅里叶变换 \(\int^\infty_a \cos(\omega x)f(x)dx\) 或 \(\int^\infty_a \sin(\omega x)f(x)dx\) 。 QAWO 的过程应用于连续的有限间隔，并且通过\(\varepsilon\) 算法的收敛加速应用于一系列积分近似。

夸瑟：

近似 \(\int^b_a w(x)f(x)dx\) 和 \(a < b\) 其中 \(w(x) = (x-a)^{\alpha}(b-x)^{\beta}v(x)\) 和 \(\alpha,\beta > -1\) ，其中 \(v(x)\) 可以是以下函数之一： \(1\) 、 \(\log(x-a)\) 、 \(\log(b-x)\) 、 \(\log(x-a)\log(b-x)\) 。

用户指定\(\alpha\),\(\beta\)和函数的类型\(v\)。应用全局自适应细分策略，并对包含以下内容的子区间进行修改后的 Clenshaw-Curtis 积分：a或者b.

夸斯：

计算 \(\int^b_a f(x) / (x-c)dx\)，其中对于用户指定的 \(c\) 和 \(f\) ，积分必须解释为柯西主值积分。该策略具有全局适应性。修改后的 Clenshaw-Curtis 积分用于包含点 \(x = c\) 的那些区间。

参考：

[1]

罗伯特·皮森斯； de Doncker-Kapenga，伊莉斯；克里斯托夫·W·尤伯胡贝尔；大卫·卡哈纳 (1983)。 QUADPACK：用于自动集成的子程序包。 Springer-Verlag。 ISBN 978-3-540-12553-2。

例子：

计算

\[\int^{1}_{-0.15} \int^{0.8}_{0.13} \int^{1}_{-1} \int^{1}_{0} f(x_0, x_1, x_2, x_3) \,dx_0 \,dx_1 \,dx_2 \,dx_3 ,\]

其中

\[\begin{split}f(x_0, x_1, x_2, x_3) = \begin{cases} x_0^2+x_1 x_2-x_3^3+ \sin{x_0}+1 & (x_0-0.2 x_3-0.5-0.25 x_1 > 0) \\ x_0^2+x_1 x_2-x_3^3+ \sin{x_0}+0 & (x_0-0.2 x_3-0.5-0.25 x_1 \leq 0) \end{cases} .\end{split}\]

>>> import numpy as np
>>> from scipy import integrate
>>> func = lambda x0,x1,x2,x3 : x0**2 + x1*x2 - x3**3 + np.sin(x0) + (
...                                 1 if (x0-.2*x3-.5-.25*x1>0) else 0)
>>> def opts0(*args, **kwargs):
...     return {'points':[0.2*args[2] + 0.5 + 0.25*args[0]]}
>>> integrate.nquad(func, [[0,1], [-1,1], [.13,.8], [-.15,1]],
...                 opts=[opts0,{},{},{}], full_output=True)
(1.5267454070738633, 2.9437360001402324e-14, {'neval': 388962})

计算

\[\int^{t_0+t_1+1}_{t_0+t_1-1} \int^{x_2+t_0^2 t_1^3+1}_{x_2+t_0^2 t_1^3-1} \int^{t_0 x_1+t_1 x_2+1}_{t_0 x_1+t_1 x_2-1} f(x_0,x_1, x_2,t_0,t_1) \,dx_0 \,dx_1 \,dx_2,\]

其中

\[\begin{split}f(x_0, x_1, x_2, t_0, t_1) = \begin{cases} x_0 x_2^2 + \sin{x_1}+2 & (x_0+t_1 x_1-t_0 > 0) \\ x_0 x_2^2 +\sin{x_1}+1 & (x_0+t_1 x_1-t_0 \leq 0) \end{cases}\end{split}\]

和\((t_0, t_1) = (0, 1)\)。

>>> def func2(x0, x1, x2, t0, t1):
...     return x0*x2**2 + np.sin(x1) + 1 + (1 if x0+t1*x1-t0>0 else 0)
>>> def lim0(x1, x2, t0, t1):
...     return [t0*x1 + t1*x2 - 1, t0*x1 + t1*x2 + 1]
>>> def lim1(x2, t0, t1):
...     return [x2 + t0**2*t1**3 - 1, x2 + t0**2*t1**3 + 1]
>>> def lim2(t0, t1):
...     return [t0 + t1 - 1, t0 + t1 + 1]
>>> def opts0(x1, x2, t0, t1):
...     return {'points' : [t0 - t1*x1]}
>>> def opts1(x2, t0, t1):
...     return {}
>>> def opts2(t0, t1):
...     return {}
>>> integrate.nquad(func2, [lim0, lim1, lim2], args=(0,1),
...                 opts=[opts0, opts1, opts2])
(36.099919226771625, 1.8546948553373528e-07)