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

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

用法:

scipy.integrate.quad(func, a, b, args=(), full_output=0, epsabs=1.49e-08, epsrel=1.49e-08, limit=50, points=None, weight=None, wvar=None, wopts=None, maxp1=50, limlst=50, complex_func=False)#

计算一个定积分。

使用 Fortran 库 QUADPACK 中的技术将 func 从 a 集成到 b(可能是无限区间)。

参数 ：：

func： {函数，scipy.LowLevelCallable}

要集成的 Python 函数或方法。如果 func 接受许多参数，它会沿对应于第一个参数的轴进行积分。

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

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

user_data 是包含在 scipy.LowLevelCallable 中的数据。在带有 xx 的调用表单中，n 是包含 xx[0] == x 的 xx 数组的长度，其余项是包含在 quad 的 args 参数中的数字。

此外，为了向后兼容，支持某些 ctypes 调用签名，但不应在新代码中使用这些调用签名。

a： 浮点数

积分的下限(使用 -numpy.inf 表示 -infinity)。

b： 浮点数

积分上限(使用 numpy.inf 表示 +infinity)。

args： 元组，可选

传递给 func 的额外参数。

full_output： 整数，可选

非零返回集成信息字典。如果非零，则警告消息也会被抑制，并且消息会附加到输出元组。

complex_func： 布尔型，可选

指示函数的 (函数) 返回类型为实数 (complex_func=False：默认)或复杂(complex_func=True)。在这两种情况下，函数的参数都是实数。如果 full_output 也不为零，则信息词典,信息， 和解释实数和复数分量在带有键 “real output” 和 “imag output” 的字典中返回。

返回 ：：

y： 浮点数

func 从 a 到 b 的积分。

abserr： 浮点数

结果中绝对误差的估计。

infodict： dict

包含附加信息的字典。

信息

收敛消息。

解释

仅附加 ‘cos’ 或 ‘sin’ 权重和无限积分限制，它包含 infodict[‘ierlst’] 中代码的解释

其他参数 ：：

epsabs： float 或 int，可选

绝对容错。默认值为 1.49e-8。quad试图获得准确度abs(i-result) <= max(epsabs, epsrel*abs(i))其中i= 积分函数从a到b， 和result是数值近似。看埃普雷尔以下。

epsrel： float 或 int，可选

相对误差容限。默认值为 1.49e-8。如果epsabs <= 0,埃普雷尔必须大于 5e-29 和50 * (machine epsilon).看易胜宝更多。

limit： float 或 int，可选

自适应算法中使用的子区间数的上限。

points： (浮点数，整数序列)，可选

有界积分区间中的一系列断点，其中可能出现被积函数的局部困难(例如，奇异点、不连续点)。序列不必排序。请注意，此选项不能与 weight 一起使用。

weight： float 或 int，可选

表示权重函数的字符串。对此和其余论点的完整解释可以在下面找到。

wvar： 可选的

与加权函数一起使用的变量。

wopts： 可选的

用于重用切比雪夫矩的可选输入。

maxp1： float 或 int，可选

切比雪夫矩数的上限。

limlst： 整数，可选

用于正弦加权和无限 end-point 的周期数上限 (>=3)。

注意：

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

quad() 输入和输出的额外信息

如果 full_output 不为零，则第三个输出参数 (infodict) 是一个字典，其条目如下表所示。对于无限限制，范围转换为 (0,1)，并针对此转换范围给出可选输出。令 M 为输入参数限制，令 K 为 infodict[‘last’]。条目是：

‘neval’：

函数评估的数量。

‘last’：

细分过程中产生的子区间数 K。

‘alist’：

一个长度为 M 的 rank-1 数组，其前 K 个元素是积分范围分区中子区间的左端点。

‘blist’：

一个长度为 M 的 rank-1 数组，其中的前 K 个元素是子区间的右端点。

‘rlist’：

长度为 M 的 rank-1 数组，其中的前 K 个元素是子区间上的积分近似值。

‘elist’：

长度为 M 的 rank-1 数组，其中的前 K 个元素是子区间上的绝对误差估计的模。

‘iord’：

一个长度为 M 的 rank-1 整数数组，如果为 K<=M/2+2 或 L=M+1-K，则其中的前 L 个元素是指向带有 L=K 的子间隔上的错误估计的指针。让 I 为序列 infodict['iord'] ，让 E 为序列 infodict['elist'] 。然后E[I[1]], ..., E[I[L]] 形成一个递减序列。

如果提供了输入参数点(即，它不是无)，则以下附加输出将放置在输出字典中。假设点序列的长度为 P。

‘pts’：

长度为 P+2 的 rank-1 数组，包含积分限制和按升序排列的区间断点。这是一个数组，给出了积分将发生的子区间。

‘level’：

长度为 M (=limit) 的 rank-1 整数数组，包含子区间的细分级别，即如果 (aa,bb) 是 (pts[1], pts[2]) 的子区间，其中 pts[0], pts[2] 是 infodict['pts'] 的相邻元素，则(aa,bb) 如果 |bb-aa| = |pts[2]-pts[1]| * 2**(-l) 具有级别 l。

‘ndin’：

长度为 P+2 的 rank-1 整数数组。在对区间 (pts[1], pts[2]) 进行第一次积分后，一些区间的误差估计可能已经人为地增加了，以推进它们的细分。该数组在与发生这种情况的子间隔相对应的插槽中包含一个。

加权被积函数

输入变量 weight 和 wvar 用于通过选择的函数列表对被积函数进行加权。不同的积分方法用于计算这些加权函数的积分，这些方法不支持指定断点。权重的可能值和相应的权重函数是。

