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{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_datascipy.LowLevelCallable 中包含的数据。在使用 xx 的调用形式中,nxx 数组的长度,该数组包含 xx[0] == x,其余项是 quad 的 args 参数中包含的数字。

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

a浮点数

积分下限(-无穷大请使用 -numpy.inf)。

b浮点数

积分上限(+无穷大请使用 numpy.inf)。

args元组, 可选

要传递给 func 的额外参数。

full_output整数, 可选

非零时返回一个包含积分信息的字典。如果为非零,还会抑制警告消息,并将消息附加到输出元组中。

complex_func布尔值, 可选

指示函数 (func) 的返回类型是实数 (complex_func=False: 默认) 还是复数 (complex_func=True)。在两种情况下,函数的参数都是实数。如果 full_output 也为非零,则实部和虚部的 infodictmessageexplain 将在以“real output”和“imag output”为键的字典中返回。

返回:
y浮点数

func 从 ab 的积分。

abserr浮点数

结果的绝对误差估计。

infodict字典

包含附加信息的字典。

message

收敛消息。

explain

仅在使用“cos”或“sin”加权和无限积分限制时附加,它包含对 infodict['ierlst'] 中代码的解释。

其他参数:
epsabs浮点数或整数, 可选

绝对误差容差。默认值为 1.49e-8。quad 试图达到精度 abs(i-result) <= max(epsabs, epsrel*abs(i)),其中 i = funcab 的积分,result 是数值近似。参见下面的 epsrel

epsrel浮点数或整数, 可选

相对误差容差。默认值为 1.49e-8。如果 epsabs <= 0,则 epsrel 必须大于 5e-29 和 50 * (machine epsilon)。参见上面的 epsabs

limit浮点数或整数, 可选

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

points(浮点数或整数序列), 可选

有界积分区间内的一系列断点,在这些断点处,被积函数可能出现局部困难(例如,奇异点、不连续点)。该序列无需排序。请注意,此选项不能与 weight 一起使用。

weight浮点数或整数, 可选

指示加权函数的字符串。此参数和其余参数的完整解释可在下面找到。

wvar可选

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

wopts可选

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

maxp1浮点数或整数, 可选

切比雪夫矩数量的上限。

limlst整数, 可选

与正弦加权和无限端点一起使用的循环次数(>=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 的一维数组,其前 K 个元素是积分范围划分中子区间的左端点。

‘blist’

一个长度为 M 的一维数组,其前 K 个元素是子区间的右端点。

‘rlist’

一个长度为 M 的一维数组,其前 K 个元素是子区间上的积分近似值。

‘elist’

一个长度为 M 的一维数组,其前 K 个元素是子区间上绝对误差估计的模。

‘iord’

一个长度为 M 的一维整数数组,其前 L 个元素是指向子区间误差估计的指针,如果 K<=M/2+2L=K,否则 L=M+1-K。令 I 为序列 infodict['iord'],令 E 为序列 infodict['elist']。则 E[I[1]], ..., E[I[L]] 构成一个递减序列。

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

‘pts’

一个长度为 P+2 的一维数组,包含积分限制和区间的断点,按升序排列。这是一个给出将进行积分的子区间的数组。

‘level’

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

‘ndin’

一个长度为 P+2 的一维整数数组。在对区间 (pts[1], pts[2]) 进行第一次积分后,为了推进其细分,某些区间上的误差估计可能被人为地增加了。此数组在发生这种情况的子区间对应的槽位中包含 1。

对被积函数进行加权

输入变量 weightwvar 用于通过选定的函数列表对被积函数进行加权。使用不同的积分方法来计算带有这些加权函数的积分,并且这些方法不支持指定断点。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。在这些表达式中,a 和 b 是积分限制。

对于“cos”和“sin”加权,提供了额外的输入和输出。

对于具有有限积分限制的加权积分,使用 Clenshaw-Curtis 方法执行积分,该方法使用切比雪夫矩。对于重复计算,这些矩会保存在输出字典中

‘momcom’

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

‘nnlog’

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

‘chebmo’

一个形状为 (25, maxp1) 的二维数组,包含计算出的切比雪夫矩。可以通过将此数组作为序列 wopts 的第二个元素,并将 infodict[‘momcom’] 作为第一个元素,将其传递给相同区间上的积分。

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

‘lst’

积分所需的子区间数量(称之为 K_f)。

‘rslst’

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

‘erlst’

一个长度为 M_f 的一维数组,包含与 infodict['rslist'] 中相同位置的区间对应的误差估计。

‘ierlst’

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

QUADPACK 级别例程的详细信息

quad 调用 FORTRAN 库 QUADPACK 中的例程。本节提供了调用每个例程的条件以及每个例程的简短描述。调用的例程取决于 weightpoints 以及积分限制 ab

QUADPACK 例程

加权

points

无限边界

qagse

None

qagie

None

qagpe

None

qawoe

‘sin’, ‘cos’

qawfe

‘sin’, ‘cos’

a 或 b

qawse

‘alg*’

qawce

‘cauchy’

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

qagse

是一个基于全局自适应区间细分结合外推的积分器,它将消除多种被积函数奇异性的影响。积分使用每个子区间内的 21 点高斯-克朗罗德求积法执行。

qagie

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

qagpe

服务于与 QAGS 相同的目的,但还允许用户提供有关问题点位置和类型的明确信息,即内部奇异点、不连续点和被积函数其他困难的横坐标。

qawoe

是一个积分器,用于评估有限区间 [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)\) 中的奇异性。

qawfe

计算用户提供的 \(\omega\)\(f\) 的傅里叶变换 \(\int^\infty_a \cos(\omega x)f(x)dx\)\(\int^\infty_a \sin(\omega x)f(x)dx\)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\) 的类型。采用全局自适应细分策略,并在包含 ab 的子区间上使用改进的 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)