3D空间任意三角形上f=1/||x-y||的Python高效积分实现

核心结论 #

SciPy完全适用，但需要根据奇点x是否在三角形上选择不同的积分策略：

• 若x在三角形外部：直接用普通数值积分即可收敛
• 若x在三角形内部/边上：需处理奇点，通过坐标变换消除奇异性后再积分

方法详解与代码实现 #

第一步：三角形参数化

将三角形内任意点y用参数u, v表示：

def triangle_point(u, v, v1, v2, v3):
    return v1 + u*(v2 - v1) + v*(v3 - v1)

对应的面积元缩放因子（雅可比行列式的模长）为三角形两条边叉乘的模长：

def jacobian(v1, v2, v3):
    edge1 = v2 - v1
    edge2 = v3 - v1
    return np.linalg.norm(np.cross(edge1, edge2))

情况1：奇点在三角形外部

此时被积函数在积分区域内连续，直接使用scipy.integrate.dblquad计算二重积分：

import numpy as np
from scipy.integrate import dblquad

x = np.array([0, 0, 1])
v1 = np.array([0, 1, 0])
v2 = np.array([1, 0, 0])
v3 = np.array([0, 0, 0])

# 定义被积函数（带雅可比因子）
def integrand(v, u, x, v1, v2, v3):
    y = triangle_point(u, v, v1, v2, v3)
    return jacobian(v1, v2, v3) / np.linalg.norm(x - y)

# 积分限：u从0到1，v从0到1-u
integral, error = dblquad(
    integrand,
    0, 1,
    lambda u: 0, lambda u: 1 - u,
    args=(x, v1, v2, v3)
)

print(f"积分结果：{integral:.6f}，误差估计：{error:.6e}")

情况2：奇点在三角形内部/边上

此时被积函数在积分区域内有奇点，直接积分会导致精度差或不收敛，需要做极坐标变换消除奇异性：

1. 将坐标原点移到奇点x，把三角形转化为相对于x的多边形
2. 用极坐标(r, θ)替换原参数，此时被积函数的奇异性被消除，只需对θ和r积分

以下是实现代码：

import numpy as np
from scipy.integrate import quad

def integrate_singular_triangle(x, v1, v2, v3):
    # 将三角形顶点转换为相对于x的坐标
    p1 = v1 - x
    p2 = v2 - x
    p3 = v3 - x

    # 计算三角形的边向量与顶点
    edges = [p2-p1, p3-p2, p1-p3]
    vertices = [p1, p2, p3]

    # 计算所有顶点的极角并统一到[0, 2π)区间
    thetas = [np.arctan2(p[1], p[0]) if (p[0] !=0 or p[1] !=0) else 0 for p in vertices]
    thetas = [theta % (2*np.pi) for theta in thetas]

    # 按极角排序顶点与对应角度
    sorted_indices = np.argsort(thetas)
    sorted_thetas = [thetas[i] for i in sorted_indices]
    sorted_vertices = [vertices[i] for i in sorted_indices]

    total_integral = 0.0

    # 遍历每个极角区间计算积分
    for i in range(3):
        theta1 = sorted_thetas[i]
        theta2 = sorted_thetas[(i+1)%3]
        p_a = sorted_vertices[i]
        p_b = sorted_vertices[(i+1)%3]

        # 推导当前边的极坐标r(θ)表达式
        a = p_b[1] - p_a[1]
        b = p_a[0] - p_b[0]
        c = p_a[0]*p_b[1] - p_a[1]*p_b[0]

        def r(theta):
            denom = a * np.cos(theta) + b * np.sin(theta)
            return c / denom if denom != 0 else np.inf

        # 处理θ跨越2π的情况
        if theta2 < theta1:
            part1, _ = quad(r, theta1, 2*np.pi)
            part2, _ = quad(r, 0, theta2)
            segment_integral = part1 + part2
        else:
            segment_integral, _ = quad(r, theta1, theta2)

        total_integral += segment_integral

    return total_integral

# 测试奇点在三角形内部的情况
x_inside = np.array([0.2, 0.2, 0])
integral_singular = integrate_singular_triangle(x_inside, v1, v2, v3)
print(f"奇点在内部时的积分结果：{integral_singular:.6f}")

关键说明 #

• 奇点在三角形上时，极坐标变换将原奇异积分转化为普通积分，避免了数值发散问题
• dblquad和quad返回的误差估计可用于判断结果可靠性
• 预计算三角形几何参数（边向量、叉乘模长）可减少重复计算，提升效率

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