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

本文将介绍二维TE波CNDG-FDTD方法的matlab实现，其中将引入高斯源和PML边界条件，并计算数值解误差以显示收敛精度阶数。

首先，我们需要了解CNDG-FDTD方法的基本原理。该方法是通过将有限差分时间域（FDTD）方法与非协调网格（CNDG）方法相结合，从而实现高效且精确的电磁仿真。在二维TE波问题中，我们需要求解电磁场的横向电场分量Ez，其满足以下波动方程：

$$\frac{\partial^2E_z}{\partial x^2}+\frac{\partial^2E_z}{\partial y^2}+\omega^2\mu\epsilon E_z=0$$

其中，ω是角频率，μ和ε分别是介质的磁导率和电介质常数。

为了引入高斯源和PML边界条件，我们需要对上述方程进行一些修改。具体来说，我们可以将高斯源表示为：

$$E_z(x,y,t)=E_0exp(-\frac{(x-x_s)^2+(y-y_s)^2}{2\sigma^2})sin(\omega t)$$

其中，E0是振幅，σ是高斯源的展宽系数，xs和ys是高斯源的中心点坐标。

对于PML边界条件，我们可以通过在电场分量Ez中引入一个修正项来实现。具体来说，我们可以将Ez表示为：

$$E_z=E_{z0}exp(-\frac{\sigma_p}{\delta_p}x)$$

其中，Ez0是未经过修正的电场分量，σp和δp分别是PML层的展宽系数和厚度。

在进行数值模拟时，我们需要对空间进行离散化，将其划分为有限的网格。在本文中，我们将采用中心差分格式进行离散化，即：

$$\frac{\partial^2E_z}{\partial x^2}\approx\frac{E_{i+1,j}-2E_{i,j}+E_{i-1,j}}{\Delta x^2}$$

$$\frac{\partial^2E_z}{\partial y^2}\approx\frac{E_{i,j+1}-2E_{i,j}+E_{i,j-1}}{\Delta y^2}$$

其中，Δx和Δy分别是网格的横向和纵向尺寸。

通过将以上公式代入波动方程，我们可以得到以下差分方程：

$$\frac{E_{i+1,j}-2E_{i,j}+E_{i-1,j}}{\Delta x^2}+\frac{E_{i,j+1}-2E_{i,j}+E_{i,j-1}}{\Delta y^2}+\omega^2\mu\epsilon E_{i,j}=0$$

该方程可以通过时间步进法进行求解。在本文中，我们将采用Leapfrog法进行时间步进，即：

$$E_{i,j}^{n+1}=2E_{i,j}^{n}-E_{i,j}^{n-1}+\frac{c^2}{\omega^2\mu\epsilon}\left(\frac{E_{i+1,j}-2E_{i,j}+E_{i-1,j}}{\Delta x^2}+\frac{E_{i,j+1}-2E_{i,j}+E_{i,j-1}}{\Delta y^2}\right)$$

其中，c是光速。

现在我们已经得到了求解二维TE波问题的CNDG-FDTD方法的差分方程。接下来，我们可以使用matlab来实现该方法，并计算数值解误差以显示收敛精度阶数。

%% CNDG-FDTD method for 2D TE wave with Gaussian source and PML boundary conditions
clc;
clear all;
close all;

%% simulation parameters
c=3e8;                      % speed of light
mu=4*pi*1e-7;               % permeability of free space
epsilon=8.85e-12;           % permittivity of free space
omega=2*pi*10e9;            % angular frequency
lambda=2*pi*c/omega;        % wavelength
dx=lambda/20;               % grid size in x direction
dy=lambda/20;               % grid size in y direction
dt=dx/c/sqrt(2);            % time step
sigma_p=4;                  % PML width
delta_p=dx*sigma_p;         % PML thickness
sigma=2;                    % Gaussian source width
E0=1;                       % Gaussian source amplitude
xs=0;                       % Gaussian source center in x direction
ys=0;                       % Gaussian source center in y direction
nmax=1000;                  % maximum number of time steps
x_max=2*lambda;             % simulation region size in x direction
y_max=2*lambda;             % simulation region size in y direction
nx=ceil(x_max/dx);          % number of grid points in x direction
ny=ceil(y_max/dy);          % number of grid points in y direction

%% initialization
Ez=zeros(nx,ny);            % electric field along z direction
Ez_n=zeros(nx,ny);          % electric field at time step n
Ez_nm1=zeros(nx,ny);        % electric field at time step n-1
Hx=zeros(nx,ny);            % magnetic field along x direction
Hy=zeros(nx,ny);            % magnetic field along y direction
Hx_n=zeros(nx,ny);          % magnetic field at time step n
Hy_n=zeros(nx,ny);          % magnetic field at time step n
Jx=zeros(nx,ny);            % current density along x direction
Jy=zeros(nx,ny);            % current density along y direction
sigma_x=zeros(nx,ny);       % PML conductivity along x direction
sigma_y=zeros(nx,ny);       % PML conductivity along y direction
sigma_xy=zeros(nx,ny);      % PML conductivity along xy direction
Bx=zeros(nx+1,ny);          % auxiliary magnetic field along x direction
By=zeros(nx,ny+1);          % auxiliary magnetic field along y direction

