我正在尝试使用n阶高斯积分计算双积分(在(0,0),(0,1),(1,0)处的节点上的三角形).但是,跑步
import scipy.integrate as integrate
f = lambda x,y: x+y
inside = lambda x: integrate.fixed_quad(f, 0, 1-x, args=(x,), n=5)[0]
outside = integrate.fixed_quad(inside, 0, 1, n=5)
给
Traceback (most recent call last):
File “”, line 1, in
File “/Users/username/anaconda/lib/python3.5/site-packages/scipy/integrate/quadrature.py”, line 82, in fixed_quad
return (b-a)/2.0 * np.sum(w*func(y, *args), axis=0), None
File “”, line 1, in
File “/Users/username/anaconda/lib/python3.5/site-packages/scipy/integrate/quadrature.py”, line 78, in fixed_quad
if np.isinf(a) or np.isinf(b):ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()
这是问题Can scipy.integrate.fixed_quad compute integral with functional boundaries?的第二部分.
最佳答案 在某些情况下,您的问题的答案是肯定的.
出于演示目的,我首先选择不同的界限(11而不是1 – x).
通常,可以使用dblquad解决这些类型的积分:
area_dblquad = integrate.dblquad(lambda x, y: x + y, 0, 1, lambda x: 0, lambda x: 11)[0]
这里返回66.这不是你在评论中提到的选项.
现在可以逐步进行这种集成,它适用于quad和fixed_quad:
def integrand(x, y):
return x + y
def fint_quad(x):
return integrate.quad(integrand, 0, 11, args=(x, ))[0]
def fint_fixed_quad(x):
return integrate.fixed_quad(integrand, 0, 11, args=(x, ), n=5)[0]
res_quad = integrate.quad(lambda x: fint_quad(x), 0, 1)
res_fixed_quad = integrate.fixed_quad(lambda x: fint_fixed_quad(x), 0, 1, n=5)
正如预期的那样,他们都回归66.这表明它可以使用scipy.integrate.fixed_quad计算双积分.
但是,当一个人现在将上限更改回你所拥有的那个(从11到1 – x)时,它仍适用于quad但是为fixed_quad崩溃:
area_dblquad = integrate.dblquad(lambda x, y: x + y, 0, 1, lambda x: 0, lambda x: 1 - x)[0]
res_quad = integrate.quad(lambda x: fint_quad(x), 0, 1)
两者都返回0.333333 …,使用fixed_quad调用会导致您收到错误.通过查看源代码可以理解错误:
x, w = _cached_roots_legendre(n)
x = np.real(x)
if np.isinf(a) or np.isinf(b):
raise ValueError("Gaussian quadrature is only available for "
"finite limits.")
y = (b-a)*(x+1)/2.0 + a
return (b-a)/2.0 * np.sum(w*func(y, *args), axis=-1), None
当一个打印a和b一个得到:
a: 0
b: 1
a: 0
b: [ 0.95308992 0.76923466 0.5 0.23076534 0.04691008]
因此,对于使用1-x的调用,b实际上是具有n个元素的numpy数组,并且无法将数组与无穷大进行比较,这就是它崩溃的原因.无论是预期的行为还是错误,我无法回答;可能值得在github上打开一个问题.