二维偏微分形式麦克斯韦方程组紧差分格式数值解与解析解误差并显示收敛精度阶数matlab

二维偏微分形式的麦克斯韦方程组可以写成以下形式：

$$ \begin{aligned} \frac{\partial E_x}{\partial t} &= \frac{1}{\epsilon_0} \left(\frac{\partial H_z}{\partial y} - \frac{\partial H_y}{\partial z}\right) \\ \frac{\partial E_y}{\partial t} &= \frac{1}{\epsilon_0} \left(\frac{\partial H_x}{\partial z} - \frac{\partial H_z}{\partial x}\right) \\ \frac{\partial E_z}{\partial t} &= \frac{1}{\epsilon_0} \left(\frac{\partial H_y}{\partial x} - \frac{\partial H_x}{\partial y}\right) \\ \frac{\partial H_x}{\partial t} &= \frac{1}{\mu_0} \left(\frac{\partial E_z}{\partial y} - \frac{\partial E_y}{\partial z}\right) \\ \frac{\partial H_y}{\partial t} &= \frac{1}{\mu_0} \left(\frac{\partial E_x}{\partial z} - \frac{\partial E_z}{\partial x}\right) \\ \frac{\partial H_z}{\partial t} &= \frac{1}{\mu_0} \left(\frac{\partial E_y}{\partial x} - \frac{\partial E_x}{\partial y}\right) \end{aligned} $$

其中，$\epsilon_0$ 和 $\mu_0$ 分别为真空介电常数和真空磁导率。为了数值求解，我们采用了紧差分格式。具体来说，我们采用了二阶中心差分格式，即：

$$ \frac{\partial u}{\partial x} \approx \frac{u_{i+1,j} - u_{i-1,j}}{2\Delta x}, \quad \frac{\partial u}{\partial y} \approx \frac{u_{i,j+1} - u_{i,j-1}}{2\Delta y} $$

将其代入上面的麦克斯韦方程组中，得到：

$$ \begin{aligned} \frac{E_{x,i,j}^{n+1} - E_{x,i,j}^{n}}{\Delta t} &= \frac{1}{\epsilon_0} \left(\frac{H_{z,i,j+1}^{n+1} - H_{z,i,j-1}^{n+1} - H_{z,i,j+1}^{n} + H_{z,i,j-1}^{n}}{2\Delta y} - \frac{H_{y,i,j+1}^{n+1} - H_{y,i,j-1}^{n+1} - H_{y,i,j+1}^{n} + H_{y,i,j-1}^{n}}{2\Delta z}\right) \\ \frac{E_{y,i,j}^{n+1} - E_{y,i,j}^{n}}{\Delta t} &= \frac{1}{\epsilon_0} \left(\frac{H_{x,i,j+1}^{n+1} - H_{x,i,j-1}^{n+1} - H_{x,i,j+1}^{n} + H_{x,i,j-1}^{n}}{2\Delta z} - \frac{H_{z,i+1,j}^{n+1} - H_{z,i-1,j}^{n+1} - H_{z,i+1,j}^{n} + H_{z,i-1,j}^{n}}{2\Delta x}\right) \\ \frac{E_{z,i,j}^{n+1} - E_{z,i,j}^{n}}{\Delta t} &= \frac{1}{\epsilon_0} \left(\frac{H_{y,i+1,j}^{n+1} - H_{y,i-1,j}^{n+1} - H_{y,i+1,j}^{n} + H_{y,i-1,j}^{n}}{2\Delta x} - \frac{H_{x,i+1,j}^{n+1} - H_{x,i-1,j}^{n+1} - H_{x,i+1,j}^{n} + H_{x,i-1,j}^{n}}{2\Delta y}\right) \\ \frac{H_{x,i,j}^{n+1} - H_{x,i,j}^{n}}{\Delta t} &= \frac{1}{\mu_0} \left(\frac{E_{z,i,j+1}^{n+1} - E_{z,i,j-1}^{n+1} - E_{z,i,j+1}^{n} + E_{z,i,j-1}^{n}}{2\Delta y} - \frac{E_{y,i+1,j}^{n+1} - E_{y,i-1,j}^{n+1} - E_{y,i+1,j}^{n} + E_{y,i-1,j}^{n}}{2\Delta z}\right) \\ \frac{H_{y,i,j}^{n+1} - H_{y,i,j}^{n}}{\Delta t} &= \frac{1}{\mu_0} \left(\frac{E_{x,i+1,j}^{n+1} - E_{x,i-1,j}^{n+1} - E_{x,i+1,j}^{n} + E_{x,i-1,j}^{n}}{2\Delta z} - \frac{E_{z,i,j+1}^{n+1} - E_{z,i,j-1}^{n+1} - E_{z,i,j+1}^{n} + E_{z,i,j-1}^{n}}{2\Delta x}\right) \\ \frac{H_{z,i,j}^{n+1} - H_{z,i,j}^{n}}{\Delta t} &= \frac{1}{\mu_0} \left(\frac{E_{y,i+1,j}^{n+1} - E_{y,i-1,j}^{n+1} - E_{y,i+1,j}^{n} + E_{y,i-1,j}^{n}}{2\Delta x} - \frac{E_{x,i,j+1}^{n+1} - E_{x,i,j-1}^{n+1} - E_{x,i,j+1}^{n} + E_{x,i,j-1}^{n}}{2\Delta y}\right) \end{aligned} $$

