给定一个函数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)
，在用户提供(x,y)输入时，用scipy.integrate.fixed_quad计算其对z的积分。输入(x,y)可以是单个浮点数，也可以是浮点数数组(当两者都是数组时，它们的形状相同)。
scipy.integrate.fixed_quad支持向量值函数积分。为此，函数func必须返回一个相应的高维数组:如果对一个向量值函数进行积分，返回的数组必须具有形状为(..., len(x))的数组;(来自文档).
因此，我的问题是如何生成lambda_func的相应输出数组(可以使用特殊用途的class来实现)。
编辑:为了帮助理解我的问题，这里是一个工作的实现，但不是在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浮点数或数组)。
另一个编辑:为了通过scipy.integrate.fixed_quad实现，被积函数必须采用形状为(nz)的z的1D数组。输入x和y必须一起广播，当它们广播的形状可以是任何形状，比如(n0,n1,..,nk)，那么从func返回的形状一定是(n0,n1,..,nk,nz)——我是怎么生成的?
这似乎是一个矢量值函数，矢量值必须在第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))
经过缺陷的(不完整)回答和阅读numpy广播后，我找到了一个解决方案。我很乐意了解这是否仍然可以改进和/或如果这真的是正确的，即适用于任何有效的输入(它为我的测试到目前为止)。
重要的是调整x和y的形状，使
-
func(x,y,z)工作正常，即x、y和z是可联合广播的; -
对最后一个维度(
z)的func输出求和后，得到x和y的联合广播形状。
这是我的解决方案:
def integral_over_z(func,x,y,zmin=0,zmax=1,n=16):
xe = x
ye = y
if type(xe) is np.ndarray or type(ye) is np.ndarray:
xe,ye = np.broadcast_arrays(x,y) # replace x,y by their joint broadcast
xe = np.expand_dims(xe, xe.ndim) # expand by an extra dimension for z
ye = np.expand_dims(ye, ye.ndim) # expand by an extra dimension for z
return scipy.integrate.fixed_quad(lambda z : func(xe,ye,z), a=zmin, b=zmax, n=n)