二维TE波ADI-FDTD方法，引入PML吸收边界条件，采用紧差分格式matlab实现

二维TE波ADI-FDTD方法是一种用于求解电磁波传播的数值方法，ADI-FDTD方法结合了时间分步法和空间差分法的优点，具有较高的数值稳定性和计算效率。PML（吸收边界条件）是一种可以在计算区域外部完全吸收电磁波的技术，常用于防止计算结果出现反射和波漏等现象。

以下是使用matlab实现二维TE波ADI-FDTD方法，引入PML吸收边界条件的代码：

%定义计算区域的大小和网格大小
Lx = 1; %计算区域的长度
Ly = 1; %计算区域的宽度
dx = 0.01; %x方向上的网格大小
dy = 0.01; %y方向上的网格大小
Nx = round(Lx/dx); %x方向上的网格数
Ny = round(Ly/dy); %y方向上的网格数

%定义时间和空间步长
c = 3e8; %电磁波速度
f0 = 10e9; %中心频率
lambda = c/f0; %波长
dt = dx/(2*c); %时间步长

%定义电场和磁场
Ex = zeros(Nx,Ny); %x方向电场
Ey = zeros(Nx,Ny); %y方向电场
Hx = zeros(Nx,Ny); %x方向磁场
Hy = zeros(Nx,Ny); %y方向磁场

%定义PML吸收边界的参数
npml = 10; %PML层数
sigma_max = 3*(npml+1)*(-log(1e-6)/(2*lambda*(npml*dx)^2)); %最大吸收率
kappa_max = 1; %最大耗散率
sigma_x = zeros(Nx,Ny); %x方向吸收率
sigma_y = zeros(Nx,Ny); %y方向吸收率
kappa_x = zeros(Nx,Ny); %x方向耗散率
kappa_y = zeros(Nx,Ny); %y方向耗散率
for i = 1:npml
    sigma_x(i,:) = sigma_max*((npml-i+1)/npml)^3;
    sigma_x(Nx-i+1,:) = sigma_max*((npml-i+1)/npml)^3;
    kappa_x(i,:) = 1+(kappa_max-1)*((npml-i+1)/npml)^3;
    kappa_x(Nx-i+1,:) = 1+(kappa_max-1)*((npml-i+1)/npml)^3;
    sigma_y(:,i) = sigma_max*((npml-i+1)/npml)^3;
    sigma_y(:,Ny-i+1) = sigma_max*((npml-i+1)/npml)^3;
    kappa_y(:,i) = 1+(kappa_max-1)*((npml-i+1)/npml)^3;
    kappa_y(:,Ny-i+1) = 1+(kappa_max-1)*((npml-i+1)/npml)^3;
end

%定义计算系数
mEy = c*dt/dy; %Ey系数
mEx = c*dt/dx; %Ex系数
mHx = dt/(mu*dx); %Hx系数
mHy = dt/(mu*dy); %Hy系数