weight	使用的权重函数	wvar
‘cos’	cos(w*x)	wvar = w
‘sin’	罪(w*x)	wvar = w
‘alg’	g(x) = ((x-a)**alpha)*((b-x)**beta)	wvar =(阿尔法，贝塔)
‘alg-loga’	g(x)*log(x-a)	wvar =(阿尔法，贝塔)
‘alg-logb’	g(x)*log(b-x)	wvar =(阿尔法，贝塔)
‘alg-log’	g(x)*log(x-a)*log(b-x)	wvar =(阿尔法，贝塔)
‘cauchy’	1/(x-c)	wvar = c

wvar 保存参数 w、(alpha, beta) 或 c，具体取决于所选的权重。在这些表达式中，a 和 b 是积分极限。

对于‘cos’ 和‘sin’ 加权，可以使用额外的输入和输出。

对于有限积分限制，积分使用Clenshaw-Curtis 方法执行，该方法使用切比雪夫矩。对于重复计算，这些时刻保存在输出字典中：

‘momcom’：

已计算的最大切比雪夫矩水平，即，如果 M_c 为 infodict['momcom']，则已计算长度为 |b-a| * 2**(-l) 、 l=0,1,...,M_c 的区间的矩。

‘nnlog’：

长度为 M(=limit) 的 rank-1 整数数组，包含子区间的细分级别，即，如果相应的子区间为 |b-a|* 2**(-l)，则此数组的元素等于 l。

‘chebmo’：

包含计算出的 Chebyshev 矩的 rank-2 形状数组 (25, maxp1)。通过将此数组作为序列 wopts 的第二个元素传递并将 infodict[‘momcom’] 作为第一个元素传递，这些可以在相同的时间间隔内传递给积分。

如果积分极限之一是无穷大，则计算傅里叶积分(假设 w neq 0)。如果 full_output 为 1 并且遇到数字错误，除了附加到输出元组的错误消息之外，还会将一个字典附加到输出元组，该字典将数组 info['ierlst'] 中的错误代码转换为英语消息。输出信息字典包含以下条目，而不是‘last’, ‘alist’, ‘blist’, ‘rlist’和‘elist’：

‘lst’：

集成所需的子区间数(称为 K_f )。

‘rslst’：

长度为 M_f=limlst 的 rank-1 数组，其第一个 K_f 元素包含区间 (a+(k-1)c, a+kc) 的积分贡献，其中 c = (2*floor(|w|) + 1) * pi / |w| 和 k=1,2,...,K_f 。

‘erlst’：

长度为 M_f 的 rank-1 数组，包含与 infodict['rslist'] 中相同位置的间隔对应的误差估计。

‘ierlst’：

长度为 M_f 的 rank-1 整数数组，包含与 infodict['rslist'] 中相同位置的间隔相对应的错误标志。有关代码的含义，请参阅解释字典(输出元组中的最后一个条目)。

QUADPACK关卡例程的详细信息

quad调用 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\) 的那些区间。

实变量的复函数积分

实数变量的复值函数 \(f\) 可以写为 \(f = g + ih\) 。类似地，\(f\) 的积分可以写为

\[\int_a^b f(x) dx = \int_a^b g(x) dx + i\int_a^b h(x) dx\]

假设 \(g\) 和 \(h\) 的积分存在于区间 \([a,b]\) [2] 上。因此，quad通过分别集成实部和虚部来集成complex-valued函数。

参考：

[1]

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

[2]

麦卡洛，托马斯；菲利普斯·基思 (1973)。复平面分析基础。霍尔特·莱因哈特·温斯顿。国际标准书号 0-03-086370-8

例子：

计算\(\int^4_0 x^2 dx\)并与解析结果进行比较

>>> from scipy import integrate
>>> import numpy as np
>>> x2 = lambda x: x**2
>>> integrate.quad(x2, 0, 4)
(21.333333333333332, 2.3684757858670003e-13)
>>> print(4**3 / 3.)  # analytical result
21.3333333333

计算\(\int^\infty_0 e^{-x} dx\)

>>> invexp = lambda x: np.exp(-x)
>>> integrate.quad(invexp, 0, np.inf)
(1.0, 5.842605999138044e-11)

计算 \(a = 1, 3\) 的 \(\int^1_0 a x \,dx\)

>>> f = lambda x, a: a*x
>>> y, err = integrate.quad(f, 0, 1, args=(1,))
>>> y
0.5
>>> y, err = integrate.quad(f, 0, 1, args=(3,))
>>> y
1.5

使用 ctypes 计算 \(\int^1_0 x^2 + y^2 dx\)，保持 y 参数为 1：

testlib.c =>
    double func(int n, double args[n]){
        return args[0]*args[0] + args[1]*args[1];}
compile to library testlib.*

from scipy import integrate
import ctypes
lib = ctypes.CDLL('/home/.../testlib.*') #use absolute path
lib.func.restype = ctypes.c_double
lib.func.argtypes = (ctypes.c_int,ctypes.c_double)
integrate.quad(lib.func,0,1,(1))
#(1.3333333333333333, 1.4802973661668752e-14)
print((1.0**3/3.0 + 1.0) - (0.0**3/3.0 + 0.0)) #Analytic result
# 1.3333333333333333

请注意，与积分间隔的大小相比，脉冲形状和其他尖锐特征可能无法使用此方法正确积分。此限制的一个简化示例是将 y 轴反射阶跃函数与积分范围内的许多零值积分。

>>> y = lambda x: 1 if x<=0 else 0
>>> integrate.quad(y, -1, 1)
(1.0, 1.1102230246251565e-14)
>>> integrate.quad(y, -1, 100)
(1.0000000002199108, 1.0189464580163188e-08)
>>> integrate.quad(y, -1, 10000)
(0.0, 0.0)