高阶FDTD的MARLAB 代码

这里提供一个使用MATLAB实现的高阶FDTD代码示例：

% 定义常数
c0 = 3e8; % 光速
dx = 5e-8; % 空间离散步长
dt = dx / (2 * c0); % 时间离散步长
tmax = 2e-14; % 最大模拟时间
l = 1.55e-6; % 光波长
f = c0 / l; % 光频率
w = 2 * pi * f; % 光角频率
eps0 = 8.85e-12; % 真空介电常数
mu0 = 4 * pi * 1e-7; % 真空磁导率
epsr = 12; % 介质相对介电常数
sigma = 0; % 介质电导率

% 定义计算区域
N = 200; % 区域长度
M = round(N * l / dx); % 区域宽度
x = 0:dx:(N - 1) * dx;
y = 0:dx:(M - 1) * dx;

% 初始化电磁场
Ey = zeros(N, M);
Hz = zeros(N, M);

% 计算系数
eps = eps0 * epsr * ones(N, M);
mu = mu0 * ones(N, M);
sigmam = sigma * ones(N, M); % 磁场电导率
sigmae = sigma * ones(N, M); % 电场电导率
a = dt / eps0;
b = dt / mu0;
cme = (1 - sigmam * dt / (2 * mu)) ./ (1 + sigmam * dt / (2 * mu));
cmm = (1 - sigmae * dt / (2 * eps)) ./ (1 + sigmae * dt / (2 * eps));
dte = 1 ./ (eps .* (1 + sigmae * dt / (2 * eps)));
dtm = 1 ./ (mu .* (1 + sigmam * dt / (2 * mu)));

% 定义高阶FDTD算法
for t = 0:dt:tmax
    % 更新电场
    for j = 1:M
        for i = 1:N-1
            Ey(i,j) = cme(i,j) * Ey(i,j) + cmm(i,j) * (Hz(i,j) - Hz(i+1,j)) / dx;
        end
    end
    
    % 更新磁场
    for j = 1:M-1
        for i = 1:N
            Hz(i,j) = cme(i,j) * Hz(i,j) - cmm(i,j) * (Ey(i,j) - Ey(i,j+1)) / dx;
        end
    end
    
    % 添加光源
    Ey(1,round(M/2)) = sin(w*t);
    
    % 绘图
    imagesc(y, x, Ey');
    caxis([-1, 1]);
    colormap(gray);
    xlabel('y');
    ylabel('x');
    title(['t = ' num2str(t)]);
    drawnow;
end

该代码实现了一个二维高阶FDTD算法，使用了Yee网格，采用了一阶Mur吸收边界条件。在计算区域中央添加了一个正弦波光源。运行代码可以得到随时间变化的电场分布图像。