高斯立方体电通量计算异常：Riemann积分结果不符问题排查

核心错误点 #

1. 微分面积向量da定义错误 #

立方体顶面（z=1）的面元法向量沿+z方向，微分面积向量应为np.array([0, 0, dx*dy])，你代码中用x、y、z缩放的形式完全违背了面元向量的物理定义，直接导致点积结果失效。

2. 电场强度公式错误 #

原点处点电荷的电场公式应为：
$$\vec{E} = \frac{q \vec{r}}{r^3}$$
其中$r = \sqrt{x^2 + y^2 + z^{2}$，你仅除以$r}2$，遗漏了分母的$r$，导致电场强度量级偏差，这是结果数值异常的关键原因。

3. 积分时重复计算面积元素 #

代码中flux += dot(Evector(x,y,1), da(x,y,1,dx,dy,0)) * dx * dy，但正确的da本身已包含面元面积$dx*dy$，重复相乘会让结果多一个面积平方量级，进一步放大误差。

修正后的代码 #

import numpy as np

q = 1  # 1 statocoulomb

def dot(vec0, vec1):
    return vec0 @ vec1

# 立方体顶面（z=1）的微分面积向量，法向量沿+z方向
def da_top(dx, dy):
    return np.array([0, 0, dx * dy])

def Evector(x, y, z):
    r = np.sqrt(x**2 + y**2 + z**2)
    r3 = r ** 3
    return np.array([x, y, z]) * (q / r3)

# 顶面积分范围
xmin = -1
xmax = 1
ymin = -1
ymax = 1

# 分割数（增大可提高精度）
xstepmax = 1000
ystepmax = 1000

dx = (xmax - xmin) / xstepmax
dy = (ymax - ymin) / ystepmax

flux = 0.0

for xsteps in range(xstepmax):
    # 中点采样避免边界误差
    x = xmin + dx * (xsteps + 0.5)
    for ysteps in range(ystepmax):
        y = ymin + dy * (ysteps + 0.5)
        flux += dot(Evector(x, y, 1), da_top(dx, dy))

print(f"顶面电通量计算值: {flux}")
print(f"高斯定律理论值(4π/6): {4 * np.pi / 6}")

结果说明 #

修正后，当分割数足够大时（如1000×1000），计算结果会趋近于理论值$4\pi/6 \approx 2.094$，符合预期。

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