其中，$E_{x,i,j}^n$ 表示在 $x=i\Delta x$、$y=j\Delta y$ 和 $t=n\Delta t$ 时刻的电场 $E_x$ 值，$H_{z,i,j}^{n+1}$ 表示在 $x=i\Delta x$、$y=j\Delta y$ 和 $t=(n+1)\Delta t$ 时刻的磁场 $H_z$ 值，等等。

我们可以使用 MATLAB 实现该差分格式，具体代码如下：

% 定义计算区域
Lx = 1; % 区域长度（x方向）
Ly = 1; % 区域长度（y方向）
Nx = 50; % 离散点数（x方向）
Ny = 50; % 离散点数（y方向）
dx = Lx / Nx; % 离散步长（x方向）
dy = Ly / Ny; % 离散步长（y方向）

% 定义常数
c = 3e8; % 光速
mu0 = pi * 4e-7; % 真空磁导率
eps0 = 1 / (mu0 * c^2); % 真空介电常数

% 初始化电场和磁场
Ex = zeros(Nx+1, Ny+1); % 电场（x分量）
Ey = zeros(Nx+1, Ny+1); % 电场（y分量）
Ez = zeros(Nx+1, Ny+1); % 电场（z分量）
Hx = zeros(Nx+1, Ny+1); % 磁场（x分量）
Hy = zeros(Nx+1, Ny+1); % 磁场（y分量）
Hz = zeros(Nx+1, Ny+1); % 磁场（z分量）

% 定义模拟参数
dt = 0.99 / (c * sqrt(1/dx^2 + 1/dy^2)); % 时间步长
T = 500 * dt; % 总模拟时间
Nt = round(T / dt); % 时间步数

% 定义边界条件
Ex(1,:,:) = 0;
Ex(end,:,:) = 0;
Ex(:,1,:) = 0;
Ex(:,end,:) = 0;
Ey(1,:,:) = 0;
Ey(end,:,:) = 0;
Ey(:,1,:) = 0;
Ey(:,end,:) = 0;
Ez(1,:,:) = 0;
Ez(end,:,:) = 0;
Ez(:,1,:) = 0;
Ez(:,end,:) = 0;
Hx(1,:,:) = 0;
Hx(end,:,:) = 0;
Hx(:,1,:) = 0;
Hx(:,end,:) = 0;
Hy(1,:,:) = 0;
Hy(end,:,:) = 0;
Hy(:,1,:) = 0;
Hy(:,end,:) = 0;
Hz(1,:,:) = 0;
Hz(end,:,:) = 0;
Hz(:,1,:) = 0;
Hz(:,end,:) = 0;

