如何正确使用 scipy tplquad？答案

【问题标题】：how to use scipy tplquad properly?如何正确使用 scipy tplquad？
【发布时间】：2019-12-31 12:33:06
【问题描述】：

我正在尝试使用 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。

【问题讨论】：

标签： python scipy numerical-integration

【解决方案1】：

大多数数值积分例程（例如tplquad）使用多项式逼近积分。如果函数是平滑的，这很好用。不幸的是，特征函数但是平滑的，因为它们具有不连续的边界。这就是tplquad 失败的原因。

如果您想近似域的体积，一种合理的方法是为其创建三角形（2D）或四面体（3D）网格，然后添加单纯形的体积。网格生成器的示例是 pygmsh 和 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，并返回精确的结果。

【讨论】：

【解决方案2】：

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

解决方法是将集成域定义为单位球：

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 可能是解决更复杂问题的正确方法...

【讨论】：