基于Python的多面体非矩形域三重积分计算及误差优化

一、嵌套一维积分实现（对应R语言方案） #

使用scipy.integrate.quad（Python常用一维积分函数）嵌套实现，逻辑和R版本一致。需用numpy.vectorize处理内层积分的矢量化调用，适配外层积分的变量遍历要求：

from scipy.integrate import quad
import numpy as np

def f(x, y, z):
    return x + y * z

# 最内层：对z积分
def inner_z(x, y):
    val, _ = quad(f, -10, 6 - x - y, args=(x, y))
    return val

# 中间层：对y积分，矢量化内层函数以适配quad输入
def middle_y(x):
    vec_inner = np.vectorize(lambda y: inner_z(x, y))
    val, _ = quad(vec_inner, -5, 3 - x)
    return val

# 最外层：对x积分
final_val, outer_err = quad(middle_y, -5, 4)
print(f"嵌套积分结果: {final_val:.1f}，外层积分估计误差: {outer_err}")

该方案和R版本存在相同局限：仅能获取最外层积分的误差估计，前两层积分的误差未被纳入总误差计算，整体误差结果可靠性不足。

二、带可靠误差估计的解决方案：使用nquad #

scipy.integrate.nquad是专为多重积分设计的函数，会自动处理各层积分的误差传播，输出整体可信的误差估计。只需按积分层级定义变量范围（内层变量范围可依赖外层变量）：

from scipy.integrate import nquad

# 调整参数顺序为(z, y, x)，匹配积分层级（z最内层，x最外层）
def integrand(z, y, x):
    return x + y * z

# 定义各变量积分范围：从内层到外层依次为z、y、x
limits = [
    [-10, lambda y, x: 6 - x - y],  # z的上限依赖y和x
    [-5, lambda x: 3 - x],          # y的上限依赖x
    [-5, 4]                         # x的固定范围
]

result, total_error = nquad(integrand, limits)
print(f"nquad积分结果: {result:.1f}，整体估计误差: {total_error}")

nquad的核心优势：它会将每一层积分的误差按统计学方法传播到外层，最终输出的total_error是覆盖所有层级误差的整体估计，比嵌套quad的误差结果更可靠。

内容的提问来源于stack exchange，提问作者Stéphane Laurent