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)[source]#

计算定积分。

使用 Fortran 库 QUADPACK 中的技术，对 func 从 a 到 b（可能是无限区间）进行积分。

参数:

func{function, 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 数组的长度，其中包含 xx[0] == x，其余项是 quad 的 args 参数中包含的数字。

此外，出于向后兼容性的考虑，支持某些 ctypes 调用签名，但新代码中不应使用这些签名。

afloat

积分下限（对于负无穷大，使用 -numpy.inf）。

bfloat

积分上限（对于正无穷大，使用 numpy.inf）。

argstuple, optional

要传递给 func 的额外参数。

full_outputint, optional

非零值表示返回一个包含积分信息的字典。如果非零，则还会抑制警告消息，并将该消息附加到输出元组。

complex_funcbool, optional

指示函数的 (func) 返回类型是实数 (complex_func=False: 默认）还是复数 (complex_func=True)。在这两种情况下，函数的参数都是实数。如果 full_output 也非零，则会将实部和虚部的 infodict、message 和 explain 返回到一个包含键“real output”和“imag output”的字典中。

返回:

yfloat: func 从 a 到 b 的积分。
abserrfloat: 对结果中绝对误差的估计。
infodictdict: 包含附加信息的字典。
message: 收敛消息。
explain: 仅在使用“cos”或“sin”加权和无限积分限时附加，其中包含对 infodict[‘ierlst’] 中代码的解释。

其他参数:

epsabsfloat or int, optional: 绝对误差容差。默认为 1.49e-8。 quad 尝试获得 abs(i-result) <= max(epsabs, epsrel*abs(i)) 的精度，其中 i = func 从 a 到 b 的积分，而 result 是数值近似值。请参阅下面的 epsrel。
epsrelfloat or int, optional: 相对误差容差。默认为 1.49e-8。如果 epsabs <= 0，则 epsrel 必须大于 5e-29 和 50 * (machine epsilon)。请参阅上面的 epsabs。
limitfloat or int, optional: 自适应算法中使用的子区间的数量上限。
points(sequence of floats,ints), optional: 有界积分区间中的断点序列，在这些断点处可能会出现被积函数的局部困难（例如，奇点、不连续点）。该序列不必排序。请注意，此选项不能与 weight 结合使用。
weightfloat or int, optional: 指示加权函数的字符串。此选项以及其余参数的完整解释可在下面找到。
wvaroptional: 与加权函数一起使用的变量。
woptsoptional: 用于重用 Chebyshev 矩的可选输入。
maxp1float or int, optional: Chebyshev 矩的数量上限。
limlstint, optional: 在使用正弦加权和无限端点时，用于循环次数的上限（>=3）。

另请参阅

dblquad: 二重积分
tplquad: 三重积分
nquad: n 维积分（递归使用 quad）
fixed_quad: 固定阶高斯求积
simpson: 用于采样数据的积分器
romb: 用于采样数据的积分器
scipy.special: 用于正交多项式的系数和根

注释

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

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

如果 full_output 非零，则第三个输出参数 (infodict) 将是一个字典，其中包含如下所列的条目。对于无限限，范围将转换为 (0,1)，并且可选输出将针对此转换后的范围给出。令 M 为输入参数 limit，令 K 为 infodict[‘last’]。条目为

‘neval’: 函数评估次数。
‘last’: 在细分过程中产生的子区间数量 K。
‘alist’: 长度为 M 的秩 1 数组，其中前 K 个元素是积分范围分区的子区间的左端点。
‘blist’: 长度为 M 的秩 1 数组，其中前 K 个元素是子区间的右端点。
‘rlist’: 长度为 M 的秩 1 数组，其中前 K 个元素是子区间上的积分近似值。
‘elist’: 长度为 M 的秩 1 数组，其中前 K 个元素是子区间上绝对误差估计值的模。
‘iord’: 长度为 M 的秩 1 整数数组，其中前 L 个元素是指向子区间上误差估计值的指针，其中 L=K 如果 K<=M/2+2 或者 L=M+1-K 否则。令 I 为序列 infodict['iord']，令 E 为序列 infodict['elist']。则 E[I[1]], ..., E[I[L]] 形成一个递减序列。

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

‘pts’: 长度为 P+2 的秩 1 数组，包含积分限和区间的断点，按升序排列。这是一个数组，给出了将执行积分的子区间。
‘level’: 一个长度为 M（=limit）的秩 1 整数数组，包含子区间的细分级别，即，如果 (aa,bb) 是 (pts[1], pts[2]) 的子区间，其中 pts[0] 和 pts[2] 是 infodict['pts'] 的相邻元素，那么 (aa,bb) 的级别为 l，如果 |bb-aa| = |pts[2]-pts[1]| * 2**(-l)。
‘ndin’: 一个长度为 P+2 的秩 1 整数数组。在对区间 (pts[1], pts[2]) 进行第一次积分后，某些区间的误差估计可能被人工提高，以使其细分向前推进。此数组在对应于发生此情况的子区间的插槽中具有 1。

对被积函数进行加权

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

`weight`	使用的权重函数	`wvar`
‘cos’	cos(w*x)	wvar = w
‘sin’	sin(w*x)	wvar = w
‘alg’	g(x) = ((x-a)*alpha)((b-x)**beta)	wvar = (alpha, beta)
‘alg-loga’	g(x)*log(x-a)	wvar = (alpha, beta)
‘alg-logb’	g(x)*log(b-x)	wvar = (alpha, beta)
‘alg-log’	g(x)log(x-a)log(b-x)	wvar = (alpha, beta)
‘cauchy’	1/(x-c)	wvar = c

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

对于 'cos' 和 'sin' 加权，可以使用其他输入和输出。

对于有限积分限，积分使用 Clenshaw-Curtis 方法执行，该方法使用 Chebyshev 矩。对于重复计算，这些矩保存在输出字典中

‘momcom’: 已计算的 Chebyshev 矩的最大级别，即，如果 M_c 是 infodict['momcom']，则已为长度为 |b-a| * 2**(-l) 的区间计算了矩，l=0,1,...,M_c。
‘nnlog’: 一个长度为 M（=limit）的秩 1 整数数组，包含子区间的细分级别，即，如果相应的子区间为 |b-a|* 2**(-l)，则此数组中的一个元素等于 l。
‘chebmo’: 一个形状为 (25, maxp1) 的秩 2 数组，包含计算出的 Chebyshev 矩。这些可以传递到同一区间上的积分，方法是将此数组作为 wopts 序列的第二个元素传递，并将 infodict[‘momcom’] 作为第一个元素传递。

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

‘lst’: 积分所需的子区间数（称为 K_f）。
‘rslst’: 一个长度为 M_f=limlst 的秩 1 数组，其前 K_f 个元素包含区间 (a+(k-1)c, a+kc) 上的积分贡献，其中 c = (2*floor(|w|) + 1) * pi / |w| 且 k=1,2,...,K_f。
‘erlst’: 一个长度为 M_f 的秩 1 数组，包含与 infodict['rslist'] 中相同位置的区间对应的误差估计。
‘ierlst’: 一个长度为 M_f 的秩 1 整数数组，包含与 infodict['rslist'] 中相同位置的区间对应的错误标志。有关代码含义，请参见解释字典（输出元组中的最后一个条目）。

QUADPACK 级别例程的详细信息

quad 调用 FORTRAN 库 QUADPACK 中的例程。本节提供了有关调用每个例程的条件以及每个例程的简短描述的详细信息。调用的例程取决于 weight、points 以及积分限 a 和 b。

QUADPACK 例程	weight	points	无限边界
qagse	无	否	否
qagie	无	否	是
qagpe	无	是	否
qawoe	‘sin’、‘cos’	否	否
qawfe	‘sin’、‘cos’	否	a 或 b
qawse	‘alg*’	否	否
qawce	‘cauchy’	否	否

以下内容从 [1] 中提供了每个例程的简短描述。

qagse

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

qagie

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

qagpe

与 QAGS 的目的相同，但还允许用户提供有关疑难点的明确信息，即被积函数的内部奇点、不连续点和其他困难的横坐标。

qawoe

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

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

qawfe

计算傅里叶变换 \(\int^\infty_a \cos(\omega x)f(x)dx\) 或 \(\int^\infty_a \sin(\omega x)f(x)dx\)，对于用户提供的 \(\omega\) 和 \(f\)。在连续的有限区间上应用 QAWO 的过程，并通过 \(\varepsilon\)-算法对积分近似值序列进行收敛加速。

qawse

近似 \(\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\) 的类型。应用全局自适应细分策略，在包含 a 或 b 的那些子区间上使用修改后的 Clenshaw-Curtis 积分。

qawce

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

实变量复函数的积分

实变量的复值函数 \(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 通过分别积分实部和虚部来积分复值函数。

参考文献

[1]

Piessens, Robert; de Doncker-Kapenga, Elise; Überhuber, Christoph W.; Kahaner, David (1983). QUADPACK: A subroutine package for automatic integration. Springer-Verlag. ISBN 978-3-540-12553-2.

[2]

McCullough, Thomas; Phillips, Keith (1973). Foundations of Analysis in the Complex Plane. Holt Rinehart Winston. ISBN 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)

计算 \(\int^1_0 a x \,dx\)，其中 \(a = 1, 3\)

>>> 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)