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

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

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

其中，$E_x$、$E_y$、$H_x$、$H_y$、$E_z$、$H_z$ 分别表示电场和磁场在 $x$ 和 $y$ 方向的分量以及在 $z$ 方向的分量。

为了求解这个方程组的数值解，我们可以采用紧差分格式。具体来说，我们可以将时间轴上的 $t$ 分成若干个小时间步长 $\Delta t$，对于每个时间步长，我们采用中心差分格式来计算 $E_x$、$E_y$、$H_x$、$H_y$、$E_z$、$H_z$ 在该时间步长的值。具体来说，我们可以采用以下的差分格式：

$$
\begin{aligned}
\frac{E_x^{n+1}_{i,j}-E_x^{n}_{i,j}}{\Delta t} &= \frac{1}{\epsilon}\left(\frac{H_z^{n+1/2}_{i,j+1/2}-H_z^{n+1/2}_{i,j-1/2}}{\Delta y}\right) \\
\frac{E_y^{n+1}_{i,j}-E_y^{n}_{i,j}}{\Delta t} &= -\frac{1}{\epsilon}\left(\frac{H_z^{n+1/2}_{i+1/2,j}-H_z^{n+1/2}_{i-1/2,j}}{\Delta x}\right) \\
\frac{H_x^{n+1/2}_{i+1/2,j}-H_x^{n+1/2}_{i-1/2,j}}{\Delta t} &= -\frac{1}{\mu}\left(\frac{E_z^{n+1}_{i,j+1}-E_z^{n+1}_{i,j-1}}{\Delta y}\right) \\
\frac{H_y^{n+1/2}_{i,j+1/2}-H_y^{n+1/2}_{i,j-1/2}}{\Delta t} &= \frac{1}{\mu}\left(\frac{E_z^{n+1}_{i+1,j}-E_z^{n+1}_{i-1,j}}{\Delta x}\right) \\
\frac{E_z^{n+1}_{i,j}-E_z^{n}_{i,j}}{\Delta t} &= -\frac{1}{\epsilon}\left(\frac{H_y^{n+1/2}_{i+1/2,j}-H_y^{n+1/2}_{i-1/2,j}}{\Delta x}-\frac{H_x^{n+1/2}_{i,j+1/2}-H_x^{n+1/2}_{i,j-1/2}}{\Delta y}\right) \\
\frac{H_z^{n+1/2}_{i,j+1/2}-H_z^{n+1/2}_{i,j-1/2}}{\Delta t} &= \frac{1}{\mu}\left(\frac{E_y^{n+1}_{i+1,j}-E_y^{n+1}_{i-1,j}}{\Delta x}-\frac{E_x^{n+1}_{i,j+1}-E_x^{n+1}_{i,j-1}}{\Delta y}\right) \\
\end{aligned}
$$

其中，$n$ 表示时间步长的编号，$i$ 和 $j$ 分别表示 $x$ 和 $y$ 方向上的网格节点编号，$\Delta x$ 和 $\Delta y$ 分别表示 $x$ 和 $y$ 方向上的网格间距。

为了验证我们的数值解的正确性，我们可以将数值解与解析解进行比较，计算误差和显示收敛精度阶数。具体来说，我们可以采用以下的方式计算误差：

$$
\mathrm{error} = \frac{\|E_z^\mathrm{numerical}-E_z^\mathrm{analytical}\|}{\|E_z^\mathrm{analytical}\|}
$$

其中，$E_z^\mathrm{numerical}$ 表示数值解，$E_z^\mathrm{analytical}$ 表示解析解。

为了计算显示收敛精度阶数，我们可以采用如下方式：

$$
\mathrm{order} = \frac{\log_2(\mathrm{error}_1/\mathrm{error}_2)}{\log_2(\mathrm{h}_1/\mathrm{h}_2)}
$$

其中，$\mathrm{error}_1$ 和 $\mathrm{error}_2$ 分别表示两个不同网格尺寸下的误差，$\mathrm{h}_1$ 和 $\mathrm{h}_2$ 分别表示两个不同网格尺寸。

在 MATLAB 中，我们可以使用以下代码来实现这个数值算法和误差计算：