scipy.integrate.quad 关于返回 numpy 数组的函数

Question

我有一个函数返回一个需要集成的 numpy 数组。我想使用scipy.integrate.dblquad，但它需要函数返回一个float。

我尝试在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))

得到一个数组I，形状为(3,3)，入口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 的 2D 函数 - 只是参数。所以我们在 dblquad 函数中整合 x 从零到一 - 第二和第三个 args，然后在 dblquad 中整合从零到一 - 第 4 和第 5 个 args 的 y，但是以函数的形式；然后你以 args=(1.0,1.0) 的形式传递你的 a,b。其余的都是默认的东西

如果您想在某个值范围内改变参数，您可以执行以下操作：

[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)]]