使用scipy.integrate.tplquad计算插值器三重积分遇停滞问题求助

问题根源 #

NearestNDInterpolator生成的是阶梯状的非光滑函数，每个采样点对应一个Voronoi胞元，函数在胞元内取常数值。而scipy.integrate.tplquad采用的是自适应数值积分，这类积分器会不断细分区间来捕捉函数的不连续点——你的采样点多达16000个，意味着积分域内存在大量跳变，积分器会陷入无限细分的循环，导致计算时间爆炸。

解决方案 #

根据你的采样点结构（看起来是规则网格），优先选择以下优化方案：

方案1：改用规则网格插值器（最优） #

从coords输出看，x、y、z轴都是均匀间隔的规则网格（x从0.4到15.6，步长0.8；y同理；z从-15.6到15.6），这种场景下NearestNDInterpolator是大材小用，改用RegularGridInterpolator能大幅提升插值和积分效率：

import numpy as np
from scipy.interpolate import RegularGridInterpolator
from scipy.integrate import tplquad

# 提取规则网格的坐标轴
x = np.unique(coords[:, 0])
y = np.unique(coords[:, 1])
z = np.unique(coords[:, 2])

# 将一维density数组重塑为三维网格（注意维度顺序需与坐标轴匹配）
# 若coords的扁平化顺序不同，需调整reshape参数
density_3d = density.values.reshape(len(x), len(y), len(z))

# 创建规则网格最近邻插值器
interp_regular = RegularGridInterpolator((x, y, z), density_3d, method='nearest')

# 执行三重积分（tplquad参数顺序：z限→y限函数→x限函数）
result, error = tplquad(
    interp_regular,
    z1, z2,
    lambda z: y1, lambda z: y2,
    lambda z, y: x1, lambda z, y: x2
)

方案2：蒙特卡洛积分（适合不规则网格） #

如果采样点是不规则分布，无法使用规则网格插值器，蒙特卡洛积分是替代自适应积分的高效选择：通过在积分域内生成大量随机点，用采样值的平均值乘以域体积来近似积分，避免了自适应积分对不连续函数的低效细分：

import numpy as np

# 定义积分域范围
x_range = (x1, x2)
y_range = (y1, y2)
z_range = (z1, z2)

# 生成随机采样点（数量越多精度越高，1e5个点通常足够）
sample_count = 100000
x_samples = np.random.uniform(*x_range, sample_count)
y_samples = np.random.uniform(*y_range, sample_count)
z_samples = np.random.uniform(*z_range, sample_count)

# 计算插值器在采样点上的值
sample_values = interp(np.column_stack((x_samples, y_samples, z_samples)))

# 计算积分结果：平均值 × 积分域体积
domain_volume = (x_range[1]-x_range[0]) * (y_range[1]-y_range[0]) * (z_range[1]-z_range[0])
mc_result = domain_volume * np.mean(sample_values)

方案3：Voronoi胞元精确积分（复杂但精确） #

最近邻插值的积分本质是每个采样点的Voronoi胞元在积分域内的体积 × 对应density值的总和。如果需要精确结果，可以用scipy.spatial.Voronoi计算每个胞元，再裁剪到积分域内计算体积，但实现复杂度较高，仅推荐对精度要求极高的场景。

关键注意点 #

• 优先验证采样点是否为规则网格：这是效率提升最显著的优化，规则网格插值器的速度比NearestNDInterpolator快10~100倍。
• 蒙特卡洛积分的精度可通过增加采样点数量提升，1e5个点的计算时间通常在几秒内，精度足够大多数工程场景。

内容的提问来源于stack exchange，提问作者adawg