scipy.integrte.nquad не может проинтегрировать эллипсоид

Question

для работы надо интегрировать области, начал с простой фигуры - эллипсоид. Код на Python:

import scipy.integrate as integral

x0 = 0
y0 = 0
z0 = 0

a = 2
b = 1
c = 1

def f(x, y, z):
    return 1

def bounds_z(x, y):
    return [z0 - c*(1 - ((x - x0)/a)**2 - ((y - y0)/b)**2)**0.5, z0 + c*(1 - ((x - x0)/a)**2 - ((y - y0)/b)**2)**0.5]

def bounds_y(x):
    return [y0 - b*(1 - ((x - x0)/a)**2)**0.5, y0 + b*(1 - ((x - x0)/a)**2)**0.5]

def bounds_x():
    return [x0 - a, x0 + a]

res = integral.nquad(f, [bounds_z, bounds_y, bounds_x])
print(res)

Выполнение прерывается ошибкой:

TypeError: '<' not supported between instances of 'complex' and 'complex'

Логически пределы интегрирования и их порядок поставлены правильно. Посмотрел на wolfram alpha, там все решается как надо.

Пробовал в функции bounds ставить числа вместо переменных-параметров (например a), ничего не меняется. Если убрать одну размерность, например z, все считается. Если вместо эллипсоида ставить сферу все считается. Если увеличивать параметр b, то начиная где то между b = 1.9 .. 2 начинает все вычисляться правильно.

Что тут не так и как проинтегрировать все с параметрами из кода?

Answer 1

Проблема в том, что из-за округления выражение (x/a)**2 + (y/b)**2 на границах интергрирования может стать больше 1.

Вот пример:

import numpy as np

X = np.linspace(x0-a,x0+a,100)
Y = bounds_y(X)
print(1 - (((X-x0)/a)**2 + ((Y[0]-y0)/b)**2))

Результат:

[ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
  ...
  0.00000000e+00  0.00000000e+00  0.00000000e+00 -2.22044605e-16
  0.00000000e+00  1.11022302e-16  2.22044605e-16 -2.22044605e-16
 -2.22044605e-16  0.00000000e+00  0.00000000e+00 -2.22044605e-16
  ...]

Видите отрицательное число? Когда в этой точке будет вычисляться Z по вашей формуле, получится комплексное число, и интегрирование остановится с ошибкой, что для комплексных чисел операция < не определена.

Вам нужно в bounds_z проверять, что выражение под корнем неотрицательное. Для отрицательных значений возвращать [z0,z0] в качестве границ интегрирования.

import scipy.integrate as integral
import numpy as np

x0 = 0
y0 = 0
z0 = 0

a = 2
b = 1
c = 1

def f(x, y, z):
    return 1

def bounds_z(x, y):
  D = 1 - ((x - x0)/a)**2 - ((y - y0)/b)**2
  if D < 0:
    D = 0
  d = D**0.5
  return (z0 - c*d, z0 + c*d)

def bounds_y(x):
  D = 1 - ((x - x0)/a)**2
  if D<0:
    D = 0
  d = D**0.5
  return (y0 - b*d, y0 + b*d)

def bounds_x():
    return (x0 - a, x0 + a)

Проверка на комплексные z:

X = np.linspace(x0-a,x0+a,100)
Y = np.vectorize(bounds_y)(X)
Z1 = np.vectorize(bounds_z)(X, Y[0])
print(Z1)

Видно, что комплексных чисел среди Z1 нет.

Но что-то не так с вызовом nquad: integral.nquad(f, [bounds_z, bounds_y, bounds_x]) возвращает (5.640598097619295, 6.442042087200108e-08), а должен вернуть объем эллипсоида 4.0/3.0*math.pi*a*b*c, 8.377580409572781

Зато integrate.tplquad Работает как надо:

xl = x0-a
xh = x0+a
yl = lambda x: bounds_y(x)[0]
yh = lambda x: bounds_y(x)[1]
zl = lambda x,y: bounds_z(x,y)[0]
zh = lambda x,y: bounds_z(x,y)[1]

def g(z,y,x):
  return f(x,y,z)

print(integral.tplquad(g, xl, xh, yl, yh, zl, zh))

вывод

(8.377580409572793, 2.000470900043183e-09)

Для правильного ответа в nquad нужно изменить порядок вычисления граничных значений:

range_z = (-c,c)

def range_y(z):
  D = 1 - ((z-z0)/c)**2
  d = D**0.5 if D>0 else 0
  return (y0 - b*d, y0 + b*d)

def range_x(y,z):
  D = 1 - ((z-z0)/c)**2 - ((y-y0)/b)**2
  d = D**0.5 if D>0 else 0
  return (x0 - a*d, x0 + a*d)

print(integral.nquad(f, [range_x, range_y, range_z]))

Результат

(8.377580409572793, 4.000941800086366e-09)

блокнот с вычислениями

БЛОГ НА HUSL

scipy.integrte.nquad не может проинтегрировать эллипсоид

Ответы (1 шт):