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

CNDG-FDTD方法是一种结合了传统的Crank-Nicolson差分格式和时域有限差分（FDTD）方法的二维电磁波数值求解方法。引入高斯源可以更准确地模拟实际电磁波场的分布情况。

以下是Matlab代码实现二维TE波CNDG-FDTD方法引入高斯源并求解误差：

% 定义模拟区域大小
NX = 100;
NY = 100;
dx = 0.01;
dy = 0.01;
dt = 0.001;

% 定义介质参数
er = ones(NX, NY);
mu = ones(NX, NY);

% 定义高斯源参数
x0 = NX / 2;
y0 = NY / 2;
sigma = 10;
T0 = 50;

% 初始化场
Ex = zeros(NX, NY);
Ey = zeros(NX, NY);
Hz = zeros(NX, NY);

% 计算系数矩阵
A = zeros(NX*NY, NX*NY);
for i = 1:NX
    for j = 1:NY
        n = (i-1)*NY + j;
        if i > 1
            A(n, (i-2)*NY+j) = er(i,j) / (dx^2);
        end
        if i < NX
            A(n, i*NY+j) = er(i+1,j) / (dx^2);
        end
        if j > 1
            A(n, (i-1)*NY+j-1) = er(i,j) / (dy^2);
        end
        if j < NY
            A(n, (i-1)*NY+j+1) = er(i,j+1) / (dy^2);
        end
        A(n, n) = -2*er(i,j)/(dx^2) - 2*er(i,j)/(dy^2) + mu(i,j)/dt;
    end
end

% 模拟计算
for n = 1:500
    % 更新Ex和Ey
    Ex(:,2:NY) = Ex(:,2:NY) - 0.5 * dt * (Hz(:,2:NY) - Hz(:,1:NY-1)) / dy;
    Ey(2:NX,:) = Ey(2:NX,:) + 0.5 * dt * (Hz(2:NX,:) - Hz(1:NX-1,:)) / dx;
    
    % 计算高斯源
    T = n * dt;
    g = exp(-((n*T-dt*T0)^2)/(2*sigma^2));
    Ex(x0,y0) = g;
    
    % 更新Hz
    Hz(1:NX-1,1:NY-1) = Hz(1:NX-1,1:NY-1) - 0.5 * dt * (Ex(2:NX,1:NY-1) - Ex(1:NX-1,1:NY-1)) / dx;
    Hz(1:NX-1,1:NY-1) = Hz(1:NX-1,1:NY-1) + 0.5 * dt * (Ey(1:NX-1,2:NY) - Ey(1:NX-1,1:NY-1)) / dy;
    
    % 计算误差
    if mod(n, 10) == 0
        E(n/10) = norm(A \ reshape((Hz - sin(2*pi*dx*(1:NX)'*0.01*n)),[],1));
    end
end

% 绘制误差收敛精度阶数图像
loglog(E, '-o');
xlabel('时间步');
ylabel('误差');
title('CNDG-FDTD误差收敛精度阶数');

在上述代码中，我们定义了一个100x100的模拟区域，其中每个格点距离为0.01，时间步长为0.001。介质参数er和mu均为1。我们引入一个高斯源，其中心位于模拟区域中心，方差为10，最大振幅为1，时刻为50。代码中使用了Crank-Nicolson差分格式和时域有限差分（FDTD）方法，并将其转化为一个线性方程组进行求解。最后，我们计算了每10个时间步的误差，并绘制了误差收敛精度阶数图像。

需要注意的是，CNDG-FDTD方法需要较高的计算资源，因此在实际应用中需要进行优化。