%开始时间步进
for n = 1:Nt
    %更新Hx和Hy
    for i = 1:Nx-1
        for j = 1:Ny-1
            Hx(i,j) = Hx(i,j) - mHy*(Ey(i,j+1)-Ey(i,j)) + sigma_y(i,j)*Hx(i,j);
            Hy(i,j) = Hy(i,j) + mHx*(Ex(i+1,j)-Ex(i,j)) - sigma_x(i,j)*Hy(i,j);
        end
    end
    
    %在PML区域内更新Hx和Hy
    for i = 1:Nx
        for j = 1:npml
            Hx(i,j) = Hx(i,j) - mHy*(Ey(i,j+1)-Ey(i,j)) + sigma_y(i,j)*Hx(i,j);
            Hy(i,j) = Hy(i,j) + mHx*(Ex(i+1,j)-Ex(i,j)) - sigma_x(i,j)*Hy(i,j);
            Hx(i,Ny-j+1) = Hx(i,Ny-j+1) - mHy*(Ey(i,Ny-j+1)-Ey(i,Ny-j)) + sigma_y(i,Ny-j+1)*Hx(i,Ny-j+1);
            Hy(i,Ny-j+1) = Hy(i,Ny-j+1) + mHx*(Ex(i+1,Ny-j+1)-Ex(i,Ny-j+1)) - sigma_x(i,Ny-j+1)*Hy(i,Ny-j+1);
        end
    end
    for j = 1:Ny
        for i = 1:npml
            Hx(i,j) = Hx(i,j) - mHy*(Ey(i,j+1)-Ey(i,j)) + sigma_y(i,j)*Hx(i,j);
            Hy(i,j) = Hy(i,j) + mHx*(Ex(i+1,j)-Ex(i,j)) - sigma_x(i,j)*Hy(i,j);
            Hx(Nx-i+1,j) = Hx(Nx-i+1,j) - mHy*(Ey(Nx-i+1,j+1)-Ey(Nx-i+1,j)) + sigma_y(Nx-i+1,j)*Hx(Nx-i+1,j);
            Hy(Nx-i+1,j) = Hy(Nx-i+1,j) + mHx*(Ex(Nx-i+2,j)-Ex(Nx-i+1,j)) - sigma_x(Nx-i+1,j)*Hy(Nx-i+1,j);
        end
    end
    
    %更新Ex和Ey
    for j = 2:Ny
        for i = 2:Nx
            Ex(i,j) = Ex(i,j) + mEy*(Hy(i,j)-Hy(i-1,j)) - sigma_y(i,j)*Ex(i,j);
            Ey(i,j) = Ey(i,j) - mEx*(Hx(i,j)-Hx(i,j-1)) + sigma_x(i,j)*Ey(i,j);
        end
    end
    
    %在PML区域内更新Ex和Ey
    for j = 1:Ny
        for i = 1:npml
            Ex(i,j) = Ex(i,j) + mEy*(Hy(i,j)-Hy(i-1,j)) - sigma_y(i,j)*Ex(i,j);
            Ey(i,j) = Ey(i,j) - mEx*(Hx(i,j)-Hx(i,j-1)) + sigma_x(i,j)*Ey(i,j);
            Ex(Nx-i+1,j) = Ex(Nx-i+1,j) + mEy*(Hy(Nx-i+1,j)-Hy(Nx-i,j)) - sigma_y(Nx-i+1,j)*Ex(Nx-i+1,j);
            Ey(Nx-i+1,j) = Ey(Nx-i+1,j) - mEx*(Hx(Nx-i+1,j)-Hx(Nx-i+1,j-1)) + sigma_x(Nx-i+1,j)*Ey(Nx-i+1,j);
        end
    end
    for i = 1:Nx
        for j = 1:npml
            Ex(i,j) = Ex(i,j) + mEy*(Hy(i,j)-Hy(i-1,j)) - sigma_y(i,j)*Ex(i,j);
            Ey(i,j) = Ey(i,j) - mEx*(Hx(i,j)-Hx(i,j-1)) + sigma_x(i,j)*Ey(i,j);
            Ex(i,Ny-j+1) = Ex(i,Ny-j+1) + mEy*(Hy(i,Ny-j+1)-Hy(i-1,Ny-j+1)) - sigma_y(i,Ny-j+1)*Ex(i,Ny-j+1);
            Ey(i,Ny-j+1) = Ey(i,Ny-j+1) - mEx*(Hx(i,Ny-j+1)-Hx(i,Ny-j)) + sigma_x(i,Ny-j+1)*Ey(i,Ny-j+1);
        end
    end
end

其中，npml表示PML层数，sigma_max和kappa_max分别是吸收率和耗散率的最大值，sigma_x、sigma_y、kappa_x和kappa_y分别是x方向和y方向的吸收率和耗散率。在时间步进中，先更新Hx和Hy，再在PML区域内更新Hx和Hy，最后更新Ex和Ey，并在PML区域内更新Ex和Ey。