找回密码
 立即注册
首页 python python-基础 查看内容

scipy tripple矢量积分


我想整合一个三维功能。让我们将函数定义为:
def pi(x, y, z): return x + y ** 2 + z ** 3


举个简单的例子。让我们选择域名[0,1]x[0,2]x[0,3]。根据wolframalpha,所需的积分结果是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) 

所以,我很困惑。如何执行所需的集成?

解决方法


在你的Wolfram例子中,你的积分是从内到外:从0到3整合x ,然后在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)

分享至 : QQ空间
收藏

0 个回复

您需要登录后才可以回帖 登录 | 立即注册