python - 如何使用 scipy.integrate.fixed_quad 一次计算多个积分？

Question

给定一个函数func(x,y,z)，我想提供一个函数

def integral_over_z(func,x,y,zmin=0,zmax=1,n=16):
    lambda_func = z,x,y: ???
    return scipy.integrate.fixed_quad(lambda_func,a=zmin,b=zmax,args=(x,y),n=n)

使用scipy.integrate.fixed_quad计算用户提供的(x,y)输入在z上的积分。输入 (x,y) 可以是单个浮点或 float 组(当两者都是数组时，它们的形状相同)。

scipy.integrate.fixed_quad支持向量值函数的积分。为此，函数 func 必须返回相应的更高维度的数组:“如果集成向量值函数，则返回的数组必须具有形状 (..., len(x) )"(来自文档)。

因此，我的问题是如何生成 lambda_func 的相应输出数组(可以使用特殊用途的类来实现)。

编辑:为了帮助理解我的问题，这里有一个可行的实现，但没有通过 z 进行矢量化(因此不使用 scipy.integrate.fixed_quad )。

def integral_over_z(func,x,y,zmin,zmax,n=16):
    z,w = scipy.special.roots_legendre(n)
    dz = 0.5*(zmax-zmin)
    z = zmin + (np.real(z)+1) * dz
    w = np.real(w) * dz
    result = w[0] * func(x,y,z[0])
    for i in range(1,len(z)):
        result += w[i] * func(x,y,z[i])
    return result

问题是:如何对其进行矢量化，以便它适用于任何有效输入(x 和/或 y float 或数组)。

另一个编辑: 对于通过 scipy.integrate.fixed_quad 实现，被积函数必须采用形状为 (nz) 的 z 一维数组。输入x和y必须一起广播，当它们的广播形状可以是任何东西时，比如(n0,n1,..,nk) 那么从 func 返回的值必须具有形状 (n0,n1,..,nk,nz) ——我如何生成它？

Answer 1

它似乎是一个向量值函数，向量值必须位于第 0 维，并且积分参数(在您的情况下 z)必须放在最后(这就是它们的含义 ( ..., len(x))，他们的x是你的z)，我认为这来自广播规则。以下示例对我来说效果很好 - 这里的关键是 x 和 y 必须具有正确的形状才能使广播工作

import numpy as np
import scipy.integrate
def integral_over_z(func,x,y,n=16):
    lambda_func = lambda z, x, y: func(x[..., None],y[..., None],z)  # the last dimension of (x,y) needs to be size 1, but you can have as many leading dimensions as you want
    return scipy.integrate.fixed_quad(lambda_func,a=0,b=1,args=(x,y),n=n)
func = lambda x,y,z: 1 + 0*x + 0*y + 0*z  # make sure that the output has the right (broadcast) shape
x = np.zeros((5,))
y =  np.arange(5)
print(integral_over_z(func, x, y, 2))