四边形#
- 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中的技术,从 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
数组的长度,其中包含xx[0] == x
,其余项是包含在 quad 的args
参数中的数字。此外,某些 ctypes 调用签名支持向后兼容,但这些不应在新代码中使用。
- a浮动
积分下限(使用 -numpy.inf 表示负无穷大)。
- b浮动
积分上限(使用 numpy.inf 表示 +无穷大)。
- 参数tuple, 可选
传递给 func 的额外参数。
- 完整输出int, 可选
返回一个集成信息字典。如果不为零,警告信息也会被抑制,并且该信息会被附加到输出元组中。
- complex_funcbool, 可选
指示函数 (func) 的返回类型是实数 (
complex_func=False
: 默认) 还是复数 (complex_func=True
)。在这两种情况下,函数的参数都是实数。如果 full_output 也是非零的,则 infodict、message 和 explain 的实部和虚部将作为字典返回,键为 “real output” 和 “imag output”。
- 返回:
- y浮动
从 a 到 b 的 func 的积分。
- abserr浮动
结果中绝对误差的估计。
- infodictdict
包含附加信息的字典。
- 消息
一个收敛消息。
- 解释
仅附加了 ‘cos’ 或 ‘sin’ 加权和无限积分限,它包含了对 infodict[‘ierlst’] 中代码的解释
- 其他参数:
- epsabs浮点数或整数,可选
绝对误差容限。默认值为 1.49e-8。
quad
尝试获得abs(i-result) <= max(epsabs, epsrel*abs(i))
的精度,其中i
是从 a 到 b 的 func 的积分,而result
是数值近似值。参见下面的 epsrel。- epsrel浮点数或整数,可选
相对误差容限。默认值为 1.49e-8。如果
epsabs <= 0
,epsrel 必须大于 5e-29 和50 * (机器精度)
。参见上面的 epsabs。- 限制浮点数或整数,可选
自适应算法中使用的子区间数量的上限。
- 点(浮点数序列, 整数序列), 可选
在有界积分区间中可能出现被积函数局部困难(例如,奇点、不连续点)的断点序列。该序列不必排序。请注意,此选项不能与
weight
一起使用。- 重量浮点数或整数,可选
指示权重函数的字符串。关于此参数及剩余参数的完整解释可在下方找到。
- wvar可选的
用于加权函数的变量。
- wopts可选的
可选输入以重用切比雪夫矩。
- maxp1浮点数或整数,可选
切比雪夫矩数量的上限。
- limlstint, 可选
用于正弦加权和无限端点的循环次数的上限(>=3)。
参见
dblquad
双重积分
tplquad
三重积分
nquad
n 维积分(递归使用
quad
)fixed_quad
固定顺序的高斯求积
simpson
采样数据的积分器
romb
采样数据的积分器
scipy.special
对于正交多项式的系数和根
注释
为了得到有效的结果,积分必须收敛;对于发散积分的处理行为不保证。
quad() 输入和输出的额外信息
如果 full_output 非零,那么第三个输出参数 (infodict) 是一个字典,其条目如下表所示。对于无限极限,范围被转换为 (0,1),可选输出根据此转换范围给出。设 M 为输入参数 limit,设 K 为 infodict[‘last’]。条目如下:
- ‘neval’
函数评估的次数。
- ‘最后’
细分过程中产生的子区间数量,K。
- ‘alist’
长度为 M 的一维数组,其中前 K 个元素是积分范围分区的子区间的左端点。
- ‘blist’
长度为 M 的一维数组,其前 K 个元素是子区间的右端点。
- ‘rlist’
长度为 M 的一维数组,前 K 个元素是子区间上的积分近似值。
- ‘elist’
长度为 M 的一维数组,前 K 个元素是子区间上绝对误差估计的模数。
- ‘iord’
一个长度为 M 的一维整数数组,前 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 级数组,包含积分限和区间断点,按升序排列。这是一个数组,给出了将进行积分的子区间。
- ‘级别’
一个长度为 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)
‘柯西’
1/(x-c)
wvar = c
wvar 根据所选权重保存参数 w、(alpha, beta) 或 c。在这些表达式中,a 和 b 是积分限。
对于 ‘cos’ 和 ‘sin’ 加权,可以使用额外的输入和输出。
对于有限积分限,积分使用Clenshaw-Curtis方法进行,该方法使用Chebyshev矩。对于重复计算,这些矩保存在输出字典中:
- ‘momcom’
已经计算的切比雪夫矩的最大级别,即,如果
M_c
是infodict['momcom']
,那么矩已经为长度为|b-a| * 2**(-l)
的区间计算,l=0,1,...,M_c
。- ‘nnlog’
一个长度为 M(=limit) 的秩为 1 的整数数组,包含子区间的细分级别,即,该数组的元素等于 l,如果对应的子区间是
|b-a|* 2**(-l)
。- ‘切比莫’
一个形状为 (25, maxp1) 的秩-2数组,包含计算得到的切比雪夫矩。这些矩可以通过将此数组作为序列 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
‘柯西’
不
不
以下为每个例程从 [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: 一个用于自动积分的子程序包。Springer-Verlag. ISBN 978-3-540-12553-2.
[2]McCullough, Thomas; Phillips, Keith (1973). 《复平面分析基础》. 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)