二维偏微分形式麦克斯韦方程组紧差分格式matlab

二维麦克斯韦方程组的紧差分格式如下：

$$\frac{E_x^{n+1/2,m}-E_x^{n-1/2,m}}{\Delta t}-\frac{1}{\epsilon}\left(\frac{H_z^{n,m+1/2}-H_z^{n,m-1/2}}{\Delta y}\right)=0$$

$$\frac{E_y^{n+1/2,m}-E_y^{n-1/2,m}}{\Delta t}-\frac{1}{\epsilon}\left(\frac{H_z^{n+1/2,m}-H_z^{n-1/2,m}}{\Delta x}\right)=0$$

$$\frac{H_z^{n,m+1/2}-H_z^{n,m-1/2}}{\Delta t}-\frac{1}{\mu}\left(\frac{E_y^{n+1/2,m}-E_y^{n-1/2,m}}{\Delta x}-\frac{E_x^{n,m+1/2}-E_x^{n,m-1/2}}{\Delta y}\right)=0$$

其中，$E_x^{n+1/2,m}$ 表示 $E_x$ 在时间 $t^{n+1/2}$ 和位置 $(m\Delta x,(n+1/2)\Delta t)$ 处的值，$H_z^{n,m+1/2}$ 表示 $H_z$ 在时间 $t^{n}$ 和位置 $((m+1/2)\Delta x,n\Delta t)$ 处的值，$\Delta x$ 和 $\Delta y$ 分别为空间步长，$\Delta t$ 为时间步长，$\epsilon$ 和 $\mu$ 分别为介电常数和磁导率。

将上述差分格式转化为矩阵形式，得到：

$$\begin{pmatrix}
1 & 0 & 0 & -\frac{\Delta t}{\epsilon\Delta y} & 0 & 0 \\
0 & 1 & 0 & 0 & -\frac{\Delta t}{\epsilon\Delta x} & 0 \\
0 & 0 & 1 & \frac{\Delta t}{\mu\Delta x} & -\frac{\Delta t}{\mu\Delta y} & 0 \\
-\frac{\Delta t}{\epsilon\Delta y} & 0 & 0 & 1 & 0 & 0 \\
0 & -\frac{\Delta t}{\epsilon\Delta x} & 0 & 0 & 1 & 0 \\
0 & 0 & 0 & \frac{\Delta t}{\mu\Delta y} & -\frac{\Delta t}{\mu\Delta x} & 1
\end{pmatrix}\begin{pmatrix}
E_x^{n+1/2,m}\\
E_y^{n+1/2,m}\\
H_z^{n,m+1/2}\\
E_x^{n-1/2,m}\\
E_y^{n-1/2,m}\\
H_z^{n,m-1/2}
\end{pmatrix}=\begin{pmatrix}
0\\
0\\
0\\
E_x^{n+1/2,m}+E_x^{n-1/2,m}\\
E_y^{n+1/2,m}+E_y^{n-1/2,m}\\
H_z^{n,m+1/2}+H_z^{n,m-1/2}
\end{pmatrix}$$

其中，未知量为 $E_x^{n+1/2,m}$、$E_y^{n+1/2,m}$ 和 $H_z^{n,m+1/2}$，等号右侧的向量中包含已知量。

下面是使用 MATLAB 实现的代码：

% 输入参数
Nx = 100; % x方向空间步数
Ny = 100; % y方向空间步数
Nt = 500; % 时间步数
dt = 0.1; % 时间步长
dx = 0.1; % x方向空间步长
dy = 0.1; % y方向空间步长
mu = 1; % 磁导率
epsilon = 1; % 介电常数

% 初始化电场和磁场
Ex = zeros(Nx+1,Ny+1);
Ey = zeros(Nx+1,Ny+1);
Hz = zeros(Nx+1,Ny+1);

% 构造系数矩阵
A = [1, 0, 0, -dt/(epsilon*dy), 0, 0; 
     0, 1, 0, 0, -dt/(epsilon*dx), 0; 
     0, 0, 1, dt/(mu*dx), -dt/(mu*dy), 0; 
     -dt/(epsilon*dy), 0, 0, 1, 0, 0; 
     0, -dt/(epsilon*dx), 0, 0, 1, 0; 
     0, 0, 0, dt/(mu*dy), -dt/(mu*dx), 1];
 
% 构造右侧向量
b = zeros(6*(Nx+1)*(Ny+1),1);
for n = 2:Nt+1
    for m = 2:Nx
        for l = 2:Ny
            % 更新电场和磁场
            Ex(m,l) = Ex(m,l) + dt/epsilon*(Hz(m,l)-Hz(m,l-1))/dy;
            Ey(m,l) = Ey(m,l) - dt/epsilon*(Hz(m,l)-Hz(m-1,l))/dx;
            Hz(m,l) = Hz(m,l) + dt/mu*(Ey(m+1,l)-Ey(m,l)+Ex(m,l)-Ex(m,l+1))/dx;
            
            % 构造右侧向量
            b((m-2)*(Ny+1)*6+(l-2)*6+1) = Ex(m,l);
            b((m-2)*(Ny+1)*6+(l-2)*6+2) = Ey(m,l);
            b((m-2)*(Ny+1)*6+(l-2)*6+3) = Hz(m,l);
            b((m-2)*(Ny+1)*6+(l-2)*6+4) = Ex(m,l+1)+Ex(m,l-1);
            b((m-2)*(Ny+1)*6+(l-2)*6+5) = Ey(m+1,l)+Ey(m-1,l);
            b((m-2)*(Ny+1)*6+(l-2)*6+6) = Hz(m,l+1)+Hz(m,l-1);
        end
    end
    
    % 求解线性方程组
    x = A\b;
    
    % 更新电场和磁场
    for m = 2:Nx
        for l = 2:Ny
            Ex(m,l) = x((m-2)*(Ny+1)*6+(l-2)*6+1);
            Ey(m,l) = x((m-2)*(Ny+1)*6+(l-2)*6+2);
            Hz(m,l) = x((m-2)*(Ny+1)*6+(l-2)*6+3);
        end
    end
    
    % 绘制电场和磁场的图像
    surf(Ex')
    title(['Ex at t = ',num2str((n-1)*dt)])
    xlabel('x')
    ylabel('y')
    zlabel('Ex')
    pause(0.01)
end