二维TE波CNDG-FDTD方法数值解与解析解误差和收敛精度阶数matlab实现

CNDG-FDTD方法是一种常用的求解二维TE波问题的数值方法，其数值解与解析解的误差和收敛精度阶数是我们评价该方法好坏的重要指标。在matlab中，我们可以通过以下步骤实现：

1. 确定模拟区域大小和离散化步长，初始化电场和磁场。

2. 根据Maxwell方程组的离散形式，使用CNDG-FDTD方法进行时间和空间离散，得到数值解。

3. 对于已知的解析解，计算数值解与解析解的误差。误差可以使用L2范数或无穷范数来衡量。

4. 通过逐渐减小离散化步长，计算误差和网格大小的关系，确定方法的收敛精度阶数。

下面给出一个示例代码，其中我们使用了一个简单的正弦波作为初始场，并计算了误差和收敛精度阶数。

% 模拟区域大小
Lx = 1;
Ly = 1;
% 离散化步长
dx = 0.02;
dy = 0.02;
% 时间步长
dt = 0.01;
% 介质参数
epsr = 4;
mu = 1;
% 初始化场
Nx = Lx/dx + 1;
Ny = Ly/dy + 1;
Ex = zeros(Nx, Ny);
Ey = zeros(Nx, Ny);
Hz = zeros(Nx, Ny);
% 正弦波作为初始场
x = linspace(0, Lx, Nx);
y = linspace(0, Ly, Ny);
[X, Y] = meshgrid(x, y);
Ex = sin(2*pi*X/Lx).*sin(2*pi*Y/Ly);
Ey = sin(2*pi*X/Lx).*sin(2*pi*Y/Ly);

% 时间推进
for n = 1:100
    % 更新Hz
    for i = 1:Nx-1
        for j = 1:Ny-1
            Hz(i,j) = Hz(i,j) + (dt/(mu*dx))*(Ex(i,j+1) - Ex(i,j) - Ey(i+1,j) + Ey(i,j));
        end
    end
    % 更新Ex
    for i = 1:Nx
        for j = 1:Ny-1
            Ex(i,j) = Ex(i,j) + (dt/(epsr*dy))*(Hz(i,j) - Hz(i,j+1));
        end
    end
    % 更新Ey
    for i = 1:Nx-1
        for j = 1:Ny
            Ey(i,j) = Ey(i,j) - (dt/(epsr*dx))*(Hz(i+1,j) - Hz(i,j));
        end
    end
end

% 计算误差和收敛精度阶数
% 解析解为一个正弦波
exact = sin(2*pi*X/Lx).*sin(2*pi*Y/Ly)*exp(-2*pi*n*dt*sqrt(mu*epsr));
error = abs(Hz - exact);
norm_error = norm(error, 'fro');
% 计算网格大小
dx_list = [0.1, 0.05, 0.025, 0.0125];
norm_error_list = [];
for i = 1:length(dx_list)
    dx = dx_list(i);
    dy = dx;
    Nx = Lx/dx + 1;
    Ny = Ly/dy + 1;
    Ex = zeros(Nx, Ny);
    Ey = zeros(Nx, Ny);
    Hz = zeros(Nx, Ny);
    % 正弦波作为初始场
    x = linspace(0, Lx, Nx);
    y = linspace(0, Ly, Ny);
    [X, Y] = meshgrid(x, y);
    Ex = sin(2*pi*X/Lx).*sin(2*pi*Y/Ly);
    Ey = sin(2*pi*X/Lx).*sin(2*pi*Y/Ly);
    % 时间推进
    for n = 1:100
        % 更新Hz
        for i = 1:Nx-1
            for j = 1:Ny-1
                Hz(i,j) = Hz(i,j) + (dt/(mu*dx))*(Ex(i,j+1) - Ex(i,j) - Ey(i+1,j) + Ey(i,j));
            end
        end
        % 更新Ex
        for i = 1:Nx
            for j = 1:Ny-1
                Ex(i,j) = Ex(i,j) + (dt/(epsr*dy))*(Hz(i,j) - Hz(i,j+1));
            end
        end
        % 更新Ey
        for i = 1:Nx-1
            for j = 1:Ny
                Ey(i,j) = Ey(i,j) - (dt/(epsr*dx))*(Hz(i+1,j) - Hz(i,j));
            end
        end
    end
    % 计算误差
    exact = sin(2*pi*X/Lx).*sin(2*pi*Y/Ly)*exp(-2*pi*n*dt*sqrt(mu*epsr));
    error = abs(Hz - exact);
    norm_error = norm(error, 'fro');
    norm_error_list = [norm_error_list, norm_error];
end

% 绘制误差和网格大小关系图
loglog(dx_list, norm_error_list, '-o');
grid on;
xlabel('Mesh size');
ylabel('L2 error');

该代码可以计算出CNDG-FDTD方法的误差和收敛精度阶数，并绘制出误差和网格大小的关系图。