贵州大学数学与统计学院，贵州 贵阳

收稿日期：2022年10月10日；录用日期：2022年11月6日；发布日期：2022年11月15日

摘要

针对一类时域Maxwell方程最优控制问题，借助Lagrange乘子法和Helmholtz分解理论，推导出由状态方程、对偶状态方程、一阶最优性条件构成的最优性系统，分别利用时域有限差分(FDTD)格式和四阶差分格式离散最优性系统，得到最优性系统的离散格式。最后，建立数值算例分别求解基于两种差分格式的Maxwell方程最优控制问题。为了避免求解大型耦合代数方程，采用迭代法进行计算，数值结果验证了理论分析结论的正确性。

关键词

Maxwell方程，最优控制，FDTD格式，四阶差分格式

Difference Schssemes for a Class of Maxwell Equation Optimal Control Problems

Lian Mei, Xianbing Luo

School of Mathematics and Statistics, Guizhou University, Guiyang Guizhou

Received: Oct. 10^th, 2022; accepted: Nov. 6^th, 2022; published: Nov. 15^th, 2022

ABSTRACT

Aiming at a class of time-domain Maxwell equation optimal control problems, an optimality system composed of equation of state, dual equation of state and first-order optimality condition is derived by using Lagransssge multiplier method and Helmholtz decomposition theory. The finite-difference time-domain (FDTD) scheme and fourth-order difference scheme are used to discretization the optimality system, respectively. The discrete scheme of the optimality system is obtained. Finally, numerical examples are established to solve the Maxwell equation optimal control problem based on the two difference schemes. In order to avoid solving large coupled algebraic equations, the iterative method is used for calculation, and the numerical results verify the correctness of the theoretical analysis.

Keywords:Maxwell’s Equation, Optimal Control, FDTD Scheme, Four-Order Difference Scheme

Copyright © 2022 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

电磁学在许多现代应用和技术中起着重要的作用，它们包括在能源科学、纳米技术、生命科学、磁约束聚变、磁悬浮技术、微波加热、传感器技术等领域的应用。近年来，在物理或者工程等实际领域中的很多问题都可以划归为Maxwell方程约束的最优控制问题。这种复杂的电磁过程的最优控制是具有挑战性的，需要仔细的数学和数值研究。对于Maxwell方程最优控制问题，通常采用有限元法(FEM)进行求解，即建立离散的目标泛函，通过有限元方法推导出状态变量和对偶状态变量的离散格式及其最优性条件进行求解。众所周知，求解Maxwell方程的方法除了FEM之外，还有矩量法(MoM)、边界元法(BEM)、差分法等。本文主要考虑Maxwell方程最优控制问题的FDTD法、四阶差分法，其中FDTD法在时间和空间上均为二阶收敛，四阶差分法在时间上二阶收敛，空间上四阶收敛。FDTD法是由K.S.Yee [1] 在1966年首次提出的一种计算电磁场的新方法，许多研究(如 [2] - [7])表明，在简单介质中，麦克斯韦方程的高阶FDTD法比Yee格式 [8] [9] 更准确和有效。特别是高阶空间差分格式降低了色散误差和相位速度各向异性误差。

随着Maxwell方程的数值计算得到了深入的研究，截止目前，Maxwell方程最优控制问题已有一些研究成果，如抛物型涡流方程 [10] [11] [12] [13] 的最优控制，这与电磁流测量中的应用问题有关(另见 [14])。对Maxwell方程的边界可控性的早期贡献可以在 [15] [16] 中找到。我们还提到了 [17] - [23] 对线性时谐波涡流方程的最优控制 [24]，而对于非线性的情况，2013年，I.Yousept等人 [25] 对拟线性椭圆型电磁场偏微分方程的最优控制展开了研究。2017年，S.Nicaise [26] 和V.Bommer [27] 等人研究了拟线性抛物型Maxwell方程的最优控制。对于文献 [27] 中研究的时域Maxwell方程最优控制问题，在求解时域Maxwell方程最优控制问题时采用FEM进行求解。根据文献检索可知，对于Maxwell方程最优控制问题，用差分法进行数值求解的研究工作很少，大多都是基于差分法进行Maxwell方程的数值求解。因此，本文讨论时域Maxwell方程最优控制问题的差分法。

本文的结构如下：第二节给出模型问题，采用先优化后离散的方案，利用偏微分方程最优控制理论给出对偶状态方程和一阶最优性条件，得到最优性系统；第三节给出最优性系统的标量形式，利用FDTD法和四阶差分法对最优性系统建立离散格式；第四节采用迭代法计算Maxwell方程最优控制问题，并给出了基于FDTD法和四阶差分法的数值实验结果，验证了两种差分格式的有效性；第五节是总结部分。

2. 模型描述

设 $\Omega \subset {ℝ}^2$ 是一个有连通Lipschitz边界 $\Gamma =\partial \Omega$ 的有界域。此外，我们考虑了一个具有连通的Lipschitz边界 ${\Gamma }_c=\partial {\Omega }_c$，其中控制域 ${\Omega }_c\subseteq \Omega$。我们的目标是在 ${\Omega }_c$ 中找到一个最优的电流密度和随时间变化的振幅，从而驱动电场和磁场达到所需的水平。

