寻求支持四面体三维数值积分的Python库（支持点输入）

• Scipy + 仿射变换适配
  主流积分库如Scipy的tplquad虽要求函数形式的边界，但可以通过仿射变换将任意四面体映射到标准四面体（顶点为(0,0,0), (1,0,0), (0,1,0), (0,0,1)），从而实现基于顶点输入的积分计算，步骤如下：

  1. 用四个顶点构造变换向量：取其中一个顶点为原点基准，其余三个顶点与基准点的向量作为变换基；
  2. 计算变换的雅可比行列式绝对值，作为积分的权重系数；
  3. 将原多项式函数替换为变换后的表达式，再在标准四面体的固定积分限内计算积分，最后乘以雅可比行列式得到结果。

  示例代码：

  import numpy as np
  from scipy.integrate import tplquad
  
  def tetra_integrate(p0, p1, p2, p3, func):
      # 构造变换基向量
      v1 = np.array(p1) - np.array(p0)
      v2 = np.array(p2) - np.array(p0)
      v3 = np.array(p3) - np.array(p0)
      # 计算雅可比行列式
      jacobian = abs(np.linalg.det(np.array([v1, v2, v3])))
      
      # 定义变换后的目标函数
      def transformed(u, v, w):
          x = p0[0] + u*v1[0] + v*v2[0] + w*v3[0]
          y = p0[1] + u*v1[1] + v*v2[1] + w*v3[1]
          z = p0[2] + u*v1[2] + v*v2[2] + w*v3[2]
          return func(x, y, z)
      
      # 标准四面体的积分限：u∈[0,1], v∈[0,1-u], w∈[0,1-u-v]
      result, _ = tplquad(transformed,
                          0, 1,
                          lambda u: 0, lambda u: 1 - u,
                          lambda u, v: 0, lambda u, v: 1 - u - v)
      return result * jacobian
  
  # 测试示例：积分f(x,y,z)=x+y+z，顶点为标准四面体
  test_func = lambda x,y,z: x + y + z
  p0, p1, p2, p3 = (0,0,0), (1,0,0), (0,1,0), (0,0,1)
  print(tetra_integrate(p0, p1, p2, p3, test_func))  # 预期输出：0.25（即1/4）
  

• Sympy 符号精确积分
  如果需要多项式积分的解析精确结果，可以用Sympy做符号化处理，思路同样是通过仿射变换将任意四面体转为标准区域，再计算符号积分，适合需要精确表达式的场景。

  示例代码：

  import sympy as sp
  
  # 定义符号变量
  u, v, w = sp.symbols('u v w')
  x0, y0, z0 = sp.symbols('x0 y0 z0')
  x1, y1, z1 = sp.symbols('x1 y1 z1')
  x2, y2, z2 = sp.symbols('x2 y2 z2')
  x3, y3, z3 = sp.symbols('x3 y3 z3')
  
  # 构造仿射变换
  v1 = (x1-x0, y1-y0, z1-z0)
  v2 = (x2-x0, y2-y0, z2-z0)
  v3 = (x3-x0, y3-y0, z3-z0)
  x = x0 + u*v1[0] + v*v2[0] + w*v3[0]
  y = y0 + u*v1[1] + v*v2[1] + w*v3[1]
  z = z0 + u*v1[2] + v*v2[2] + w*v3[2]
  
  # 示例多项式函数
  poly_func = x + y + z
  # 计算雅可比行列式
  jacobian = abs(sp.Matrix([[v1[0], v2[0], v3[0]], 
                            [v1[1], v2[1], v3[1]], 
                            [v1[2], v2[2], v3[2]]]).det())
  
  # 计算符号积分
  integral_expr = sp.integrate(
      sp.integrate(
          sp.integrate(poly_func * jacobian, (w, 0, 1-u-v)),
          (v, 0, 1-u)
      ),
      (u, 0, 1)
  )
  
  # 代入具体顶点值求值
  result = integral_expr.subs({
      x0:0, y0:0, z0:0,
      x1:1, y1:0, z1:0,
      x2:0, y2:1, z2:0,
      x3:0, y3:0, z3:1
  })
  print(result)  # 输出：1/4
  

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