是否可以使用 scipy.integrate.fixed_quad 计算二重积分？

Question

我正在尝试使用 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()

这是问题的第二部分。

Answer 1

你的问题的答案是，是的，在特定条件下。

出于演示目的，我首先选择与您不同的边界（11 而不是 1 - x）。

通常，可以使用 dblquad:

求解这些类型的积分

area_dblquad = integrate.dblquad(lambda x, y: x + y, 0, 1, lambda x: 0, lambda x: 11)[0]

这里 returns 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)

正如预期的那样，他们都 return 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)

return0.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.

上开一个问题

Answer 2

fixed_quad 要求 f 接受向量输入。结果应该是输入的映射值（即 map(f, xs) 之类的东西）。考虑到这一点，只需确保您的 inside 函数 returns 映射值，您就可以开始了。

import scipy.integrate as integrate
f = lambda y,x: x+y
inside = lambda xs, n: np.array([integrate.fixed_quad(f, 0, 1-x, args=(x,), n=n)[0]
                                 for x in np.array(xs)])
order = 5
outside = integrate.fixed_quad(inside, 0, 1, n=order, args=(order,))

另外，请注意被积函数的参数顺序。从您的代码 arg=(x,) 来看，您似乎希望沿 y 维度进行内积分。被积函数的第一个参数是它被积分的维度。所以它应该是lambda y,x（注意这也是dblquad期望的被积函数的形状）。

是否可以使用 scipy.integrate.fixed_quad 计算二重积分？

Is it possible to compute double integral using scipy.integrate.fixed_quad?

python

scipy

numerical-integration