主页 / 问题 / 1598615
Accepted
Слава Погодин
Asked:2024-11-02 20:25:33 +0000 UTC2024-11-02 20:25:33 +0000 UTC

scipy.integrte.nquad 无法积分椭球

  • 772

为了工作，你需要整合区域，我从一个简单的图形开始 - 一个椭球体。 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，一切都按其应有的方式解决了。

我尝试将数字放入边界函数而不是参数变量（例如a）中，没有任何变化。如果删除一维，例如z，则一切都很重要。如果您使用球体而不是椭球体，则一切都很重要。如果增加参数b，那么从b = 1.9 .. 2之间的某个位置开始，所有内容都会开始正确计算。

这里出了什么问题以及如何将所有内容与代码中的参数集成？

python

1 个回答

  1. Best Answer

    问题在于，由于舍入，(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)

    有计算的笔记本

    • 2

相关问题