вектор тройного интеграла - PullRequest
Новая
       11

вектор тройного интеграла

0 голосов
/ 05 декабря 2018

Я хочу интегрировать трехмерную функцию.Давайте определим функцию как: 

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

в качестве простого примера.Давайте выберем домен, который будет [0,1]x[0,2]x[0,3]. По словам Вольфрамальфы, желаемый результат интеграции - 18,5. Вот первое, что я попробовал.Я создаю трехмерный тензор оценок pi (x, y, z), затем делаю 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.Это приводит к неожиданному несоответствию формы, связанному с первым вызовом 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)

Итак, я в замешательстве.Как я могу выполнить желаемую интеграцию?

1 Ответ

0 голосов
/ 05 декабря 2018

В вашем примере с 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)
...