%% PML conductivity
for i=1:nx
    for j=1:ny
        if i<=sigma_p
            sigma_x(i,j)=0.25*(sigma_p-i)^3/sigma_p^2;
        elseif i>=nx-sigma_p+1
            sigma_x(i,j)=0.25*(i-nx+sigma_p)^3/sigma_p^2;
        end
        if j<=sigma_p
            sigma_y(i,j)=0.25*(sigma_p-j)^3/sigma_p^2;
        elseif j>=ny-sigma_p+1
            sigma_y(i,j)=0.25*(j-ny+sigma_p)^3/sigma_p^2;
        end
        if i<=sigma_p && j<=sigma_p
            sigma_xy(i,j)=0.25*(sigma_p-sqrt((sigma_p-i)^2+(sigma_p-j)^2))^3/sigma_p^2;
        elseif i<=sigma_p && j>=ny-sigma_p+1
            sigma_xy(i,j)=0.25*(sigma_p-sqrt((sigma_p-i)^2+(j-ny+sigma_p)^2))^3/sigma_p^2;
        elseif i>=nx-sigma_p+1 && j<=sigma_p
            sigma_xy(i,j)=0.25*(sigma_p-sqrt((i-nx+sigma_p)^2+(sigma_p-j)^2))^3/sigma_p^2;
        elseif i>=nx-sigma_p+1 && j>=ny-sigma_p+1
            sigma_xy(i,j)=0.25*(sigma_p-sqrt((i-nx+sigma_p)^2+(j-ny+sigma_p)^2))^3/sigma_p^2;
        end
    end
end

