scipy 三元向量积分

Question

我想整合一个三维函数。让我们将函数定义为：

def pi(x, y, z): return x + y ** 2 + z ** 3

作为一个简单的例子。让我们选择域为[0,1]x[0,2]x[0,3]。 According to wolframalpha, the desired integration result is 18.5. 这是我尝试的第一件事。我创建了一个 pi(x,y,z) 评估的 3D 张量，然后进行 3 次 1D 积分：

from scipy.integrate import trapz
import numpy as np
x = np.linspace(0, 1)
y = np.linspace(0, 2)
z = np.linspace(0, 3)
print(trapz(trapz(trapz(pi(x[:, None, None], y[None, :, None], z[None, None, :]), x), y), z)) # 51.51853394418992

请注意我的输出不正确。我认为这是错误的，因为我没有正确的整合顺序。接下来我尝试的是明确引用 x、y 和 z 的 3D 张量。这会导致与第一个 trapz 调用相关的意外形状不匹配：

print(trapz(trapz(trapz(pi(x[:, None, None], y[None, :, None], z[None, None, :]), x[:, None, None]), y[None, :, None]), z[None, None, :]))

Traceback (most recent call last):
  File "/usr/local/Cellar/python/3.6.5/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/numpy/lib/function_base.py", line 4523, in trapz
    ret = (d * (y[slice1] + y[slice2]) / 2.0).sum(axis)
ValueError: operands could not be broadcast together with shapes (50,1,0) (50,50,49) 

所以，我很困惑。如何执行所需的集成？

Answer 1

在您的 Wolfram 示例中，您的积分从内到外是：在 x 上从 0 到 3 积分，然后在 y 上从 0 到 2，最后在 z 从 0 到 1。但是在您的代码中，x 从 0 到 1，z 从 0 到 3。当我输入这些不同的范围时，我得到 18.502290712203248：

def pi(x, y, z):
    return x + y ** 2 + z ** 3

x = np.linspace(0, 3) # To three!
y = np.linspace(0, 2)
z = np.linspace(0, 1) # To one!

# I broke this up here, just to make it easier for me to read and debug.
x_int = trapz(pi(x[:, None, None], y[None, :, None], z[None, None, :]), x)
y_int = trapz(x_int, y)
z_int = trapz(y_int, z)
print(z_int)