我一直在研究 Python 中的 Particle-In-Cell 框架，并注意到一个问题，即 EM 脉冲的幅度随着模拟的更新而无限增加。
目前，我正在使用本文中概述的非标准有限差分算法：JLVay 2008。
这是用于在真空中推进电场和磁场的算法：
@n.njit(parallel=True)
def E_Update(E, B, J, ds, dt, R):
C = 2.99792458e8
EPSILON = 8.8541878128e-12
MU = 1.25663706212e-6
for i in n.prange(1,E.shape[0]-1):
for j in n.prange(1,E.shape[1]-1):
for k in n.prange(1,E.shape[2]-1):
# x-component update
# -> Solve DyBz
Bz_y0 = B[i,j-1,k,2]
Bz_y1 = B[i,j,k,2]
DyBz = (Bz_y1-Bz_y0)/ds
# -> Solve DzBy
By_z0 = B[i,j,k-1,1]
By_z1 = B[i,j,k,1]
DzBy = (By_z1-By_z0)/ds
# -> Solve Ex
R[i,j,k,0] = (1/EPSILON)*((DyBz-DzBy) - J[i,j,k,0])*dt + E[i,j,k,0]
# y-component update
# -> Solve DxBz
Bz_x0 = B[i-1,j,k,2]
Bz_x1 = B[i,j,k,2]
DxBz = (Bz_x1-Bz_x0)/ds
# -> Solve DzBx
Bx_z0 = B[i,j,k-1,0]
Bx_z1 = B[i,j,k,0]
DzBx = (Bx_z1-Bx_z0)/ds
# -> Solve Ey
R[i,j,k,1] = (1/EPSILON)*((-1)*(DyBz-DzBy) - J[i,j,k,1])*dt + E[i,j,k,1]
# z-component update
# -> Solve DxBy
By_x0 = B[i-1,j,k,1]
By_x1 = B[i,j,k,1]
DxBy = (By_x1-By_x0)/ds
# -> Solve DyBx
Bx_y0 = B[i,j-1,k,0]
Bx_y1 = B[i,j,k,0]
DyBx = (Bx_y1-Bx_y0)/ds
# -> Solve Ez
R[i,j,k,2] = (1/EPSILON)*((DxBy-DyBx) - J[i,j,k,2])*dt + E[i,j,k,2]
@n.njit(parallel=True)
def B_Update(E, B, ds, dt, a, b, g, R):
C = 2.99792458e8
EPSILON = 8.8541878128e-12
MU = 1.25663706212e-6
C1 = ((1-(dt/(2*MU)))/(1+(dt/(2*MU))))
C2 = (1/(1+(dt/(2*MU))))
for i in n.prange(1,E.shape[0]-1):
for j in n.prange(1,E.shape[1]-1):
for k in n.prange(1,E.shape[2]-1):
# x-component update
# -> Solve DzEy
Ey_z0 = E[i,j,k,1]
Ey_z1 = E[i,j,k+1,1]
DzEy = (Ey_z1-Ey_z0)/ds
# -> Solve S1zEy
S1zEyQ1 = E[i+1,j,k,1] # Ey(i+1,j,k)
S1zEyQ2 = E[i-1,j,k,1] # Ey(i-1,j,k)
S1zEyQ3 = E[i,j+1,k,1] # Ey(i,j+1,k)
S1zEyQ4 = E[i,j-1,k,1] # Ey(i,j-1,k)
S1zEy = S1zEyQ1 + S1zEyQ2 + S1zEyQ3 + S1zEyQ4
# -> Solve S2zEy
S2zEyQ1 = E[i+1,j+1,k,1] # Ey[i+1,j+1,k]
S2zEyQ2 = E[i+1,j-1,k,1] # Ey[i+1,j-1,k]
S2zEyQ3 = E[i-1,j+1,k,1] # Ey[i-1,j+1,k]
S2zEyQ4 = E[i-1,j-1,k,1] # Ey[i-1,j-1,k]
S2zEy = S2zEyQ1 + S2zEyQ2 + S2zEyQ3 + S2zEyQ4
# -> Solve DSzEy
DSzEy = (a+b*S1zEy+g*S2zEy)*(DzEy)
# -> Solve DyEz
Ez_y0 = E[i,j,k,2]
Ez_y1 = E[i,j+1,k,2]
DyEz = (Ez_y1-Ez_y0)/ds
# -> Solve S1yEz
S1yEzQ1 = E[i+1,j,k,2] # Ez(i+1,j,k)
S1yEzQ2 = E[i-1,j,k,2] # Ez(i-1,j,k)
S1yEzQ3 = E[i,j,k+1,2] # Ez(i,j,k+1)
S1yEzQ4 = E[i,j,k-1,2] # Ez(i,j,k-1)
S1yEz = S1yEzQ1 + S1yEzQ2 + S1yEzQ3 + S1yEzQ4
# -> Solve S2yEz
S2yEzQ1 = E[i+1,j,k+1,2] # Ez(i+1,j,k+1)
S2yEzQ2 = E[i+1,j,k-1,2] # Ez(i+1,j,k-1)
S2yEzQ3 = E[i-1,j,k+1,2] # Ez(i-1,j,k+1)
S2yEzQ4 = E[i-1,j,k-1,2] # Ez(i-1,j,k-1)
S2yEz = S2yEzQ1 + S2yEzQ2 + S2yEzQ3 + S2yEzQ4
# -> Solve DSyEz
DSyEz = (a+b*S1yEz+g*S2yEz)*(DyEz)
# -> Solve DtBx
R[i,j,k,0] = (1/MU)*((DSzEy-DSyEz))*dt + B[i,j,k,0]
# y-component update
# -> Solve DzEx
Ex_z0 = E[i,j,k,0]
Ex_z1 = E[i,j,k+1,0]
DzEx = (Ex_z1-Ex_z0)/ds
# -> Solve S1zEx
S1zExQ1 = E[i+1,j,k,0] # Ex(i+1,j,k)
S1zExQ2 = E[i-1,j,k,0] # Ex(i-1,j,k)
S1zExQ3 = E[i,j+1,k,0] # Ex(i,j+1,k)
S1zExQ4 = E[i,j-1,k,0] # Ex(i,j-1,k)
S1zEx = S1zExQ1 + S1zExQ2 + S1zExQ3 + S1zExQ4
# -> Solve S2zEx
S2zExQ1 = E[i+1,j+1,k,0] # Ex(i+1,j+1,k)
S2zExQ2 = E[i+1,j-1,k,0] # Ex(i+1,j-1,k)
S2zExQ3 = E[i-1,j+1,k,0] # Ex(i-1,j+1,k)
S2zExQ4 = E[i-1,j-1,k,0] # Ex(i-1,j-1,k)
S2zEx = S2zExQ1 + S2zExQ2 + S2zExQ3 + S2zExQ4
# -> Solve DSzEx
DSzEx = (a+b*S1zEx+g*S2zEx)*(DzEx)
# -> Solve DxEz
Ez_x0 = E[i,j,k,2]
Ez_x1 = E[i+1,j,k,2]
DxEz = (Ez_x1-Ez_x0)/ds
# -> Solve S1xEz
S1xEzQ1 = E[i,j+1,k,2] # Ez(i,j+1,k)
S1xEzQ2 = E[i,j-1,k,2] # Ez(i,j-1,k)
S1xEzQ3 = E[i,j,k+1,2] # Ez(i,j,k+1)
S1xEzQ4 = E[i,j,k-1,2] # Ez(i,j,k-1)
S1xEz = S1xEzQ1 + S1xEzQ2 + S1xEzQ3 + S1xEzQ4
# -> Solve S2xEz
S2xEzQ1 = E[i,j+1,k+1,2] # Ez(i,j+1,k+1)
S2xEzQ2 = E[i,j+1,k-1,2] # Ez(i,j+1,k-1)
S2xEzQ3 = E[i,j-1,k+1,2] # Ez(i,j-1,k+1)
S2xEzQ4 = E[i,j-1,k-1,2] # Ez(i,j-1,k-1)
S2xEz = S2xEzQ1+S2xEzQ2+S2xEzQ3+S2xEzQ4
# -> Solve DSxEz
DSxEz = (a+b*S1xEz+g*S2xEz)*(DxEz)
# -> Solve DtBy
R[i,j,k,1] = (1/MU)*((DSxEz-DSzEx))*dt + B[i,j,k,1]
# z-component update
# -> Solve DyEx
Ex_y0 = E[i,j,k,0]
Ex_y1 = E[i,j+1,k,0]
DyEx = (Ex_y1-Ex_y0)/ds
# -> Solve S1yEx
S1yExQ1 = E[i+1,j,k,0] # Ex(i+1,j,k)
S1yExQ2 = E[i-1,j,k,0] # Ex(i-1,j,k)
S1yExQ3 = E[i,j,k+1,0] # Ex(i,j,k+1)
S1yExQ4 = E[i,j,k-1,0] # Ex(i,j,k-1)
S1yEx = S1yExQ1 + S1yExQ2 + S1yExQ3 + S1yExQ4
# -> Solve S2yEx
S2yExQ1 = E[i+1,j,k+1,0] # Ex(i+1,j,k+1)
S2yExQ2 = E[i+1,j,k-1,0] # Ex(i+1,j,k-1)
S2yExQ3 = E[i-1,j,k+1,0] # Ex(i-1,j,k+1)
S2yExQ4 = E[i-1,j,k-1,0] # Ex(i-1,j,k-1)
S2yEx = S2yExQ1 + S2yExQ2 + S2yExQ3 + S2yExQ4
# -> Solve DSyEx
DSyEx = (a+b*S1yEx+g*S2yEx)*(DyEx)
# -> Solve DxEy
Ey_x0 = E[i,j,k,1]
Ey_x1 = E[i+1,j,k,1]
DxEy = (Ey_x1-Ey_x0)/ds
# -> Solve S1xEy
S1xEyQ1 = E[i,j+1,k,1] # Ey(i,j+1,k)
S1xEyQ2 = E[i,j-1,k,1] # Ey(i,j-1,k)
S1xEyQ3 = E[i,j,k+1,1] # Ey(i,j,k+1)
S1xEyQ4 = E[i,j,k-1,1] # Ey(i,j,k-1)
S1xEy = S1xEyQ1 + S1xEyQ2 + S1xEyQ3 + S1xEyQ4
# -> Solve S2xEy
S2xEyQ1 = E[i,j+1,k+1,1] # Ey(i,j+1,k+1)
S2xEyQ2 = E[i,j+1,k-1,1] # Ey(i,j+1,k-1)
S2xEyQ3 = E[i,j-1,k+1,1] # Ey(i,j-1,k+1)
S2xEyQ4 = E[i,j-1,k-1,1] # Ey(i,j-1,k-1)
S2xEy = S2xEyQ1 + S2xEyQ2 + S2xEyQ3 + S2xEyQ4
# -> Solve DSxEy
DSxEy = (a+b*S1xEy+g*S2xEy)*(DxEy)
# -> Solve DtBz
R[i,j,k,2] = (1/MU)*((DSyEx-DSxEy))*dt + B[i,j,k,2]
每个更新步骤都涉及通过，然后电场由，最后 B-Field 再次通过.
到目前为止，我已尝试通过以下方式解决此问题：
- 恢复到标准 FDTD（即将 B_Update 方法中的 a、b 和 g 分别设置为 1、0 和 0）
- 通过保持时间步长比穿过 Yee-Cell 长度的光的穿越时间低几个数量级来确保满足 CFL 条件。
- 使用自然单位进行模拟（
我发现的唯一涉及 FDTD 模拟“爆炸”问题的资源是关于算法基础的 YouTube 讲座。FDTD 方法的稳定性取决于在适当的时间和空间坐标处更新方程中使用的某些值的存在。我可能是错的，但我相信我写的更新函数满足了这个稳定性要求。
任何建议的更改或资源建议将不胜感激。