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

二维麦克斯韦方程组紧差分格式的数值解误差可以通过计算数值解与精确解之间的差异得到。但是，由于二维麦克斯韦方程组的精确解很难求解，因此常常采用比较精确的数值解作为参考，计算数值解的误差。一般来说，采用 $L^2$ 范数或 $H^1$ 范数计算误差。

假设采用 $L^2$ 范数计算误差，即计算数值解 $u_h$ 与精确解 $u$ 之间的 $L^2$ 范数的差异：

$$\epsilon_{L^2} = \sqrt{\frac{\sum_{i,j}(u_{h,ij}-u_{ij})^2}{\sum_{i,j}u_{ij}^2}}$$

其中，$u_{h,ij}$ 表示数值解在网格点 $(i,j)$ 处的值，$u_{ij}$ 表示精确解在网格点 $(i,j)$ 处的值。

可以通过调整网格大小 $h$，计算不同精度的数值解，并计算相应的误差，从而确定数值解的收敛精度阶数。通常情况下，误差与网格大小 $h$ 的关系可以近似表示为：

$$\epsilon_{L^2} \approx Ch^p$$

其中，$C$ 是一个常数，$p$ 表示收敛精度阶数。

下面是一个在 MATLAB 中计算二维麦克斯韦方程组紧差分格式数值解误差并显示收敛精度阶数的示例代码：

%% 二维麦克斯韦方程组紧差分格式数值解误差及收敛精度阶数
clc; clear; close all;

%% 定义常数
epsilon0 = 8.854e-12; % 真空介电常数
mu0 = pi * 4e-7; % 真空磁导率
c = 1 / sqrt(epsilon0 * mu0); % 光速

%% 定义计算区域
Lx = 1; Ly = 1; % 区域大小
Nx = 50; Ny = 50; % 网格数
hx = Lx / (Nx - 1); hy = Ly / (Ny - 1); % 网格步长
x = linspace(0, Lx, Nx); y = linspace(0, Ly, Ny); % 网格点坐标
[X, Y] = meshgrid(x, y);

%% 初始化场变量
Ex = zeros(Nx, Ny); Ey = zeros(Nx, Ny); % 电场
Hx = zeros(Nx, Ny); Hy = zeros(Nx, Ny); % 磁场

%% 计算时间
T = 0.5; % 计算时间
dt = 0.9 * min(hx, hy) / c; % 时间步长
Nt = round(T / dt); % 时间步数
t = linspace(0, T, Nt); % 时间点

%% 计算初始条件
Ex(:, 1) = sin(pi * X / Lx); % x 方向电场
Hy(1, :) = sin(pi * Y / Ly); % y 方向磁场

%% 计算数值解
for n = 1 : Nt
    % 更新电场
    Ex(:, 2 : Ny - 1) = Ex(:, 2 : Ny - 1) + dt / epsilon0 / hx * (Hy(:, 2 : Ny - 1) - Hy(:, 1 : Ny - 2));
    Ey(2 : Nx - 1, :) = Ey(2 : Nx - 1, :) - dt / epsilon0 / hy * (Hx(2 : Nx - 1, :) - Hx(1 : Nx - 2, :));
    % 更新磁场
    Hx(2 : Nx - 1, :) = Hx(2 : Nx - 1, :) - dt / mu0 / hy * (Ey(2 : Nx - 1, :) - Ey(1 : Nx - 2, :));
    Hy(:, 2 : Ny - 1) = Hy(:, 2 : Ny - 1) + dt / mu0 / hx * (Ex(:, 2 : Ny - 1) - Ex(:, 1 : Ny - 2));
end

%% 计算误差及收敛精度阶数
u = sin(pi * X / Lx) .* sin(pi * Y / Ly); % 精确解
uh = Ex; % 数值解
epsilon = sqrt(sum(sum((uh - u).^2)) / sum(sum(u.^2))); % 误差
fprintf('误差 = %e\n', epsilon);

h_list = [hx, hx / 2, hx / 4, hx / 8]; % 不同网格步长
epsilon_list = zeros(size(h_list)); % 不同网格步长下的误差
for i = 1 : length(h_list)
    hx = h_list(i); hy = h_list(i); % 网格步长
    Nx = round(Lx / hx) + 1; Ny = round(Ly / hy) + 1; % 网格数
    x = linspace(0, Lx, Nx); y = linspace(0, Ly, Ny); % 网格点坐标
    [X, Y] = meshgrid(x, y);
    Ex = zeros(Nx, Ny); Ey = zeros(Nx, Ny); % 电场
    Hx = zeros(Nx, Ny); Hy = zeros(Nx, Ny); % 磁场
    Ex(:, 1) = sin(pi * X / Lx); % x 方向电场
    Hy(1, :) = sin(pi * Y / Ly); % y 方向磁场
    for n = 1 : Nt
        % 更新电场
        Ex(:, 2 : Ny - 1) = Ex(:, 2 : Ny - 1) + dt / epsilon0 / hx * (Hy(:, 2 : Ny - 1) - Hy(:, 1 : Ny - 2));
        Ey(2 : Nx - 1, :) = Ey(2 : Nx - 1, :) - dt / epsilon0 / hy * (Hx(2 : Nx - 1, :) - Hx(1 : Nx - 2, :));
        % 更新磁场
        Hx(2 : Nx - 1, :) = Hx(2 : Nx - 1, :) - dt / mu0 / hy * (Ey(2 : Nx - 1, :) - Ey(1 : Nx - 2, :));
        Hy(:, 2 : Ny - 1) = Hy(:, 2 : Ny - 1) + dt / mu0 / hx * (Ex(:, 2 : Ny - 1) - Ex(:, 1 : Ny - 2));
    end
    uh = Ex; % 数值解
    u = sin(pi * X / Lx) .* sin(pi * Y / Ly); % 精确解
    epsilon_list(i) = sqrt(sum(sum((uh - u).^2)) / sum(sum(u.^2))); % 误差
end
p_list = log(epsilon_list(1 : end - 1) ./ epsilon_list(2 : end)) ./ log(h_list(1 : end - 1) ./ h_list(2 : end)); % 收敛精度阶数
fprintf('收敛精度阶数 = %f\n', p_list(end));
figure;
loglog(h_list, epsilon_list, '-o');
xlabel('网格步长'); ylabel('误差');
title('紧差分格式收敛精度阶数');

在上述代码中，首先定义了常数、计算区域、时间和初始条件等参数。然后，通过循环更新电场和磁场，得到数值解。接着，计算数值解与精确解之间的误差，并通过调整网格步长计算不同精度的数值解，从而确定收敛精度阶数。最后，通过绘制误差与网格步长的对数-log图，展示收敛精度阶数。