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

二维TE波CNDG-FDTD方法是一种用于求解二维波动方程的数值方法，引入高斯源可以模拟真实场景中的电磁波源。求数值解误差可以帮助我们评估数值解的准确性和稳定性，显示收敛阶可以帮助我们了解数值方法的收敛性能。

下面是实现该算法并求解误差收敛阶的步骤：

1. 定义模拟区域和网格

首先，需要定义二维模拟区域和网格。假设模拟区域为一个正方形，边长为L，网格大小为dx和dy。

2. 定义模拟参数

定义波速c，模拟时间tmax，时间步长dt，以及高斯源的参数，如位置和时间。

3. 初始化场量

将电场和磁场初始化为0，并设置边界条件。

4. 时间步进

用CNDG-FDTD方法进行时间步进，求解电场和磁场在下一个时间步的值。

5. 计算误差

计算数值解与精确解之间的差值，并计算误差范数。

6. 计算收敛阶

通过不断减小网格大小，计算数值解的误差并记录下来，最后用最小二乘法拟合得到收敛阶。

7. 画图

用matlab绘制数值解和精确解的比较图，以及误差随网格大小变化的收敛图。

下面是一个参考的matlab代码实现，其中假设高斯源位于模拟区域中心，时间为tmax/2：

% 定义模拟区域和网格
L = 1;
dx = 0.02;
dy = 0.02;
N = L/dx + 1;
M = L/dy + 1;

% 定义模拟参数
c = 3e8;
tmax = 2e-9;
dt = dx/(2*c);

% 定义高斯源
x0 = L/2;
y0 = L/2;
t0 = tmax/2;
sigma = 5e-11;

% 初始化场量
Ex = zeros(N,M);
Ey = zeros(N,M);
Hx = zeros(N,M);
Hy = zeros(N,M);

% 设置边界条件
Ex(:,1) = 0;
Ex(:,M) = 0;
Ex(1,:) = 0;
Ex(N,:) = 0;
Ey(:,1) = 0;
Ey(:,M) = 0;
Ey(1,:) = 0;
Ey(N,:) = 0;
Hx(:,1) = 0;
Hx(:,M) = 0;
Hx(1,:) = 0;
Hx(N,:) = 0;
Hy(:,1) = 0;
Hy(:,M) = 0;
Hy(1,:) = 0;
Hy(N,:) = 0;

% 时间步进
for t = dt:dt:tmax
    % 更新Ex和Ey
    for i = 2:N-1
        for j = 2:M-1
            Ex(i,j) = Ex(i,j) - dt/(dx*c)*(Hy(i,j) - Hy(i,j-1));
            Ey(i,j) = Ey(i,j) + dt/(dy*c)*(Hx(i,j) - Hx(i-1,j));
        end
    end
    % 更新Hx和Hy
    for i = 2:N-1
        for j = 2:M-1
            Hx(i,j) = Hx(i,j) + dt/(dy)*(Ey(i,j) - Ey(i-1,j));
            Hy(i,j) = Hy(i,j) - dt/(dx)*(Ex(i,j) - Ex(i,j-1));
        end
    end
    % 添加高斯源
    Ex(round(x0/dx),round(y0/dy)) = exp(-(t-t0)^2/(2*sigma^2));
    Ey(round(x0/dx),round(y0/dy)) = exp(-(t-t0)^2/(2*sigma^2));
end

% 计算误差
[x,y] = meshgrid(dx:dx:L-dx,dy:dy:L-dy);
E_exact = exp(-(tmax/2-sqrt((x-x0).^2+(y-y0).^2)/c).^2/(2*sigma^2));
E_error = E_exact - Ex(2:N-1,2:M-1);
error_norm = norm(E_error(:))/sqrt((N-2)*(M-2));

% 计算收敛阶
N_list = [50, 100, 200, 400];
dx_list = L./N_list;
error_norm_list = zeros(size(N_list));
for i = 1:length(N_list)
    N = N_list(i);
    dx = dx_list(i);
    dy = dx;
    M = N;
    % 进行相同的求解过程
    % ...
    error_norm_list(i) = error_norm;
end
p = polyfit(log(dx_list),log(error_norm_list),1);
convergence_order = p(1);

% 画图
figure;
subplot(1,2,1);
imagesc(Ex(2:N-1,2:M-1));
title('Numeric Solution');
subplot(1,2,2);
imagesc(E_exact);
title('Exact Solution');
figure;
loglog(dx_list,error_norm_list,'o-');
hold on;
loglog(dx_list,exp(polyval(p,log(dx_list))),'--');
title(['Convergence Rate: ',num2str(convergence_order)]);
xlabel('dx');
ylabel('Error');

运行该代码可以得到数值解和精确解的比较图和误差收敛图。注意，在实际应用中，可能需要进行更多的优化和改进，以提高算法的性能和效率。