基于FDM的1D电荷电势泊松方程Python求解故障排查

问题描述 #

用Python有限差分法(FDM)求解一维电荷电势分布，泊松方程形式为A×C=Rho：

• A应为三对角矩阵（主对角线元素为-2，上下次对角线为1，实际离散后需除以dX²）
• Rho为含点电荷的源项数组

使用np.linalg.solve(A, Rho)得到的结果绘图异常，调整dX、X范围增加计算点数无效；仅将A设为整数类型时结果看似正常，但会出现取整偏差（如-199替代-200），怀疑Rho或矩阵构造存在问题。

完整原始代码 #

import numpy as np
import matplotlib.pyplot as plt

X = 30
dX = 0.1
N = int(X / dX)

x = np.arange(0, N)
x = x * dX

Rho = np.zeros(N)
Rho[int(10 / dX)] = -(1.6e-19) / dX ** 2
Rho[int(20 / dX)] = -(-1.6e-19) / dX ** 2

A = np.zeros((N, N), int)
np.fill_diagonal(A, -2 / dX ** 2)
A.ravel()[1::(2 * N + 2)] = 1 / dX ** 2
A.ravel()[2 + N::(2 * N + 2)] = 1 / dX ** 2
A.ravel()[N::(2 * N + 2)] = 1 / dX ** 2
A.ravel()[2 * N + 1::(2 * N + 2)] = 1 / dX ** 2
print(A)
print(Rho)
C = np.linalg.solve(A, Rho)

fig, ax = plt.subplots()
ax.plot(x, C, linewidth=1.0)
plt.show()

故障排查方向 #

1. A矩阵构造逻辑完全错误 #

原始代码用A.ravel()的索引切片赋值的方式完全不符合三对角矩阵的构造逻辑，导致矩阵结构混乱，这是求解结果异常的核心原因。正确的三对角矩阵构造应该直接操作对角线：

• 主对角线：-2 / dX²
• 下对角线（行i>0，列i-1）：1 / dX²
• 上对角线（列j>0，行j-1）：1 / dX²

2. 矩阵类型强制转换导致的隐性错误 #

原始代码初始化A时指定了int类型，赋值-2/dX²（如dX=0.1时为-200）会被强制转整数，但此时错误的构造逻辑可能误打误撞让矩阵接近预期结构；改为float类型后，错误的索引赋值会让矩阵完全偏离三对角结构，导致求解结果彻底异常。

3. Rho的符号与物理量纲错误 #

泊松方程的标准形式为∇²φ = -ρ/ε₀，有限差分离散后需匹配物理关系：

• 原始代码未引入真空介电常数ε₀=8.854e-12，导致数值量级偏差极大
• Rho的符号和计算方式错误：一维情况下点电荷的电荷密度应为Q/dX（线电荷密度），而非Q/dX²，且符号需与泊松方程的推导一致

4. 边界条件未正确设置 #

原始代码未明确边界条件，默认的矩阵第一行和最后一行结构不符合实际物理场景（如无穷远电势为0的Dirichlet边界）。正确做法需将两端点的行设置为单位向量，强制电势为0，同时修改对应的Rho值。

修正后示例代码 #

import numpy as np
import matplotlib.pyplot as plt

X = 30
dX = 0.1
N = int(X / dX)
epsilon0 = 8.854e-12  # 真空介电常数

x = np.arange(0, N) * dX

# 构造Rho：一维线电荷密度（Q/dX）
Rho = np.zeros(N)
# 负电荷位置（10m处）
Rho[int(10 / dX)] = -1.6e-19 / dX
# 正电荷位置（20m处）
Rho[int(20 / dX)] = 1.6e-19 / dX

# 构造正确的三对角矩阵A
A = np.zeros((N, N), dtype=np.float64)
np.fill_diagonal(A, -2 / dX**2)
np.fill_diagonal(A[1:], 1 / dX**2)  # 下对角线
np.fill_diagonal(A[:, 1:], 1 / dX**2)  # 上对角线

# 设置Dirichlet边界：两端电势为0
A[0, :] = 0
A[0, 0] = 1
A[-1, :] = 0
A[-1, -1] = 1
Rho[0] = 0
Rho[-1] = 0

# 匹配泊松方程形式：Aφ = -Rho/epsilon0
b = -Rho / epsilon0
C = np.linalg.solve(A, b)

fig, ax = plt.subplots()
ax.plot(x, C, linewidth=1.0)
ax.set_xlabel('位置 (m)')
ax.set_ylabel('电势 (V)')
plt.show()

内容的提问来源于stack exchange，提问作者Andrzej Sołtys