% 循环计算电磁场
for n = 1:Nt
    % 计算电场
    Ex(2:end-1, 2:end-1) = Ex(2:end-1, 2:end-1) + dt/eps0 * (Hz(2:end-1, 3:end) - Hz(2:end-1, 1:end-2))/(2*dy) - dt/eps0 * (Hy(3:end, 2:end-1) - Hy(1:end-2, 2:end-1))/(2*dx);
    Ey(2:end-1, 2:end-1) = Ey(2:end-1, 2:end-1) + dt/eps0 * (Hx(3:end, 2:end-1) - Hx(1:end-2, 2:end-1))/(2*dy) - dt/eps0 * (Hz(2:end-1, 3:end) - Hz(2:end-1, 1:end-2))/(2*dx);
    Ez(2:end-1, 2:end-1) = Ez(2:end-1, 2:end-1) + dt/eps0 * (Hy(2:end-1, 3:end) - Hy(2:end-1, 1:end-2))/(2*dx) - dt/eps0 * (Hx(3:end, 2:end-1) - Hx(1:end-2, 2:end-1))/(2*dy);
    
    % 计算磁场
    Hx(2:end-1, 2:end-1) = Hx(2:end-1, 2:end-1) - dt/mu0 * (Ez(2:end-1, 3:end) - Ez(2:end-1, 1:end-2))/(2*dy) + dt/mu0 * (Ey(3:end, 2:end-1) - Ey(1:end-2, 2:end-1))/(2*dx);
    Hy(2:end-1, 2:end-1) = Hy(2:end-1, 2:end-1) - dt/mu0 * (Ex(3:end, 2:end-1) - Ex(1:end-2, 2:end-1))/(2*dy) + dt/mu0 * (Ez(2:end-1, 3:end) - Ez(2:end-1, 1:end-2))/(2*dx);
    Hz(2:end-1, 2:end-1) = Hz(2:end-1, 2:end-1) - dt/mu0 * (Ey(2:end-1, 3:end) - Ey(2:end-1, 1:end-2))/(2*dx) + dt/mu0 * (Ex(3:end, 2:end-1) - Ex(1:end-2, 2:end-1))/(2*dy);
    
    % 更新边界条件
    Ex(1,:,:) = 0;
    Ex(end,:,:) = 0;
    Ex(:,1,:) = 0;
    Ex(:,end,:) = 0;
    Ey(1,:,:) = 0;
    Ey(end,:,:) = 0;
    Ey(:,1,:) = 0;
    Ey(:,end,:) = 0;
    Ez(1,:,:) = 0;
    Ez(end,:,:) = 0;
    Ez(:,1,:) = 0;
    Ez(:,end,:) = 0;
    Hx(1,:,:) = 0;
    Hx(end,:,:) = 0;
    Hx(:,1,:) = 0;
    Hx(:,end,:) = 0;
    Hy(1,:,:) = 0;
    Hy(end,:,:) = 0;
    Hy(:,1,:) = 0;
    Hy(:,end,:) = 0;
    Hz(1,:,:) = 0;
    Hz(end,:,:) = 0;
    Hz(:,1,:) = 0;
    Hz(:,end,:) = 0;
    
    % 输出进度
    if mod(n, 10) == 0
        fprintf('Progress: %d / %d\n', n, Nt);
    end
end

% 计算解析解
x = linspace(0, Lx, Nx+1);
y = linspace(0, Ly, Ny+1);
[X, Y] = meshgrid(x, y);
T = linspace(0, T, Nt+1);
[E_analytic, H_analytic] = maxwell_analytic_2d(X, Y, T, eps0, mu0);

% 计算误差
E_error = Ex - E_analytic(:, :, end);
H_error = Hx - H_analytic(:, :, end);

% 计算收敛精度阶数
E_error_norm = zeros(1, Nt+1);
for n = 1:Nt+1
    E_error_norm(n) = norm(E_error(:, :, n), 'fro');
end
E_order = polyfit(log(dt./(1:Nt+1)), log(E_error_norm), 1);
E_order = E_order(1);

H_error_norm = zeros(1, Nt+1);
for n = 1:Nt+1
    H_error_norm(n) = norm(H_error(:, :, n), 'fro');
end
H_order = polyfit(log(dt./(1:Nt+1)), log(H_error_norm), 1);
H_order = H_order(1);

% 绘制图形
figure;
subplot(2, 2, 1);
pcolor(X, Y, Ex(:, :, end)); shading flat;
title('Electric field (x component)');
xlabel('x');
ylabel('y');
colorbar;
subplot(2, 2, 2);
pcolor(X, Y, Ey(:, :, end)); shading flat;
title('Electric field (y component)');
xlabel('x');
ylabel('y');
colorbar;
subplot(2, 2, 3);
pcolor(X, Y, Hx(:, :, end)); shading flat;
title('Magnetic field (x component)');
xlabel('x');
ylabel('y');
colorbar;
subplot(2, 2, 4);
pcolor(X, Y, Hy(:, :, end)); shading flat;
title('Magnetic field (y component)');
xlabel('x');
ylabel('y');
colorbar;

figure;
subplot(1, 2, 1);
loglog(dt./(1:Nt+1), E_error_norm, '-o');
title(sprintf('Electric field error (order = %.2f)', E_order));
xlabel('\Delta t');
ylabel('Error');
subplot(1, 2, 2);
loglog(dt./(1:Nt+1), H_error_norm, '-o');
title(sprintf('Magnetic field error (order = %.2f)', H_order));
xlabel('\Delta t');
ylabel('Error');

其中，`maxwell_analytic_2d` 是一个 MATLAB 函数，用于计算麦克斯韦方程组的解析解。具体代码如下：

function [E, H] = maxwell_analytic_2d(X, Y, T, eps0, mu0)
% 计算麦克斯韦方程组的解析解

% 计算模拟参数
c = 1 / sqrt(eps0 * mu0); % 光速
omega = 2 * pi * c / 1e-6; % 角频率

% 计算电磁场
E = zeros(size(X, 1), size(X, 2), length(T));
H = zeros(size(X, 1), size(X, 2), length(T));
for n = 1:length(T)
    t = T(n);
    E(:, :, n) = sin(omega * t) * sin(pi * X) .* sin