本文考虑如下Maxwell方程最优控制问题：

$\text{Minimize}\frac{1}{2}{\Vert E\left(T\right)-E_d\Vert }_{L_{ϵ}^2\left(\Omega \right)}^2+{\Vert H\left(T\right)-H_d\Vert }_{L_{\mu }^2\left(\Omega \right)}^2+\frac{{\lambda }_u}{2}{\Vert curlu\Vert }_{L^2\left({\Omega }_c\right)}^2+\frac{{\lambda }_a}{2}{\Vert a\Vert }_{L^2\left(0,T\right)}^2,$ (1)

函数 $u$ 和a满足状态方程：

$\{\begin{array}{ll}ϵE_t-curlH+\sigma E={\chi }_{{\Omega }_c}ua\hfill & \text{in}\Omega \times \left(0,T\right),\hfill \\ \mu H_t+curlE=0\hfill & \text{in}\Omega \times \left(0,T\right),\hfill \\ E\times n=0\hfill & \text{in}\Gamma \times \left(0,T\right),\hfill \\ E\left(·,0\right)=E^0\hfill & \text{in}\Omega ,\hfill \\ H\left(·,0\right)=H^0\hfill & \text{in}\Omega ,\hfill \end{array}$ (2)

对电流密度 $u=u\left(x\right)$ 的无散度约束(电荷守恒定律)以及对随时间变化的振幅 $a=a\left(t\right)$ 的约束为：

$\{\begin{array}{l}divu=0\text{in}{\Omega }_c,\\ a_{\min }\le a\left(t\right)\le a_{\max }\text{a}\text{.e}\text{.}\text{in}\left(0,T\right),\end{array}$ (3)

以上描述的是一个受Maxwell方程约束的最优控制问题，假设Maxwell方程的解满足完美导体(PEC)边界条件。在该问题中， $E$ 表示电场， $H$ 表示磁场， $ϵ$ 、 $\mu$ 、 $\sigma$ ： $\Omega \to {ℝ}^{2\times 2}$ 分别表示介电常数、磁导率、电导率， $n$ 为边界 $\Gamma$ 的外法向量， $n_c$ 为边界 ${\Gamma }_c$ 的外法向量，向量场 $E_d$ 、 $H_d$ 分别表示最终时刻 $T\in {ℝ}^{+}$ 时所期望的电场和磁场，两者均为给定函数。常数 $a_{\min }\le a_{\max }$ 表示过程中所允许振幅的最小值和最大值， ${\lambda }_u$ 、 ${\lambda }_a\in {ℝ}^{+}$ 为控制成本参数。目标泛函(1)中加权范数 ${\Vert \cdot \Vert }_{L_{ϵ}^2\left(\Omega \right)}^2={\left(ϵ\cdot ,\cdot \right)}_{L^2\left(\Omega \right)}$，${\Vert \cdot \Vert }_{L_{\mu }^2\left(\Omega \right)}^2={\left(\mu \cdot ,\cdot \right)}_{L^2\left(\Omega \right)}$。

引理1.1 [17]最优控制问题(1)存在一个最优解。

接下来引入集合：

$H\left(div;{\Omega }_c\right)=\left\{q\in L^2\left({\Omega }_c\right)|divq\in L^2\left({\Omega }_c\right)\right\},$

$H_0\left(curl;{\Omega }_c\right)=\left\{q\in H\left(curl;{\Omega }_c\right)|q\times n_c=0\text{on}{\Gamma }_c\right\},$

$H\left(div=0;{\Omega }_c\right)=\left\{q\in H\left(div;{\Omega }_c\right)|divq=0\text{in}{\Omega }_c\right\},$

${\mathcal{U}}^{feas}=H_0\left(curl;{\Omega }_c\right)\cap H\left(div=0;{\Omega }_c\right),$

${\mathcal{A}}^{feas}=\left\{a\in L^2\left(0,T\right)|a_{\min }\le a\left(t\right)\le a_{\max }\text{a}\text{.e}\text{.}\text{in}\left(0,T\right)\right\},$

利用Lagrange乘子法，以及Helmholtz分解理论，有以下结论成立：

定理1.2令 $\left(u,a\right)\in {\mathcal{U}}^{feas}\times {\mathcal{A}}^{feas}$，问题(1)~(3)存在唯一的最优解 $\left\{E,H,u,a\right\}$，当且仅当存在对偶状态变量 $\left\{K,Q\right\}$，使得 $\left\{E,H,K,Q,u,a\right\}$ 满足如下的最优性系统：

$\{\begin{array}{ll}ϵE_t-curlH+\sigma E={\chi }_{{\Omega }_c}ua\hfill & \text{in}\Omega \times \left(0,T\right),\hfill \\ \mu H_t+curlE=0\hfill & \text{in}\Omega \times \left(0,T\right),\hfill \\ E\times n=0\hfill & \text{in}\Gamma \times \left(0,T\right),\hfill \\ E\left(·,0\right)=E_0\hfill & \text{in}\Omega ,\hfill \\ H\left(·,0\right)=H_0\hfill & \text{in}\Omega ,\hfill \end{array}$ (4)

$\{\begin{array}{ll}ϵK_t-curlQ-\sigma K=0\hfill & \text{in}\Omega \times \left(0,T\right),\hfill \\ \mu Q_t+curlK=0\hfill & \text{in}\Omega \times \left(0,T\right),\hfill \\ K\times n=0\hfill & \text{in}\Gamma \times \left(0,T\right),\hfill \\ K\left(T\right)=E\left(T\right)-E_d\hfill & \text{in}\Omega ,\hfill \\ Q\left(T\right)=H\left(T\right)-H_d\hfill & \text{in}\Omega ,\hfill \end{array}$ (5)

$\{\begin{array}{ll}curlcurlu+\nabla \phi =-{\lambda }_u^{-1}{\left(K,a\right)}_{L^2(0,T)}\hfill & \text{in}{\Omega }_c,\hfill \\ divu=0\hfill & \text{in}{\Omega }_c,\hfill \\ u\times n_c=0\hfill & \text{on}{\Gamma }_c\hfill \end{array}$ (6)

$a\left(t\right)={\mathbb{P}}_{\left[a_{\min },a_{\max }\right]}\left(-\frac{1}{{\lambda }_a}{\displaystyle {\int }_{{\Omega }_c}K\left(x,y,t\right)\cdot u\left(x,y\right)\text{d}x\text{d}y}\right),\forall t\in \left[0,T\right].$ (7)

证明：定义Lagrange函数

$\begin{array}{l}L\left(E,H,u,a,K,Q\right)\\ =\frac{1}{2}{\Vert E\left(T\right)-E_d\Vert }_{L_{ϵ}^2\left(\Omega \right)}^2+{\Vert H\left(T\right)-H_d\Vert }_{L_{\mu }^2\left(\Omega \right)}^2+\frac{{\lambda }_u}{2}{\Vert curlu\Vert }_{L^2\left({\Omega }_c\right)}^2+\frac{{\lambda }_a}{2}{\Vert a\Vert }_{L^2\left(0,T\right)}^2\\ -{\displaystyle {\iint }_{\Omega \times \left(0,T\right)}\left(ϵE_t-curlH+\sigma E-{\chi }_{{\Omega }_c}ua\right)\cdot K\text{d}x\text{d}y\text{d}t}-{\displaystyle {\iint }_{\Omega \times \left(0,T\right)}\left(\mu H_t+curlE\right)\cdot Q\text{d}x\text{d}y\text{d}t}\\ =\frac{1}{2}{\Vert E\left(T\right)-E_d\Vert }_{L_{ϵ}^2\left(\Omega \right)}^2+{\Vert H\left(T\right)-H_d\Vert }_{L_{\mu }^2\left(\Omega \right)}^2+\frac{{\lambda }_u}{2}{\Vert curlu\Vert }_{L^2\left({\Omega }_c\right)}^2+\frac{{\lambda }_a}{2}{\Vert a\Vert }_{L^2\left(0,T\right)}^2\\ -{\displaystyle {\iint }_{\Omega \times \left(0,T\right)}ϵE_t\cdot K\text{d}x\text{d}y\text{d}t}-{\displaystyle {\iint }_{\Omega \times \left(0,T\right)}\sigma E\cdot K\text{d}x\text{d}y\text{d}t}-{\displaystyle {\iint }_{\Omega \times \left(0,T\right)}\left(-curlH-{\chi }_{{\Omega }_c}ua\right)\cdot K\text{d}x\text{d}y\text{d}t}\end{array}$

$\begin{array}{l}-{\displaystyle {\iint }_{\Omega \times \left(0,T\right)}curlE\cdot Q\text{d}x\text{d}y\text{d}t}-{\displaystyle {\iint }_{\Omega \times \left(0,T\right)}\mu H_t\cdot Q\text{d}x\text{d}y\text{d}t}\\ =\frac{1}{2}{\Vert E\left(T\right)-E_d\Vert }_{L_{ϵ}^2\left(\Omega \right)}^2+{\Vert H\left(T\right)-H_d\Vert }_{L_{\mu }^2\left(\Omega \right)}^2+\frac{{\lambda }_u}{2}{\Vert curlu\Vert }_{L^2\left({\Omega }_c\right)}^2+\frac{{\lambda }_a}{2}{\Vert a\Vert }_{L^2\left(0,T\right)}^2\\ -{\displaystyle {\int }_{\Omega }ϵE\left(T\right)\cdot K\left(T\right)\text{d}x\text{d}y}+{\displaystyle {\iint }_{\Omega \times \left(0,T\right)}ϵE\cdot K_t\text{d}x\text{d}y\text{d}t}-{\displaystyle {\iint }_{\Omega \times \left(0,T\right)}\sigma E\cdot K\text{d}x\text{d}y\text{d}t}\\ -{\displaystyle {\iint }_{\Omega \times \left(0,T\right)}\left(-curlH-{\chi }_{{\Omega }_c}ua\right)\cdot K\text{d}x\text{d}y\text{d}t}-{\displaystyle {\iint }_{\Omega \times \left(0,T\right)}curlQ\cdot E\text{d}x\text{d}y\text{d}t}\\ -{\displaystyle {\iint }_{\Gamma \times \left(0,T\right)}E\times n\cdot Q\text{d}x\text{d}y\text{d}t}-{\displaystyle {\iint }_{\Omega \times \left(0,T\right)}\mu H_t\cdot Q\text{d}x\text{d}y\text{d}t},\end{array}$

最后一个等式由Green公式得到，且由边界条件可知，最后一个等式的最后一项为零。

$\begin{array}{l}{D}_{E}L\left(E,H,u,a,K,Q\right)E\\ =ϵ{\displaystyle {\int }_{\Omega }\left(E\left(T\right)-E_d\right)\cdot E\left(T\right)\text{d}x\text{d}y}-ϵ{\displaystyle {\int }_{\Omega }E\left(T\right)\cdot K\left(T\right)\text{d}x\text{d}y}-{\displaystyle {\iint }_{\Omega \times \left(0,T\right)}ϵK_t\cdot E\text{d}x\text{d}y\text{d}t}\\ -{\displaystyle {\iint }_{\Omega \times \left(0,T\right)}curlQ\cdot E\text{d}x\text{d}y\text{d}t}-{\displaystyle {\iint }_{\Omega \times \left(0,T\right)}\sigma K\cdot E\text{d}x\text{d}y\text{d}t}\\ ={\displaystyle {\int }_{\Omega }\left[\left(E\left(T\right)-E_d\right)-K\left(T\right)\right]\cdot ϵE\left(T\right)\text{d}x\text{d}y}-{\displaystyle {\iint }_{\Omega \times \left(0,T\right)}\left(ϵK_t-curlQ-\sigma K\right)\cdot E\text{d}x\text{d}y\text{d}t}\\ =0,\end{array}$

则有

$K\left(T\right)=E\left(T\right)-E_d\text{in}\Omega ,$

$ϵK_t-curlQ-\sigma K=0\text{in}\Omega \times \left(0,T\right).$

类似地，

$\begin{array}{l}{D}_{H}L\left(E,H,u,a,K,Q\right)H\\ ={\displaystyle {\int }_{\Omega }\left[\left(H\left(T\right)-H_d\right)-Q\left(T\right)\right]\cdot \mu H\left(T\right)\text{d}x\text{d}y}+{\displaystyle {\iint }_{\Omega \times \left(0,T\right)}\left(\mu Q_t+curlK\right)\cdot H\text{d}x\text{d}y\text{d}t}\\ +{\displaystyle {\iint }_{\Gamma \times \left(0,T\right)}K\times n\cdot H\text{d}x\text{d}y\text{d}t}\\ =0,\end{array}$

则有

$Q\left(T\right)=H\left(T\right)-H_d\text{in}\Omega ,$

$\mu Q_t+curlK=0\text{in}\Omega \times \left(0,T\right),$

$K\times n=0\text{in}\Gamma \times \left(0,T\right).$

(6)、(7)式由Lagrange函数对控制 $u$ 和a求导，以及Helmholtz理论得到，具体证明见文献 [27]。

证毕。

3. 最优性系统的离散格式

3.1. 最优性系统

在二维直角坐标系下，设所有的物理量均与z轴坐标无关，则定理1.2中最优性系统(4)~(7)式可以表示成如下标量形式：

$\{\begin{array}{ll}ϵ\frac{\partial E_x}{\partial t}+\sigma E_x=\frac{\partial H}{\partial y}+{\chi }_{{\Omega }_c}{u}_{x}a,\hfill & \left(x,y,t\right)\in \Omega \times \left(0,T\right),\hfill \\ ϵ\frac{\partial E_y}{\partial t}+\sigma E_y=-\frac{\partial H}{\partial x}+{\chi }_{{\Omega }_c}{u}_{y}a,\hfill & \left(x,y,t\right)\in \Omega \times \left(0,T\right),\hfill \\ \mu \frac{\partial H}{\partial t}=\frac{\partial E_x}{\partial y}-\frac{\partial E_y}{\partial x},\hfill & \left(x,y,t\right)\in \Omega \times \left(0,T\right),\hfill \\ E_x\left(x,c,t\right)=E_x\left(x,d,t\right)=E_y\left(a,y,t\right)=E_y\left(b,y,t\right)=0,\hfill & \left(x,y,t\right)\in \Gamma \times \left(0,T\right),\hfill \\ E_x\left(x,y,0\right)=E_x^0\left(x,y\right),E_y\left(x,y,0\right)=E_y^0\left(x,y\right),H\left(x,y,0\right)=H^0\left(x,y\right),\hfill & \left(x,y\right)\in \Omega ,\hfill \end{array}$ (8)

$\{\begin{array}{ll}ϵ\frac{\partial K_x}{\partial t}-\sigma K_x=\frac{\partial Q}{\partial y},\hfill & \left(x,y,t\right)\in \Omega \times \left(0,T\right),\hfill \\ ϵ\frac{\partial K_y}{\partial t}-\sigma K_y=-\frac{\partial Q}{\partial x},\hfill & \left(x,y,t\right)\in \Omega \times \left(0,T\right),\hfill \\ \mu \frac{\partial Q}{\partial t}=\frac{\partial K_x}{\partial y}-\frac{\partial K_y}{\partial x},\hfill & \left(x,y,t\right)\in \Omega \times \left(0,T\right),\hfill \\ K_x\left(x,c,t\right)=K_x\left(x,d,t\right)=K_y\left(a,y,t\right)=K_y\left(b,y,t\right)=0,\hfill & \left(x,y,t\right)\in \Gamma \times \left(0,T\right),\hfill \\ K_x\left(x,y,T\right)=E_x^T\left(x,y\right)-E_d\left(x,y\right),K_y\left(x,y,T\right)=E_y^T\left(x,y\right)-E_d\left(x,y\right),\hfill & \left(x,y\right)\in \Omega ,\hfill \\ Q\left(x,y,T\right)=H^T\left(x,y\right)-H_d\left(x,y\right),\hfill & \left(x,y\right)\in \Omega ,\hfill \end{array}$ (9)

$\{\begin{array}{ll}\Delta u_x+\frac{\partial \phi }{\partial x}=-{\lambda }_u^{-1}{\left(K_x,a\right)}_{L^2\left(0,T\right)},\hfill & \left(x,y\right)\in {\Omega }_c,\hfill \\ \Delta u_y+\frac{\partial \phi }{\partial y}=-{\lambda }_u{}^{-1}{\left(K_y,a\right)}_{L^2\left(0,T\right)},\hfill & \left(x,y\right)\in {\Omega }_c,\hfill \\ \frac{\partial \phi }{\partial x}+\frac{\partial \phi }{\partial y}=0,\hfill & \left(x,y\right)\in {\Omega }_c,\hfill \\ u_x\left(x,c\right)=u_x\left(x,d\right)=u_y\left(a,y\right)=u_y\left(b,y\right)=0,\hfill & (x,y)\in {\Gamma }_c,\hfill \end{array}$ (10)

$a\left(t\right)={\mathbb{P}}_{\left[a_{\min },a_{\max }\right]}\left(-\frac{1}{{\lambda }_a}{\displaystyle {\int }_{{\Omega }_c}{K}_x\left(x,y,t\right)u_x\left(x,y\right)+K_y\left(x,y,t\right)u_y\left(x,y\right)\text{d}x\text{d}y}\right),\forall t\in \left[0,T\right],$

这里 $\left(E_x\left(x,y,t\right),E_y\left(x,y,t\right)\right)=E$ 是二维电场强度， $H\left(x,y,t\right)\equiv H_z\left(x,y,t\right)$ 表示z方向的磁场强度，电流密度 $u=\left(u_x\left(x,y\right),u_y\left(x,y\right)\right)$。 $T>0$ 为总时间，取 $\Omega =\left[a,b\right]\times \left[c,d\right]$，这里的a和b为正常数。初始条件中 $E_x^0\left(x,y\right)$，$E_y^0\left(x,y\right)$，$H^0\left(x,y\right)$ 为给定函数。

对区域 $\Omega$ 进行均匀剖分， $a=x_0<x_1<\cdots <x_{N_x}=b$，$c=y_0<y_1<\cdots <y_{N_y}=d$，对时间区域 $\left[0,T\right]$ 剖分为 $N_t$ 个网格，则离散的时间 $t_k=k\Delta t$，$\Delta t=\frac{T}{N_t}$，$k=0,1,\cdots ,N_t$。x方向网格点 $x_i=i\Delta x$，$\Delta x=\frac{b-a}{N_x}$，$i=0,1,\cdots ,N_x$，y方向网格点 $y_j=j\Delta y$，$\Delta y=\frac{d-c}{N_y}$，$j=0,1,\cdots ,N_y$。这里的Δx和Δy可取不同的值。

在接下来的两小节中讨论了最优性系统的FDTD格式和四阶差分格式，即对状态方程和对偶状态方程分别采用FDTD格式和四阶差分格式，对最优性系统中的静磁鞍点方程(10)采用MAC格式，MAC格式类似FDTD格式在空间上的离散。对于未知函数 $E_x$ 、 $K_x$ 和 $u_x$ 在x方向的中点取值， $E_y$ 、 $K_y$ 和 $u_y$ 在y方向的中点取值，H、Q和 $\phi$ 在x和y方向的中点取值， $E_x$ 、 $E_y$ 、 $K_x$ 和 $K_y$ 在时间整节点上取值，H和Q在时间半节点上取值。

3.2. 最优性系统的FDTD格式

记 $v_{i,j}^n=v\left(x_i,y_j,t_n\right)$，则

${\delta }_x{v}_{i,j}^n=\frac{v_{i+\frac{1}{2},j}^n-v_{i-\frac{1}{2},j}^n}{\Delta x},{\delta }_y{v}_{i,j}^n=\frac{v_{i,j+\frac{1}{2}}^n-v_{i,j-\frac{1}{2}}^n}{\Delta y},$

${\delta }_t{v}_{i,j}^n=\frac{v_{i,j}^{n+\frac{1}{2}}-v_{i,j}^{n-\frac{1}{2}}}{\Delta t},{\bar{\delta }}_t{v}_{i,j}^n=\frac{v_{i,j}^{n+\frac{1}{2}}+v_{i,j}^{n-\frac{1}{2}}}{\Delta t},$

因此最优性系统的离散格式为：

$\{\begin{array}{l}ϵ\frac{E_{x_{i+\frac{1}{2},j}}^{n+1}-E_{x_{i+\frac{1}{2},j}}^n}{\Delta t}+\sigma \frac{E_{x_{i+\frac{1}{2},j}}^{n+1}+E_{x_{i+\frac{1}{2},j}}^n}{2}={\delta }_y{H}_{i+\frac{1}{2},j}^{n+\frac{1}{2}}+{\chi }_{{\Omega }_c}{u}_x{}_{i+\frac{1}{2},j}{a}^n,\\ ϵ\frac{E_{y_{i,j+\frac{1}{2}}}^{n+1}-E_{y_{i,j+\frac{1}{2}}}^n}{\Delta t}+\sigma \frac{E_{y_{i,j+\frac{1}{2}}}^{n+1}+E_{y_{ij+\frac{1}{2}}}^n}{2}=-{\delta }_x{H}_{i,j+\frac{1}{2}}^{n+\frac{1}{2}}+{\chi }_{{\Omega }_c}{u}_y{}_{i,j+\frac{1}{2}}{a}^n,\\ \mu \frac{H_{i+\frac{1}{2},j+\frac{1}{2}}^{n+\frac{3}{2}}-H_{i+\frac{1}{2},j+\frac{1}{2}}^{n+\frac{1}{2}}}{\Delta t}=-{\delta }_x{E}_{y_{i+\frac{1}{2},j+\frac{1}{2}}}^{n+1}+{\delta }_y{E}_{x_{i+\frac{1}{2},j+\frac{1}{2}}}^{n+1},\\ E_{x_{i+\frac{1}{2},0}}^n=E_{x_{i+\frac{1}{2},N_y}}^n=E_{y_{0,j+\frac{1}{2}}}^n=E_{y_{N_x,j+\frac{1}{2}}}^n=0,\\ E_{x_{i+\frac{1}{2},j}}^0=E_x^0\left(x_{i+\frac{1}{2}},y_j\right),E_{y_{i,j+\frac{1}{2}}}^0=E_y^0\left(x_i,y_{j+\frac{1}{2}}\right),H_{i+\frac{1}{2},j+\frac{1}{2}}^{\frac{1}{2}}=H^0\left(x_{i+\frac{1}{2}},y_{j+\frac{1}{2}}\right),\end{array}$ (12)

$\{\begin{array}{l}ϵ\frac{K_{x_{i+\frac{1}{2},j}}^{n+1}-K_{x_{i+\frac{1}{2},j}}^n}{\Delta t}-\sigma \frac{K_{x_{i+\frac{1}{2},j}}^{n+1}+K_{x_{i+\frac{1}{2},j}}^n}{2}={\delta }_y{Q}_{i+\frac{1}{2},j}^{n+\frac{1}{2}},\\ ϵ\frac{K_{y_{i,j+\frac{1}{2}}}^{n+1}-K_{y_{i,j+\frac{1}{2}}}^n}{\Delta t}-\sigma \frac{K_{y_{i,j+\frac{1}{2}}}^{n+1}+K_{y_{ij+\frac{1}{2}}}^n}{2}=-{\delta }_x{Q}_{i,j+\frac{1}{2}}^{n+\frac{1}{2}},\\ \mu \frac{Q_{i+\frac{1}{2},j+\frac{1}{2}}^{n+\frac{3}{2}}-Q_{i+\frac{1}{2},j+\frac{1}{2}}^{n+\frac{1}{2}}}{\Delta t}=-{\delta }_x{K}_{y_{i+\frac{1}{2},j+\frac{1}{2}}}^{n+1}+{\delta }_y{K}_{x_{i+\frac{1}{2},j+\frac{1}{2}}}^{n+1},\\ K_{x_{i+\frac{1}{2},0}}^n=K_{x_{i+\frac{1}{2},N_y}}^n=K_{y_{0,j+\frac{1}{2}}}^n=K_{y_{N_x,j+\frac{1}{2}}}^n=0,\\ K_{x_{i+\frac{1}{2},j}}^T=E_x^T\left(x_{i+\frac{1}{2}},y_j\right)-E_{x_d}\left(x_{i+\frac{1}{2}},y_j\right),k_{y_{i,j+\frac{1}{2}}}^T=E_y^T\left(x_i,y_{j+\frac{1}{2}}\right)-E_{y_d}\left(x_i,y_{j+\frac{1}{2}}\right),\\ Q_{i+\frac{1}{2},j+\frac{1}{2}}^T=H^T\left(x_{i+\frac{1}{2}},y_{j+\frac{1}{2}}\right)-H_d\left(x_{i+\frac{1}{2}},y_{j+\frac{1}{2}}\right),\end{array}$ (13)

$\{\begin{array}{l}\frac{4u_{x_{i+\frac{1}{2},j}}-u_{x_{i-\frac{1}{2},j}}-u_{x_{i+\frac{3}{2},j}}-u_{x_{i+\frac{1}{2},j-1}}-u_{x_{i+\frac{1}{2},j+1}}}{\Delta x\Delta y}+\frac{{\phi }_{i+\frac{1}{2},j+\frac{1}{2}}-{\phi }_{i-\frac{1}{2},j+\frac{1}{2}}}{\Delta x}=-\frac{1}{{\lambda }_u}{\left(K_{x_{i+\frac{1}{2},j}}^n,a^n\right)}_{L^2\left(0,T\right)},\\ \frac{4u_{y_{i,j+\frac{1}{2}}}-u_{y_{i-1,j+\frac{1}{2}}}-u_{y_{i+1,j+\frac{1}{2}}}-u_{y_{i,j-\frac{1}{2}}}-u_{y_{i,j+\frac{3}{2}}}}{\Delta x\Delta y}+\frac{{\phi }_{i+\frac{1}{2},j+\frac{1}{2}}-{\phi }_{i+\frac{1}{2},j-\frac{1}{2}}}{\Delta y}=-\frac{1}{{\lambda }_u}{\left(K_{y_{i,j+\frac{1}{2}}}^n,a^n\right)}_{L^2\left(0,T\right)},\\ \frac{u_{x_{i+\frac{1}{2},j}}-u_{x_{i+\frac{3}{2},j}}}{\Delta x}+\frac{u_{y_{i,j+\frac{1}{2}}}-u_{y_{i,j-\frac{1}{2}}}}{\Delta y}=0,\\ u_{x_{i+\frac{1}{2},0}}=u_{x_{i+\frac{1}{2},N_y}}=u_{y_{0,j+\frac{1}{2}}}=u_{y_{N_x,j+\frac{1}{2}}}=0,\end{array}$ (14)

$a^n={\mathbb{P}}_{\left[a_{\min },a_{\max }\right]}\left(\frac{1}{{\lambda }_a}{\displaystyle {\int }_{{\Omega }_c}{K}_{x_{i+\frac{1}{2},j}}^n{u}_{x_{i+\frac{1}{2},j}}+K_{y_{i,j+\frac{1}{2}}}^n{u}_{y_{i,j+\frac{1}{2}}}\text{d}x\text{d}y}\right).$ (15)

3.3. 最优性系统的四阶差分格式

在Maxwell方程中，我们使用四阶差分格式来近似这些偏导数 $\frac{\partial }{\partial x}$ ：

${\delta }_x^{\left(4\right)}{v}_{i,j}=\frac{27\left(v_{i+\frac{1}{2},j}-v_{i-\frac{1}{2},j}\right)-\left(v_{i+\frac{3}{2},j}-v_{i-\frac{3}{2},j}\right)}{24}={\frac{\partial v}{\partial x}|}_{i,j}+{\rm O}\left(h_x^4\right),$

这里 $v_{i,j}$ 表示v在 $\left(x_i,y_j\right)$ 处的近似，在计算靠近边界的节点时，使用单边差分近似：

${\delta }_x^{\left(4\right)}{v}_{1,j}=\frac{-22v_{\frac{1}{2},j}+17v_{\frac{3}{2},j}+9v_{\frac{5}{2},j}-5v_{\frac{7}{2},j}+v_{\frac{9}{2},j}}{24},$

${\delta }_x^{\left(4\right)}{v}_{N_x-1,j}=\frac{-v_{N_x-\frac{9}{2},j}+5v_{N_x-\frac{7}{2},j}-9v_{N_x-\frac{5}{2},j}-17v_{N_x-\frac{3}{2},j}+22v_{N_x-\frac{1}{2},j}}{24},$

相似地可由 ${\delta }_y^{\left(4\right)}$ 近似 $\frac{\partial }{\partial y}$，这里算子可以移半个网格点，即在上述公式中，i和j可以移 $\frac{1}{2}$。时间导数项仍采用二阶差分，最后得到最优性系统四阶差分的离散格式：

$\{\begin{array}{l}ϵ\frac{E_{x_{i+\frac{1}{2},j}}^{n+1}-E_{x_{i+\frac{1}{2},j}}^n}{\Delta t}+\sigma \frac{E_{x_{i+\frac{1}{2},j}}^{n+1}+E_{x_{i+\frac{1}{2},j}}^n}{2}=\frac{1}{h_y}{\delta }_y^{\left(4\right)}{H}_{i+\frac{1}{2},j}^{n+\frac{1}{2}}+{\chi }_{{\Omega }_c}{u}_x{}_{i+\frac{1}{2},j}{a}^n,\\ ϵ\frac{E_{y_{i,j+\frac{1}{2}}}^{n+1}-E_{y_{i,j+\frac{1}{2}}}^n}{\Delta t}+\sigma \frac{E_{y_{i,j+\frac{1}{2}}}^{n+1}+E_{y_{ij+\frac{1}{2}}}^n}{2}=-\frac{1}{h_x}{\delta }_x^{\left(4\right)}{H}_{i,j+\frac{1}{2}}^{n+\frac{1}{2}}+{\chi }_{{\Omega }_c}{u}_y{}_{i,j+\frac{1}{2}}{a}^n,\\ \mu \frac{H_{i+\frac{1}{2},j+\frac{1}{2}}^{n+\frac{3}{2}}-H_{i+\frac{1}{2},j+\frac{1}{2}}^{n+\frac{1}{2}}}{\Delta t}=-\frac{1}{h_x}{\delta }_x^{\left(4\right)}{E}_{y_{i+\frac{1}{2},j+\frac{1}{2}}}^{n+1}+\frac{1}{h_y}{\delta }_y^{\left(4\right)}{E}_{x_{i+\frac{1}{2},j+\frac{1}{2}}}^{n+1},\\ E_{x_{i+\frac{1}{2},0}}^n=E_{x_{i+\frac{1}{2},N_y}}^n=E_{y_{0,j+\frac{1}{2}}}^n=E_{y_{N_x,j+\frac{1}{2}}}^n=0,\\ E_{x_{i+\frac{1}{2},j}}^0=E_x^0\left(x_{i+\frac{1}{2}},y_j\right),E_{y_{i,j+\frac{1}{2}}}^0=E_y^0\left(x_i,y_{j+\frac{1}{2}}\right),H_{i+\frac{1}{2},j+\frac{1}{2}}^{\frac{1}{2}}=H^0\left(x_{i+\frac{1}{2}},y_{j+\frac{1}{2}}\right),\end{array}$ (16)

$\{\begin{array}{l}ϵ\frac{K_{x_{i+\frac{1}{2},j}}^{n+1}-K_{x_{i+\frac{1}{2},j}}^n}{\Delta t}-\sigma \frac{K_{x_{i+\frac{1}{2},j}}^{n+1}+K_{x_{i+\frac{1}{2},j}}^n}{2}=\frac{1}{h_y}{\delta }_y^{\left(4\right)}{Q}_{i+\frac{1}{2},j}^{n+\frac{1}{2}},\\ ϵ\frac{K_{y_{i,j+\frac{1}{2}}}^{n+1}-K_{y_{i,j+\frac{1}{2}}}^n}{\Delta t}-\sigma \frac{K_{y_{i,j+\frac{1}{2}}}^{n+1}+K_{y_{ij+\frac{1}{2}}}^n}{2}=-\frac{1}{h_x}{\delta }_x^{\left(4\right)}{Q}_{i,j+\frac{1}{2}}^{n+\frac{1}{2}},\\ \mu \frac{Q_{i+\frac{1}{2},j+\frac{1}{2}}^{n+\frac{3}{2}}-Q_{i+\frac{1}{2},j+\frac{1}{2}}^{n+\frac{1}{2}}}{\Delta t}=-\frac{1}{h_x}{\delta }_x^{\left(4\right)}{K}_{y_{i+\frac{1}{2},j+\frac{1}{2}}}^{n+1}+\frac{1}{h_y}{\delta }_y^{\left(4\right)}{K}_{x_{i+\frac{1}{2},j+\frac{1}{2}}}^{n+1},\\ K_{x_{i+\frac{1}{2},0}}^n=K_{x_{i+\frac{1}{2},N_y}}^n=K_{y_{0,j+\frac{1}{2}}}^n=K_{y_{N_x,j+\frac{1}{2}}}^n=0,\\ K_{x_{i+\frac{1}{2},j}}^0=K_x^0\left(x_{i+\frac{1}{2}},y_j\right),K_{y_{i,j+\frac{1}{2}}}^0=K_y^0\left(x_i,y_{j+\frac{1}{2}}\right),Q_{i+\frac{1}{2},j+\frac{1}{2}}^{\frac{1}{2}}=Q^0\left(x_{i+\frac{1}{2}},y_{j+\frac{1}{2}}\right),\end{array}$ (17)

$\{\begin{array}{l}\frac{4u_{x_{i+\frac{1}{2},j}}-u_{x_{i-\frac{1}{2},j}}-u_{x_{i+\frac{3}{2},j}}-u_{x_{i+\frac{1}{2},j-1}}-u_{x_{i+\frac{1}{2},j+1}}}{\Delta x\Delta y}+\frac{{\phi }_{i+\frac{1}{2},j+\frac{1}{2}}-{\phi }_{i-\frac{1}{2},j+\frac{1}{2}}}{\Delta x}=-\frac{1}{{\lambda }_u}{\left(K_{x_{i+\frac{1}{2},j}}^n,a^n\right)}_{L^2\left(0,T\right)},\\ \frac{4u_{y_{i,j+\frac{1}{2}}}-u_{y_{i-1,j+\frac{1}{2}}}-u_{y_{i+1,j+\frac{1}{2}}}-u_{y_{i,j-\frac{1}{2}}}-u_{y_{i,j+\frac{3}{2}}}}{\Delta x\Delta y}+\frac{{\phi }_{i+\frac{1}{2},j+\frac{1}{2}}-{\phi }_{i+\frac{1}{2},j-\frac{1}{2}}}{\Delta y}=-\frac{1}{{\lambda }_u}{\left(K_{y_{i,j+\frac{1}{2}}}^n,a^n\right)}_{L^2\left(0,T\right)},\\ \frac{u_{x_{i+\frac{1}{2},j}}-u_{x_{i+\frac{3}{2},j}}}{\Delta x}+\frac{u_{y_{i,j+\frac{1}{2}}}-u_{y_{i,j-\frac{1}{2}}}}{\Delta y}=0,\\ u_{x_{i+\frac{1}{2},0}}=u_{x_{i+\frac{1}{2},N_y}}=u_{y_{0,j+\frac{1}{2}}}=u_{y_{N_x,j+\frac{1}{2}}}=0,\end{array}$ (18)