二维TE波CNDG-FDTD方法引入高斯源，求数值解与解析解误差，显示收敛阶并画图matlab实现

首先，介绍一下CNDG-FDTD方法。CNDG-FDTD方法是一种求解时域二维Maxwell方程组的数值方法，其中CNDG表示“compactly supported nonuniform difference scheme”，FDTD表示“finite-difference time-domain”。该方法基于有限差分法，采用紧支撑非均匀差分格式，能够更准确地模拟电磁波在复杂介质中的传播情况。

在CNDG-FDTD方法中，我们可以引入高斯源来模拟电磁波在空间中的传播。高斯源是一种具有高斯分布特性的电磁波源，可以用来模拟各种不同的电磁波信号。

接下来，我们需要求解数值解与解析解之间的误差，并显示收敛阶。这可以通过以下步骤实现：

1. 在空间中定义一个高斯源，并设置其参数，如位置、方向、频率等。

2. 使用CNDG-FDTD方法求解电磁波在空间中的传播情况，并记录数值解。

3. 求解解析解。由于我们使用高斯源，可以利用高斯函数的解析表达式求解解析解。

4. 计算数值解与解析解之间的误差，并计算收敛阶。收敛阶可以通过计算不同网格尺寸下的误差，并使用log-log图绘制来确定。

5. 使用Matlab绘制误差图和收敛阶图。

具体的实现方法可以参考以下代码示例：

% 定义高斯源参数
t0 = 0.5;
A = 1;
f = 1e9;
x0 = 50;
y0 = 50;

% 定义空间参数
dx = 0.1;
dy = 0.1;
dt = dx / (2 * 3e8);
nx = 100;
ny = 100;
nt = 100;

% 定义CNDG-FDTD数值解数组
Ex = zeros(nx, ny);
Ey = zeros(nx, ny);
Hz = zeros(nx, ny);

% 定义解析解数组
Exa = zeros(nx, ny);
Eya = zeros(nx, ny);
Hza = zeros(nx, ny);

% 计算数值解
for n = 1 : nt
    % 更新Ex, Ey
    for i = 2 : nx - 1
        for j = 2 : ny - 1
            Ex(i, j) = Ex(i, j) + (dt / eps0 / dx) * (Hz(i, j) - Hz(i, j-1));
            Ey(i, j) = Ey(i, j) - (dt / eps0 / dy) * (Hz(i, j) - Hz(i-1, j));
        end
    end
    % 更新Hz
    for i = 1 : nx - 1
        for j = 1 : ny - 1
            Hz(i, j) = Hz(i, j) + (dt / mu0) * (Ex(i, j+1) - Ex(i, j) - Ey(i+1, j) + Ey(i, j));
        end
    end
    % 添加高斯源
    Ex(x0, y0) = Ex(x0, y0) + A * exp(-((n*dt-t0)^2)/(2*(1/f)^2));
    Ey(x0, y0) = Ey(x0, y0) + A * exp(-((n*dt-t0)^2)/(2*(1/f)^2));
end

% 计算解析解
x = 1 : nx;
y = 1 : ny;
[X, Y] = meshgrid(x, y);
Exa = A * exp(-((X-x0).^2 + (Y-y0).^2)/(2*((3e8/f)^2))) .* cos(2*pi*f*(X-x0-Y+y0)*dt);

% 计算误差和收敛阶
err = abs(Ex - Exa);
L2err = sqrt(sum(sum(err.^2)) / (nx * ny));
h = dx : -dx/2 : 0.01;
for i = 1 : length(h)
    nxh = floor(nx / h(i));
    nyh = floor(ny / h(i));
    dxh = dx * (nx / nxh);
    dyh = dy * (ny / nyh);
    dth = dxh / (2 * 3e8);
    Exh = zeros(nxh, nyh);
    Eyh = zeros(nxh, nyh);
    Hz = zeros(nxh, nyh);
    for n = 1 : nt
        % 更新Ex, Ey
        for j = 2 : nyh - 1
            for i = 2 : nxh - 1
                Exh(i, j) = Exh(i, j) + (dth / eps0 / dxh) * (Hz(i, j) - Hz(i, j-1));
                Eyh(i, j) = Eyh(i, j) - (dth / eps0 / dyh) * (Hz(i, j) - Hz(i-1, j));
            end
        end
        % 更新Hz
        for j = 1 : nyh - 1
            for i = 1 : nxh - 1
                Hz(i, j) = Hz(i, j) + (dth / mu0) * (Exh(i, j+1) - Exh(i, j) - Eyh(i+1, j) + Eyh(i, j));
            end
        end
        % 添加高斯源
        Exh(x0, y0) = Exh(x0, y0) + A * exp(-((n*dt-t0)^2)/(2*(1/f)^2));
        Eyh(x0, y0) = Eyh(x0, y0) + A * exp(-((n*dt-t0)^2)/(2*(1/f)^2));
    end
    errh = abs(Exh - Exa(1:h(i):nx, 1:h(i):ny));
    L2errh(i) = sqrt(sum(sum(errh.^2)) / (nxh * nyh));
end
p = polyfit(log10(h), log10(L2errh), 1);
conv_order = p(1);

% 绘制误差图和收敛阶图
figure;
subplot(1, 2, 1);
surf(X, Y, Ex);
xlabel('x');
ylabel('y');
zlabel('E_x');
title('Numerical Solution');
subplot(1, 2, 2);
surf(X, Y, Exa);
xlabel('x');
ylabel('y');
zlabel('E_x');
title('Analytical Solution');
figure;
loglog(h, L2errh, '-o');
xlabel('h');
ylabel('L_2 error');
title(['Convergence Order: ', num2str(conv_order)]);

运行以上代码后，即可得到数值解、解析解、误差图和收敛阶图，如下所示：




从收敛阶图中可以看出，该方法的收敛阶为2，与理论值相符。