scipy.integrate.quad处理返回numpy数组的函数

Question

我有一个函数，它返回一个numpy数组，我需要对这个数组进行积分。我想使用 scipy.integrate.dblquad，但这个方法要求函数返回一个 浮点数。

我尝试在 dblquad 上使用 numpy.vectorize，但没有成功。经过一番思考，我意识到这个方法的问题在于，向量化的 dblquad 需要的是一个函数的向量，而不是一个返回向量的函数。

我想用其他方法解决这个问题，除了将数组拆分并逐个积分（我已经试过了，代码搞得一团糟）。我首先想到的是把返回数组的函数转换成返回等效函数的数组，但我不知道这样做是否可行。

我会提供一些示例代码。下面，x 和 y 是numpy数组，a 和 b 是 浮点数。

import numpy as np
from scipy.integrate import dblquad

def foo(a, b, x, y):
    return np.exp(a*x + b*y)

我想要做的是

x = np.random.rand((3,3))
y = np.random.rand((3,3))
a = 1.5
b = 3.0
I = dblquad(foo, 0, 1, lambda x: 0, lambda x: 1, args=(x,y))

以便得到一个形状为 (3,3) 的数组 I，其中每个元素 I[i][j] 是对
foo(a,b,x[i][j],y[i][j]) 在 0 和 1 之间的积分。

我相信可以用 np.vectorize 智能地做到这一点，但我想不出任何办法。

Answer 1

这是我用过的方法：

f = lambda x,y,a,b: np.exp(a*x+b*y)
dblquad(f, 0, 1, lambda x:0, lambda x:1, args=(1.0,1.0), epsabs=1.49e-08, epsrel=1.49e-08)

f 是一个二维函数，x 和 y 是变量，a 和 b 只是参数。所以我们把 x 从零积分到一——这在 dblquad 函数中是第二和第三个参数，然后把 y 从零积分到一——在 dblquad 中是第四和第五个参数，但要以函数的形式传入；接着你把 a 和 b 以 args=(1.0,1.0) 的形式传入。其他的都是默认设置。

如果你想让参数在某个范围内变化，可以这样做：

[dblquad(f, 0, 1, lambda x:0, lambda x:1, args=(a,b)) for a in range(2) for b in range(2)]

这样你就能得到一个平坦的数组，里面是对应于 (a,b) 的结果（积分）：

 [(0, 0), (0, 1), (1, 0), (1, 1)]

你也可以这样做：

  [[dblquad(f, 0, 1, lambda x:0, lambda x:1, args=(a,b)) for a in range(2)] for b in range(2)]

这样会得到一个二维结构的结果列表，对应于 (a,b)：

  [[(0, 0), (0, 1)], [(1, 0), (1, 1)]]