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

CNDG-FDTD方法是一种结合了传统的有限差分时域（FDTD）方法和Curl-Nodal Discontinuous Galerkin（CNDG）方法的求解电磁波问题的方法。在二维情况下，TE波表示电场只在z方向有分量，磁场只在x、y方向有分量。在本文中，我们将引入高斯源，并利用Matlab实现求解数值解与解析解误差、显示收敛阶并画图。

假设我们要求解的是一个长宽分别为Lx和Ly的矩形区域内的TE波。首先，我们需要将区域网格化，将整个区域划分为N_x \times N_y个小矩形，每个小矩形的长宽分别为\Delta x和\Delta y。我们假设在该区域内引入一个高斯源，其表达式为：

$$
S(x,y,t)=S_0 exp\left(-\frac{(x-x_0)^2+(y-y_0)^2}{2\sigma^2}\right)cos(\omega t)
$$

其中，S_0是源的振幅，(x_0,y_0)是源的位置，\sigma是源的宽度，\omega是源的频率。

接下来，我们可以利用CNDG-FDTD方法求解TE波方程：

$$
\nabla \times \nabla \times \mathbf{E} - \mu \epsilon \frac{\partial^2 \mathbf{E}}{\partial t^2} = \mu \mathbf{J}_s
$$

其中，\mathbf{E}表示电场，\mu和\epsilon分别表示介质的磁导率和电介质常数，\mathbf{J}_s表示高斯源产生的电流密度。我们可以将该方程离散化为：

$$
\mathbf{C} \frac{\partial \mathbf{E}}{\partial t} + \mathbf{G} \mathbf{E} = \mathbf{S}
$$

其中，\mathbf{C}和\mathbf{G}是系数矩阵，\mathbf{S}是源项矩阵。我们可以利用FDTD方法求解\mathbf{C}和\mathbf{G}，利用CNDG方法求解\mathbf{S}。

最后，我们可以利用解析解求出TE波的精确解，与数值解做比较，计算误差，并绘制误差随网格大小变化的图像，以显示收敛阶。具体实现可参考以下Matlab代码：

% 定义常数
Lx = 1; % 区域长度
Ly = 1; % 区域宽度
Nx = 100; % x方向网格数
Ny = 100; % y方向网格数
dx = Lx / Nx; % x方向网格间距
dy = Ly / Ny; % y方向网格间距
T = 1e-8; % 模拟总时间
dt = dx / sqrt(2); % 时间间隔
Nt = T / dt; % 时间步数
x = linspace(0, Lx, Nx+1); % x方向网格
y = linspace(0, Ly, Ny+1); % y方向网格
omega = 2 * pi * 1e9; % 频率
mu = pi * 4e-7; % 磁导率
epsilon = 8.85e-12; % 电介质常数
sigma = dx / 4; % 高斯源宽度
S0 = 1; % 高斯源振幅
x0 = Lx / 2; % 高斯源x方向位置
y0 = Ly / 2; % 高斯源y方向位置

% 初始化系数矩阵和源项矩阵
C = zeros(Nx*Ny, Nx*Ny);
G = zeros(Nx*Ny, Nx*Ny);
S = zeros(Nx*Ny, 1);

% 构造系数矩阵和源项矩阵
for i = 1:Nx
    for j = 1:Ny
        % 计算节点编号
        n = (j-1)*Nx + i;

        % 构造C矩阵
        if i == 1 || i == Nx || j == 1 || j == Ny
            C(n, n) = 1;
        else
            C(n, (j-2)*Nx+i) = 1/(dy*dy);
            C(n, j*Nx+i) = 1/(dy*dy);
            C(n, (j-1)*Nx+i-1) = 1/(dx*dx);
            C(n, (j-1)*Nx+i+1) = 1/(dx*dx);
            C(n, n) = -2/(dx*dx) - 2/(dy*dy);
        end

        % 构造G矩阵
        G(n, n) = mu * epsilon / (dt*dt);

        % 构造S矩阵
        xj = (i-0.5) * dx;
        yj = (j-0.5) * dy;
        r = sqrt((xj-x0)^2 + (yj-y0)^2);
        S(n) = S0 * exp(-r^2/(2*sigma*sigma)) * cos(omega*dt*Nt);
    end
end

% 求解电场
E = zeros(Nx*Ny, Nt);
E(:, 1) = 0;
E(:, 2) = dt * inv(C) * (S - G * E(:, 1));

for n = 3:Nt
    E(:, n) = inv(C) * (2*S - G * E(:, n-1) - E(:, n-2));
end

% 计算解析解
E_analytical = zeros(Nx, Ny, Nt);
for n = 1:Nt
    for i = 1:Nx
        for j = 1:Ny
            xj = (i-0.5) * dx;
            yj = (j-0.5) * dy;
            r = sqrt((xj-x0)^2 + (yj-y0)^2);
            E_analytical(i, j, n) = S0 * exp(-r^2/(2*sigma*sigma)) * cos(omega*dt*n);
        end
    end
end

% 计算误差
error = zeros(1, Nt);
for n = 1:Nt
    error(n) = norm(E(:, n) - reshape(E_analytical, Nx*Ny, Nt)(:, n));
end

% 计算收敛阶
order = log(error(2:end)./error(1:end-1))/log(sqrt(dx*dx+dy*dy));
order = [order(1) order];

% 绘制误差随时间变化图像
figure;
plot(dt*(0:Nt-1), error);
xlabel('时间（秒）');
ylabel('误差');

% 绘制误差随网格大小变化图像
figure;
loglog(sqrt(dx*dx+dy*dy).^(1:Nt), error);
xlabel('网格大小');
ylabel('误差');

% 绘制收敛阶随时间变化图像
figure;
plot(dt*(1:Nt), order);
xlabel('时间（秒）');
ylabel('收敛阶');

需要注意的是，本文中的TE波方程是在时域中求解的，因此需要将频率离散化为角频率，并在源项中加入cos(\omega t)项。此外，我们还需要注意在计算解析解时的节点位置和电场方向。