Julia中基于离散数组的三维数值积分实现问询

针对你拥有离散数组phi1::Vector{Float64}（对应三维均匀网格x,y,z各100个点）的纯数值积分需求，无需解析函数，直接基于网格特性实现即可，无需依赖HCubature这类面向连续函数的积分包。以下是两种实用方案：

1. 矩形积分法（黎曼和） #

这是最直接的离散积分方式，适用于均匀网格：

• 先计算单个网格单元的体积
• 将一维数组重塑为三维网格数组，确保索引与x,y,z对应
• 所有网格点值求和后乘以单元体积

# 假设x、y、z已定义为range(-1, stop=1, length=100)
dx = step(x)
dy = step(y)
dz = step(z)
dV = dx * dy * dz

# 把一维phi1转为三维数组，匹配x→y→z的网格顺序
phi3d = reshape(phi1, length(x), length(y), length(z))

# 计算积分结果
integral_rect = sum(phi3d) * dV

2. 复合梯形法则（更高精度） #

比矩形积分精度更高，适合均匀网格的数值积分，核心是给网格边界点赋予0.5的权重：

# 生成梯形法则的权重数组
function trapezoidal_weights(n)
    weights = ones(n)
    weights[1] = weights[end] = 0.5
    return weights
end

wx = trapezoidal_weights(length(x))
wy = trapezoidal_weights(length(y))
wz = trapezoidal_weights(length(z))

# 计算加权后的网格值总和
weighted_sum = sum(phi3d .* wx .* wy' .* reshape(wz, 1, 1, :))

# 最终积分结果
integral_trap = weighted_sum * dx * dy * dz

关键注意事项 #

• 确认phi1的存储顺序与x,y,z的网格遍历顺序一致，保证phi3d[i,j,k]对应phi(x[i], y[j], z[k])
• step()函数比直接取x[2]-x[1]更通用，适配非均匀range场景（此处为均匀网格）
• 若追求更高精度，可使用复合辛普森法则，但要求每个维度的网格点数为奇数（当前100为偶数，可调整网格长度或改用辛普森3/8法则）

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