%% main loop
for n=1:nmax
    % update magnetic field
    for i=1:nx
        for j=1:ny-1
            Hx(i,j)=Hx(i,j)+dt/(dy*mu)*(Ez(i,j)-Ez(i,j+1));
        end
    end
    for i=1:nx-1
        for j=1:ny
            Hy(i,j)=Hy(i,j)-dt/(dx*mu)*(Ez(i,j)-Ez(i+1,j));
        end
    end
    % PML boundary condition for magnetic field
    for i=1:nx
        for j=1:ny-1
            Bx(i,j)=Bx(i,j)+dt/(dy*mu)*(Ez(i,j)-Ez(i,j+1));
        end
    end
    for i=1:nx-1
        for j=1:ny
            By(i,j)=By(i,j)-dt/(dx*mu)*(Ez(i,j)-Ez(i+1,j));
        end
    end
    for i=1:sigma_p
        for j=1:ny-1
            Hx(i,j)=Hx(i,j)+dt*sigma_x(i,j)/(2*mu)*Bx(i,j);
        end
    end
    for i=nx-sigma_p+1:nx
        for j=1:ny-1
            Hx(i,j)=Hx(i,j)+dt*sigma_x(i,j)/(2*mu)*Bx(i,j);
        end
    end
    for i=1:nx-1
        for j=1:sigma_p
            Hy(i,j)=Hy(i,j)+dt*sigma_y(i,j)/(2*mu)*By(i,j);
        end
    end
    for i=1:nx-1
        for j=ny-sigma_p+1:ny
            Hy(i,j)=Hy(i,j)+dt*sigma_y(i,j)/(2*mu)*By(i,j);
        end
    end
    for i=1:sigma_p
        for j=1:sigma_p
            Hx(i,j)=Hx(i,j)+dt*sigma_xy(i,j)/(2*mu)*Bx(i,j);
            Hy(i,j)=Hy(i,j)+dt*sigma_xy(i,j)/(2*mu)*By(i,j);
        end
    end
    for i=nx-sigma_p+1:nx
        for j=1:sigma_p
            Hx(i,j)=Hx(i,j)+dt*sigma_xy(i,j)/(2*mu)*Bx(i,j);
            Hy(i,j)=Hy(i,j)+dt*sigma_xy(i,j)/(2*mu)*By(i,j);
        end
    end
    for i=1:sigma_p
        for j=ny-sigma_p+1:ny
            Hx(i,j)=Hx(i,j)+dt*sigma_xy(i,j)/(2*mu)*Bx(i,j);
            Hy(i,j)=Hy(i,j)+dt*sigma_xy(i,j)/(2*mu)*By(i,j);
        end
    end
    for i=nx-sigma_p+1:nx
        for j=ny-sigma_p+1:ny
            Hx(i,j)=Hx(i,j)+dt*sigma_xy(i,j)/(2*mu)*Bx(i,j);
            Hy(i,j)=Hy(i,j)+dt*sigma_xy(i,j)/(2*mu)*By(i,j);
        end
    end
    % update electric field
    for i=1:nx
        for j=2:ny-1
            Ez(i,j)=Ez(i,j)+dt/(dy*epsilon)*(Hy(i,j)-Hy(i,j-1));
        end
    end
    for i=2:nx-1
        for j=1:ny
            Ez(i,j)=Ez(i,j)-dt/(dx*epsilon)*(Hx(i,j)-Hx(i-1,j));
        end
    end
    % PML boundary condition for electric field
    for i=1:nx
        for j=2:ny-1
            Jx(i,j)=Jx(i,j)+dt*sigma_y(i,j)/epsilon*(Ez(i,j)-Ez(i,j-1));
        end
    end
    for i=2:nx-1
        for j=1:ny
            Jy(i,j)=Jy(i,j)-dt*sigma_x(i,j)/epsilon*(Ez(i,j)-Ez(i-1,j));
        end
    end
    for i=1:sigma_p
        for j=2:ny-1
            Ez(i,j)=Ez(i,j)+dt*sigma_y(i,j)/2/epsilon*(Hy(i,j)-Hy(i,j-1))-dt*Jx(i,j);
        end
    end
    for i=nx-sigma_p+1:nx
        for j=2:ny-1
            Ez(i,j)=Ez(i,j)+dt*sigma_y(i,j)/2/epsilon*(Hy(i,j)-Hy(i,j-1))-dt*Jx(i,j);
        end
    end
    for i=2:nx-1
        for j=1:sigma_p
            Ez(i,j)=Ez(i,j)-dt*sigma_x(i,j)/2/epsilon*(Hx(i,j)-Hx(i-1,j))+dt*Jy(i,j);
        end
    end
    for i=2:nx-1
        for j=ny-sigma_p+1:ny
            Ez(i,j)=Ez(i,j)-dt*sigma_x(i,j)/2/epsilon*(Hx(i,j)-Hx(i-1,j))+dt*Jy(i,j);
        end
    end
    for i=1:sigma_p
        for j=1:ny-1
            Ez(i,j)=Ez(i,j)+dt*sigma_xy(i,j)/2/epsilon*(Hy(i,j)-Hy(i,j+1))-dt*sigma_x(i,j)/2/epsilon*(Hx(i,j)-Hx(i-1,j))-dt*Jx(i,j);
        end
    end
    for i=nx-sigma_p+1:nx
        for j=1:ny-1
            Ez(i,j)=Ez(i,j)+dt*sigma_xy(i,j)/2/epsilon*(Hy(i,j)-Hy(i,j+1))-dt*sigma_x(i,j)/2/epsilon*(Hx(i,j)-Hx(i-1,j))-dt*Jx(i,j);
        end
    end
    for i=1:nx-1
        for j=1:sigma_p
            Ez(i,j)=Ez(i,j)-dt*sigma_xy(i,j)/2/epsilon*(Hx(i,j)-Hx(i+1,j))+dt*sigma_y(i,j)/2/epsilon*(Hy(i,j)-Hy(i,j-1))+dt*Jy(i,j);
        end
    end
    for i=1:nx-1
        for j=ny-sigma_p+1:ny
            Ez(i,j)=Ez(i,j)-dt*sigma_xy(i,j)/2/epsilon*(Hx(i,j)-Hx(i+1,j))+dt*sigma_y(i,j)/2/epsilon*(Hy(i,j)-Hy(i,j-1))+dt*Jy(i,j);
        end
    end
    for i=1:sigma_p
        for j=1:sigma_p
            Ez(i,j)=Ez(i,j)+dt*sigma_xy(i,j)/2/epsilon*(Hy(i,j)-Hy(i,j+1))-dt*sigma_x(i,j)/2/epsilon*(Hx(i,j)-Hx(i-1,j))-dt*sigma_y(i,j)/2/epsilon*(Ez(i,j)-Ez(i,j-1));
        end
    end
    for i=nx-sigma_p+1:nx
        for j=1:sigma_p
            Ez(i,j)=Ez(i,j)+dt*sigma_xy(i,j)/2/epsilon*(Hy(i,j)-Hy(i,j+1))-dt*sigma_x(i,j)/2/epsilon*(Hx(i,j)-Hx(i-1,j))-dt*sigma_y(i,j)/2/epsilon*(Ez(i,j)-Ez(i,j-1));
        end
    end
    for i=1:sigma_p
        for j=ny-sigma_p+1:ny
            Ez(i,j)=Ez(i,j)+dt*sigma_xy(i,j)/2/epsilon*(Hy(i,j)-Hy(i,j-1))-dt*sigma_x(i,j)/2/epsilon*(Hx(i,j)-Hx(i-1,j))-dt*sigma_y(i,j)/2/epsilon*(Ez(i,j)-Ez(i,j-1));
        end
    end
    for i=nx-sigma_p+1:nx
        for j=ny-sigma_p+1:ny
            Ez(i,j)=Ez(i,j)+dt*sigma_xy(i,j)/2/epsilon*(Hy(i,j)-Hy(i,j-1))-dt*sigma_x(i,j)/2/epsilon*(Hx(i,j)-Hx(i-1,j))-dt*sigma_y(i,j)/2/epsilon*(Ez(i,j)-Ez(i,j-1));
        end
    end
    % Gaussian source
    Ez(xs/dx+1,ys/dy+1)=E0*exp(-((n-1)*dt)^2/2/sigma^2)*sin(omega*(n-1)*dt);
    % boundary condition
    Ez(1,:)=0;
    Ez(nx,: