利用FDTD计算一维二维或者三维低频电磁波在海水中传播的电场和磁场强度

问题遇到的现象和发生背景

问题相关代码，请勿粘贴截图

import numpy as np
from matplotlib import pyplot as plt
from math import exp, cos, sin, sqrt, atan2, pi
ke = 400
ex = np.zeros(ke)
dx = np.zeros(ke)
ix = np.zeros(ke)
hy = np.zeros(ke)
ddx = 1 #size
dt = ddx/6e8# step size
boundary_low = [0, 0]
boundary_high = [0, 0]# Create Dielectric Profile
epsz = 8.854e-12
epsr = 74
sigma = 3
k_start = 5
freq_in = 30
gax = np.ones(ke)
gbx = np.zeros(ke)
gax[k_start:] = 1 / (epsr + (sigma * dt / epsz))
gbx[k_start:] = sigma * dt / epsz
nsteps = 20000# Main FDTD Loop
for time_step in range(1, nsteps + 1):# Calculate Dx
    for k in range(1, ke):
        dx[k] = dx[k] + 0.5 * (hy[k - 1] - hy[k])
    # Put a sinusoidal at the low end
    pulse = sin(2 * pi * freq_in * dt * time_step)
    dx[5] = pulse + dx[5]
    # Calculate the Ex field from Dx
    for k in range(1, ke):
        ex[k] = gax[k] * (dx[k] - ix[k])
        ix[k] = ix[k] + gbx[k] * ex[k]  # Absorbing Boundary Conditions
    ex[0] = boundary_low.pop(0)
    boundary_low.append(ex[1])
    ex[ke - 1] = boundary_high.pop(0)
    boundary_high.append(ex[ke - 2])# Calculate the Hy field
    for k in range(ke - 1):
        hy[k] = hy[k] + 0.5 * (ex[k] - ex[k + 1])

运行结果及报错内容

无法模拟低频正弦电磁波在海水中传播

我的解答思路和尝试过的方法

我想要达到的结果

利用FDTD计算一维二维或者三维低频电磁波在海水中传播的电场和磁场强度

是效果不好吗？这个需要检查一下公式等等的，图画出来就好分析了

无法模拟低频正弦电磁波在海水中传播,可能是FDTD算法实现时的个别步骤出现了偏差，可以，参考一下实现成功的算法步骤，进而对比参考修正。

看一下：
https://xueshu.baidu.com/usercenter/paper/show?paperid=9c5ab725045ab86fa17bb5665cf7dfcb

太难，mark 帮顶

计算结果与理论值对不上的话，应该是程序里面的计算公式部分有问题，需要一步步执行来看

无法模拟低频正弦电磁波在海水中传播,可能是FDTD算法实现时的个别步骤出现了偏差，可以，参考一下实现成功的算法步骤，进而对比参考修正。

你就是想画个图展现出来呗，我用的python37版本的

import numpy as np
from matplotlib import pyplot as plt
from math import exp, cos, sin, sqrt, atan2, pi
ke = 400
ex = np.zeros(ke)
dx = np.zeros(ke)
ix = np.zeros(ke)
hy = np.zeros(ke)
ddx = 1 #size
dt = ddx/6e8# step size
boundary_low = [0, 0]
boundary_high = [0, 0]# Create Dielectric Profile
epsz = 8.854e-12
epsr = 74
sigma = 3
k_start = 5
freq_in = 30
gax = np.ones(ke)
gbx = np.zeros(ke)
gax[k_start:] = 1 / (epsr + (sigma * dt / epsz))
gbx[k_start:] = sigma * dt / epsz
nsteps = 20000# Main FDTD Loop
for time_step in range(1, nsteps + 1):# Calculate Dx
    for k in range(1, ke):
        dx[k] = dx[k] + 0.5 * (hy[k - 1] - hy[k])
    # Put a sinusoidal at the low end
    pulse = sin(2 * pi * freq_in * dt * time_step)
    dx[5] = pulse + dx[5]
    # Calculate the Ex field from Dx
    for k in range(1, ke):
        ex[k] = gax[k] * (dx[k] - ix[k])
        ix[k] = ix[k] + gbx[k] * ex[k]  # Absorbing Boundary Conditions
    ex[0] = boundary_low.pop(0)
    boundary_low.append(ex[1])
    ex[ke - 1] = boundary_high.pop(0)
    boundary_high.append(ex[ke - 2])# Calculate the Hy field
    for k in range(ke - 1):
        hy[k] = hy[k] + 0.5 * (ex[k] - ex[k + 1])
    plt.plot(dx,hy)
plt.show()

https://download.csdn.net/download/jiebing2020/21875456?spm=1005.2026.3001.5635&utm_medium=distribute.pc_relevant_ask_down.none-task-download-2~default~OPENSEARCH~Rate-2.pc_feed_download_top3ask&depth_1-utm_source=distribute.pc_relevant_ask_down.none-task-download-2~default~OPENSEARCH~Rate-2.pc_feed_download_top3ask

https://xueshu.baidu.com/usercenter/paper/show?paperid=61aa9b373abc0f74c59aab59476f3e93