从二维电场分量求标量电势：Python与Fortran结果不符问题排查

我有两个二维电场分量数组：Eₓ（尺寸为256×512）和E_z（尺寸为256×512）。已知标量电势的梯度等于电场矢量，因此需要通过梯度逆运算计算标量电势Φ（代码中用phi表示）。

我尝试使用以下Python循环计算电势：

nx = 256
nz = 512
phi = np.zeros([nx, nz])

for i in range(1, nz):
    phi[:, i] = -dx * Ex[:, i-1] + phi[:, i-1]
    for j in range(1, nx):
        phi[j, :] = -dz * Ez[j-1, :] + phi[j-1, :]

得到的Φ结果存在异常竖线，不符合预期；而同事使用Fortran代码计算的结果符合预期，其Fortran代码如下：

allocate(phi(nx, nz))

phi(1, 1) = 0.0

do i = 2, nx
    phi(i, 1) = -dx * ex(i-1, 1) + phi(i-1, 1)
enddo

do i = 1, nx
    do j = 2, nz
        phi(i, j) = -dz * ez(i, j-1) + phi(i, j-1)
    enddo
enddo

错误排查 #

你的Python代码与Fortran代码的计算逻辑完全不一致，这是结果异常的核心原因：

• Fortran代码采用先x方向初始化、再z方向逐列递推的顺序：
  1. 先填充第一列（z=1）的所有行，利用Eₓ分量从phi(1,1)递推出phi(2:nx,1)
  2. 再对每一行，从左到右（z从2到nz）利用E_z分量递推后续列的电势值
• 你的Python代码在每一次z方向迭代中都重复执行完整的x方向全列更新：
  1. 每计算完一列z方向的电势，就立刻对所有列执行一次x方向的全量更新
  2. 这种重复更新会干扰后续z列的计算，引入错误的累积值，最终导致结果出现异常竖线

修正后的Python代码 #

完全对齐Fortran的计算逻辑，先处理第一列的x方向递推，再逐列处理z方向递推：

import numpy as np

nx = 256
nz = 512
phi = np.zeros([nx, nz])

# 第一步：计算第一列（z=0）的所有x方向电势值
for j in range(1, nx):
    phi[j, 0] = -dx * Ex[j-1, 0] + phi[j-1, 0]

# 第二步：对每一行，从左到右递推z方向的电势值
for j in range(nx):
    for i in range(1, nz):
        phi[j, i] = -dz * Ez[j, i-1] + phi[j, i-1]

更优的向量化实现 #

利用numpy的向量化运算替代显式循环，大幅提升计算效率，结果与Fortran代码完全一致：

import numpy as np

nx = 256
nz = 512
phi = np.zeros([nx, nz])

# 第一列x方向递推：累加Ex的前j项（从0到j-1）
phi[:, 0] = -dx * np.cumsum(Ex[:, 0], axis=0)

# 每一行z方向递推：累加Ez的前i项（从0到i-1），叠加第一列的初始值
phi[:, 1:] = phi[:, 0, np.newaxis] - dz * np.cumsum(Ez[:, :-1], axis=1)

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