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

二维TE波（横向电场）的CNDG-FDTD方法数值解与解析解误差和显示收敛精度阶数的Matlab实现如下：

首先，我们需要定义模拟区域的大小和时间步长，以及波的频率和传播速度等参数。假设模拟区域大小为 $L_x \times L_y$，时间步长为 $\Delta t$，波的频率为 $\omega$，传播速度为 $c$，则可以定义如下变量：

Lx = 1; % 模拟区域长度
Ly = 1; % 模拟区域宽度
dx = 0.01; % 空间步长
dy = 0.01;
dt = 0.0001; % 时间步长
omega = 2*pi*1e9; % 波的频率
c = 3e8; % 传播速度

接下来，我们需要初始化场变量和介质参数。假设介质是均匀的，并且在模拟区域的中心有一个点源产生波，可以定义如下变量：

N = round(Lx/dx)+1; % 空间网格数
M = round(Ly/dy)+1;
Ez = zeros(N,M); % 横向电场
epsilon = ones(N,M); % 介质电容率
sigma = zeros(N,M); % 介质电导率
source_x = round(N/2); % 波源位置
source_y = round(M/2);

接下来，就是实现 CNDG-FDTD 方法的核心代码了。这里我们采用了较为简单的一阶中心差分格式，实现起来比较容易。具体实现过程如下：

% 一阶中心差分格式，计算 Ez 的时间导数
dEzdt = zeros(N,M);
dEzdt(2:end-1,2:end-1) = (Ez(3:end,2:end-1)-Ez(1:end-2,2:end-1))/(2*dx) + (Ez(2:end-1,3:end)-Ez(2:end-1,1:end-2))/(2*dy);
% 更新 Ez
Ez = Ez + dt*c^2./epsilon.*dEzdt;
% 在波源位置施加脉冲
Ez(source_x,source_y) = sin(omega*dt);

接下来，我们可以计算出数值解和解析解，并计算它们之间的误差。这里我们采用了解析解为平面波的形式，具体代码如下：

% 计算解析解
x = linspace(0,Lx,N);
y = linspace(0,Ly,M);
[X,Y] = meshgrid(x,y);
r = sqrt((X-source_x*dx).^2 + (Y-source_y*dy).^2);
analytical = sin(omega*dt - 2*pi*r/c) ./ r;
analytical(isnan(analytical)) = 0;
% 计算误差
error = abs(Ez-analytical);

最后，我们还可以计算出 CNDG-FDTD 方法的显示收敛精度阶数。这里我们可以采用不同的网格数，分别计算出数值解和解析解的误差，然后通过线性回归计算出收敛精度阶数。具体代码如下：

% 计算不同网格下的误差
dx_values = logspace(-3,-1,20);
error_values = zeros(size(dx_values));
for i=1:length(dx_values)
    dx = dx_values(i);
    dy = dx;
    N = round(Lx/dx)+1;
    M = round(Ly/dy)+1;
    x = linspace(0,Lx,N);
    y = linspace(0,Ly,M);
    [X,Y] = meshgrid(x,y);
    epsilon = ones(N,M);
    sigma = zeros(N,M);
    Ez = zeros(N,M);
    source_x = round(N/2);
    source_y = round(M/2);
    dEzdt = zeros(N,M);
    analytical = zeros(N,M);
    for n=1:1000
        dEzdt(2:end-1,2:end-1) = (Ez(3:end,2:end-1)-Ez(1:end-2,2:end-1))/(2*dx) + (Ez(2:end-1,3:end)-Ez(2:end-1,1:end-2))/(2*dy);
        Ez = Ez + dt*c^2./epsilon.*dEzdt;
        Ez(source_x,source_y) = sin(omega*n*dt);
        r = sqrt((X-source_x*dx).^2 + (Y-source_y*dy).^2);
        analytical = sin(omega*n*dt - 2*pi*r/c) ./ r;
        analytical(isnan(analytical)) = 0;
        error_values(i) = norm(Ez-analytical,inf);
    end
end
% 计算收敛精度阶数
p = polyfit(log(dx_values),log(error_values),1);
convergence_order = -p(1);

以上就是二维TE波CNDG-FDTD方法数值解与解析解误差和显示收敛精度阶数的Matlab实现。