二维TE波CNDG-FDTD方法引入高斯源和PML边界条件，并求数值解误差，显示收敛精度阶数matlab实现

CNDG-FDTD方法是一种数值求解二维时域电磁波动方程的方法。在此基础上，引入高斯源和PML边界条件可以提高数值解的精度和稳定性。下面是matlab实现代码。

% 二维TE波CNDG-FDTD方法引入高斯源和PML边界条件
% 定义常数
c = 3e8;
mu0 = 4 * pi * 1e-7;
eps0 = 1 / (c^2 * mu0);
delta = 1e-3;
dt = delta / c;

% 定义网格和时间步长
Nx = 100;
Ny = 100;
tmax = 2e-9;
Nt = round(tmax / dt);

% 初始化场
Ex = zeros(Nx + 1, Ny);
Ey = zeros(Nx, Ny + 1);
Hz = zeros(Nx, Ny);

% 定义高斯源
xc = Nx / 2;
yc = Ny / 2;
sigma = 0.1e-9;
T0 = 0.5e-9;

% 定义PML参数
npml = 20;
kappa_max = 1;
alpha_max = 0.1;
sigma_max = -(kappa_max - 1) * eps0 * c / (delta * npml) * log(1e-6);

% 初始化PML参数
kappa_x = ones(Nx + 1, Ny);
kappa_y = ones(Nx, Ny + 1);
alpha_x = zeros(Nx + 1, Ny);
alpha_y = zeros(Nx, Ny + 1);
sigma_x = zeros(Nx + 1, Ny);
sigma_y = zeros(Nx, Ny + 1);

% 计算PML参数
for i = 1:npml
    x = delta * (i - 0.5);
    kappa_x(i, :) = kappa_max * (x / (delta * npml))^2 + (1 - (x / (delta * npml))^2);
    kappa_x(Nx + 2 - i, :) = kappa_max * (x / (delta * npml))^2 + (1 - (x / (delta * npml))^2);
    alpha_x(i, :) = alpha_max * ((npml - i + 0.5) / npml)^2;
    alpha_x(Nx + 2 - i, :) = alpha_max * ((npml - i + 0.5) / npml)^2;
    sigma_x(i, :) = sigma_max * ((npml - i + 0.5) / npml)^3;
    sigma_x(Nx + 2 - i, :) = sigma_max * ((npml - i + 0.5) / npml)^3;
end
for j = 1:npml
    y = delta * (j - 0.5);
    kappa_y(:, j) = kappa_max * (y / (delta * npml))^2 + (1 - (y / (delta * npml))^2);
    kappa_y(:, Ny + 2 - j) = kappa_max * (y / (delta * npml))^2 + (1 - (y / (delta * npml))^2);
    alpha_y(:, j) = alpha_max * ((npml - j + 0.5) / npml)^2;
    alpha_y(:, Ny + 2 - j) = alpha_max * ((npml - j + 0.5) / npml)^2;
    sigma_y(:, j) = sigma_max * ((npml - j + 0.5) / npml)^3;
    sigma_y(:, Ny + 2 - j) = sigma_max * ((npml - j + 0.5) / npml)^3;
end

% 时间步进
for n = 1:Nt
    % 更新Ex
    Ex(:, 2:Ny) = kappa_x(:, 2:Ny) .* Ex(:, 2:Ny) - alpha_x(:, 2:Ny) .* (Hz(:, 2:Ny) - Hz(:, 1:Ny-1)) - sigma_x(:, 2:Ny) .* Ex(:, 2:Ny) + kappa_x(:, 2:Ny) .* (Ey(:, 2:Ny) - Ey(:, 1:Ny-1)) * dt / eps0;
    % 添加高斯源
    Ex(xc, yc) = exp(-0.5 * ((n * dt - T0) / sigma)^2);
    % 更新Ey
    Ey(2:Nx, :) = kappa_y(2:Nx, :) .* Ey(2:Nx, :) - alpha_y(2:Nx, :) .* (Hz(2:Nx, :) - Hz(1:Nx-1, :)) + sigma_y(2:Nx, :) .* Ex(2:Nx, 1:Ny-1) * dt / eps0;
    % 更新Hz
    Hz(1:Nx-1, 1:Ny-1) = Hz(1:Nx-1, 1:Ny-1) - mu0 * (Ex(2:Nx, 1:Ny-1) - Ex(1:Nx-1, 1:Ny-1) + Ey(1:Nx-1, 2:Ny) - Ey(1:Nx-1, 1:Ny-1)) * dt / delta;
    % 更新PML边界
    Hz(1:npml, :) = Hz(1:npml, :) - mu0 * Ex(2:npml+1, :) * dt / delta;
    Hz(Nx-npml+1:Nx-1, :) = Hz(Nx-npml+1:Nx-1, :) - mu0 * Ex(Nx-npml:Nx-2, :) * dt / delta;
    Hz(:, 1:npml) = Hz(:, 1:npml) + mu0 * Ey(:, 2:npml+1) * dt / delta;
    Hz(:, Ny-npml+1:Ny-1) = Hz(:, Ny-npml+1:Ny-1) + mu0 * Ey(:, Ny-npml:Ny-2) * dt / delta;
end

% 计算解析解
[x, y] = meshgrid(0:delta:Nx*delta, 0:delta:Ny*delta);
t = Nt * dt;
E0 = 1;
Er = exp(-0.5 * ((t - T0) / sigma)^2);
Ea = E0 * Er * sin(pi * x / (Nx*delta)) .* sin(pi * y / (Ny*delta));

% 计算误差和收敛精度阶数
Ee = sqrt((Ex(:, 1:Ny-1) - Ex(:, 2:Ny)).^2 + (Ey(1:Nx-1, :) - Ey(2:Nx, :)).^2);
Ea = Ea(1:Nx, 1:Ny-1);
err = norm(Ee - Ea, 2) / norm(Ea, 2);
p = log(err(2:end) ./ err(1:end-1)) / log(2);

% 显示结果
figure;
subplot(1, 2, 1);
imagesc(0:delta:Nx*delta, 0:delta:(Ny-1)*delta, Ee);
title('数值解');
xlabel('x (m)');
ylabel('y (m)');
colorbar;
subplot(1, 2, 2);
imagesc(0:delta:Nx*delta, 0:delta:(Ny-1)*delta, Ea);
title('解析解');
xlabel('x (m)');
ylabel('y (m)');
colorbar;
figure;
plot(1:length(p), p);
title('收敛精度阶数');
xlabel('网格数');
ylabel('阶数');

代码中，我们首先定义了常数和时间步长，然后初始化场和PML参数。在时间步进的过程中，我们先更新Ex，然后添加高斯源，接着更新Ey和Hz。最后，我们更新PML边界。在计算完数值解后，我们计算解析解并计算误差和收敛精度阶数。最后，我们显示数值解、解析解和收敛精度阶数。