如何在Python中利用Cubature高效积分六边形区域？

你的问题核心在于当前极坐标方法中大量采样点落在六边形外，被积函数返回0，浪费了计算资源。以下是几种高效的改进方案：

方案一：利用对称性+坐标变换（最优） #

六边形具有6次旋转对称性，我们可以只计算1/6区域的积分再乘以6，同时通过坐标变换将不规则的积分限转为规则的单位正方形，让所有采样点都落在有效区域内，彻底消除无效计算。

原理 #

将原极坐标变量 (k, theta) 转换为新变量 (u, v)：

• theta = u * π/3（将theta范围从[0, π/3]映射到u∈[0,1]）
• k = v * rfun(theta)（将r的范围从[0, rfun(theta)]映射到v∈[0,1]）

同时计算积分变换的雅可比因子，确保积分值的正确性。

代码实现 #

import numpy as np
from cubature import cubature

# 六边形边界函数
def rfun(theta):
    pt = np.pi / 3
    phi = theta % pt
    denom = np.cos(phi) + np.cos(pt - phi)
    return (3 * s) / (2 * denom)

# 坐标变换后的被积函数（支持向量化）
def transformed_fun(uv_arr, para):
    # uv_arr是形状为(n,2)的数组，对应多组(u,v)采样点
    u = uv_arr[:, 0]
    v = uv_arr[:, 1]
    
    theta = u * np.pi / 3
    r = rfun(theta)
    k = v * r
    
    # --------------------------
    # 替换为你的实际被积函数计算
    # 示例：假设原函数实部和虚部的计算逻辑
    real_part = ...  # 用numpy数组操作实现向量化计算
    imag_part = ...
    # --------------------------
    
    # 计算雅可比因子（包含极坐标面积元k dk dtheta的变换）
    jacobian = v * (np.pi / 3) * r**2
    
    # 返回乘过雅可比因子的结果
    return np.column_stack([real_part * jacobian, imag_part * jacobian])

def integralfunction(para):
    global s  # 若s是外部传入参数，可将其合并到para中
    # 积分区域：u∈[0,1], v∈[0,1]
    uv_min = np.array([0.0, 0.0], np.float64)
    uv_max = np.array([1.0, 1.0], np.float64)
    
    # 开启vectorized=True，批量处理采样点
    val, err = cubature(transformed_fun, 2, 2, uv_min, uv_max, args=(para,),
                        adaptive='h', relerr=1e-5, maxEval=1e6, vectorized=True)
    
    # 乘以6得到整个六边形的积分结果
    ireal_total = val[0] * 6
    iimag_total = val[1] * 6
    
    fresult = somefunction(ireal_total, iimag_total)
    return fresult

方案二：向量化原函数（快速小改） #

如果不想做坐标变换，可以先优化原函数的向量化能力，减少Python循环开销：

def myfunction(kk_arr, para):
    # kk_arr是形状为(n,2)的数组，每行是(k, theta)
    k = kk_arr[:, 0]
    theta = kk_arr[:, 1]
    cond = rfun(theta)
    
    # 用numpy掩码处理边界判断，避免循环
    mask = k <= cond
    real_part = np.zeros_like(k)
    imag_part = np.zeros_like(k)
    
    # 仅对有效区域计算被积函数
    real_part[mask] = ...  # 你的实部计算逻辑（针对mask选中的元素）
    imag_part[mask] = ...  # 你的虚部计算逻辑
    
    return np.column_stack([real_part, imag_part])

def integralfunction(para):
    kmin = np.array([0.0, 0.0], np.float64)
    kmax = np.array([s, 2*np.pi], np.float64)
    
    # 开启vectorized=True
    val, err = cubature(myfunction, 2, 2, kmin, kmax, args=(para,),
                        adaptive='h', relerr=1e-5, maxEval=1e6, vectorized=True)
    
    ireal, iimag = val
    fresult = somefunction(ireal, iimag)
    return fresult

关键优化点 #

• 消除无效采样：方案一通过坐标变换让所有采样点都落在六边形内，彻底避免了无效计算，这是提升速度的核心。
• 向量化计算：开启vectorized=True，让函数一次性处理多组采样点，大幅减少Python函数调用的开销。
• 对称性利用：将积分区域缩小为1/6，直接减少了5/6的计算量。

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