如何正确使用 scipy tplquad？

Question

我正在尝试使用 scipy 的 tplquad 在物理学中做一些简单的三重积分。例如，我尝试在单位立方体上积分单位球 (func d) 的恒定质量密度。这是行不通的。但是，如果我在单位立方体上积分单位立方体 (func f) 的恒定质量密度，我会很快得到结果。

我认为问题在于提供常量积分限制作为常量而不是函数。我使用 lambda 修复了这个问题，但我仍然无法获得积分。

from scipy import integrate

''' returns the mass density at a point (x,y,z)'''
def d(z, x, y):
  return int(x**2 + y**2 + z**2 <= 1) # unit ball with constant density = 1 , here all orthogonal axes are principal 

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

integrate.tplquad(d, -1, 1, lambda x: -1, lambda x: 1,lambda x, y: -1, lambda x, y: 1) # doesn't work / too slow
integrate.tplquad(f, -1, 1, lambda x: -1, lambda x: 1,lambda x, y: -1, lambda x, y: 1) # works fine

我希望 d 在给定范围内的积分为 4/3*pi。

Answer 1

不连续函数的积分是一个困难的数值问题。

解决方法是将积分域定义为单位球：

from scipy import integrate
import numpy as np

''' returns the mass density at a point (x,y,z)'''
def d(z, x, y):
      return 1

integrate.tplquad(d, -1, 1,
                  lambda x: -np.sqrt(1-x**2), lambda x: np.sqrt(1-x**2),
                  lambda x, y: -np.sqrt(1-x**2-y**2), lambda x, y: np.sqrt(1-x**2-y**2))
# (4.188790204786397, 2.000470900043183e-09)
# 4/3*np.pi = 4.1887902047863905

另一个不理想的解决方案是人为地平滑函数。例如使用 Logistic 函数，请参阅 Analytic approximations of Heaviside_step_function:

''' returns the mass density at a point (x,y,z)'''
def d(z, x, y):
    r2 = x**2 + y**2 + z**2
    smoothing_length = 0.1 # same unit as r2
    d = 1 - 1/(1 + np.exp(-2*(r2-1)/smoothing_length))
    return d

integrate.tplquad(d, -1, 1, lambda x: -1, lambda x: 1,lambda x, y: -1, lambda x, y: 1)
# (4.182852937567993, 3.021537155780628e-08)

必须谨慎选择值 smoothing_length。

Monte Carlo integration 可能是解决更复杂问题的正确方法...

Answer 2

大多数数值积分例程（例如tplquad）使用多项式近似积分。如果功能流畅，则效果很好。不幸的是，特征函数是其他所有但光滑的，因为它们具有不连续的边界。这就是 tplquad 失败的原因。

如果您想估计域的体积，合理的方法是为其创建三角形（2D）或四面体（3D）网格，然后添加单纯形的体积。网格生成器的一个示例是 pygmsh and pygalmesh（我的一个项目），但还有其他的。

如果你真的想集成球上的功能，你应该看看 quadpy（我的另一个项目）。它具有针对不同领域的许多集成方案，其中包括球。这个

import numpy
import quadpy

scheme = quadpy.ball.hammer_stroud_14_3()
val = scheme.integrate(
    lambda x: numpy.ones_like(x[0]),  # function to integrate
    [0.0, 0.0, 0.0],  # center
    1.0,  # radius
    )
print(val)  # 4.1887902047863905 == 4*pi/3

将使用 5 次方案在单位球上对函数 1 求积分，returns 正是结果。

如何正确使用 scipy tplquad？

how to use scipy tplquad properly?

python

scipy

numerical-integration