Maxwell方程是电磁场理论中非常重要的偏微分方程组, 揭示了电与磁的相互转化、相互依赖的关系, 已应用于现代科学和工程技术的很多领域^[1⇓-3]. 有限差分时域(FDTD)方法(也称为Yee方法)^[4]是电磁场计算中非常有效的数值方法, 被广泛应用和推广.

非线性Maxwell方程在科学工程相关领域普遍存在. 其中, 具有依赖于场源的非线性电导率的Maxwell方程应用在多个物理模型中, 例如第二型超导体的本构定律、NbSe3的电荷密度波态等^[5-6]. 关于非线性Maxwell方程弱解的存在唯一性以及当$t\rightarrow\infty$时解的渐近性, 很多学者已经做了重要的研究工作^[7⇓⇓-10]. 非线性Maxwell方程的数值求解也得到了很多学者的关注. Durand等^[11]针对一般情况的非线性Maxwell方程, 证明了弱解的存在性,$\!\!$并提出了时间上采用向后Euler离散方法空间上采用混合有限元法的数值方法. $\!\!$Yao等^[12]讨论了非线性Maxwell方程的数值解法, 时间上采用三阶线性向后差分方法(BDF), 空间上采用Nédélec有限元方法, 得到了$L^2$ -范数意义下的最优误差估计. Bokil等^[13]发展了一种新的空间高阶的能量稳定的FDTD方法求解一维非线性Maxwell方程. $\!\!$Huang等^[14]分析一种有限元方法来求解非线性Kerr介质中的Maxwell 方程, 采用了非迭代的方法求解此非线性问题, 并得到空间上的最优收敛阶, 时间上的二阶收敛性.

非一致介质Maxwell方程中的介电常数和磁导率在跨过不同的介质时会产生跳跃, 从而使解产生间断. 电磁场间断解问题的数值计算, 有一定的难度和挑战性, 促使很多学者从事这方面的研究并提出了一些有效的数值方法, 如高阶嵌入FDTD和导数匹配法、匹配和沉浸交界面(IIM)法、间断Galerkin(DG)方法等^{[15⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓-26]}. Zhao等^[18]讨论了非一致介质中具有间断解的Maxwell方程的数值解法, 提出了高阶的隐式分层的导数匹配方法(HDM-FDTD). 此方法在介质交界面处增加虚拟点和附加的连接条件, 利用高阶跳跃条件来得到高精度, 如采用多于$(2m-1)$阶的物理跳跃条件和标准的中心差分逼近方法处理$2m$个跳跃条件. 这样使得数值格式变得稍微有些复杂.

谱方法数值求解非线性系统的理论不断发展完善起来^{[27⇓⇓⇓-31]}. Ma^[32-33]讨论了非线性守恒方程的Chebyshev-Legendre耦合谱方法, 提出谱粘性(CLSV)和超谱粘性的方法. Ma和Sun^[34]提出了一种广义的Legendre-Petrov-Galerkin谱方法求解KdV方程, 对非线性项采用了Chebyshev谱配置方法, 并证明了半离散和全离散格式$L^2$ -范数下的最优误差估计. 多区域谱方法是求解非一致介质问题的有效方法之一. 这种方法把整个区域分成多个子区域, 采用分片多项式逼近每个子区域上的解. Ma等^[35]提出了多区域的Legendre-Galerkin Chebyshev配置方法(MLGCC)求解非一致介质中具有间断解的一维Maxwell方程, 保持了谱精度.

基于上述研究工作, 本文考虑了非线性非一致介质二维Maxwell方程的数值求解方法, 提出了leap-frog Crank-Nicolson(LFCN)多区域Legendre-tau配置(MLTC)谱方法. 本数值方法空间上达到谱精度, 时间上是二阶精度. 此数值方法采用tau方法思想, 针对电场和磁场的不同特性, 定义不同次多项式空间来逼近电磁场. 利用显隐数值格式, 对线性和非线性项采用不同的处理方式, 对线性项采用Legendre谱方法隐式处理, 而对于非线性项部分采用谱配置法来显式处理, 可以是Legendre-Gauss-Lobatto(LGL)配置点, 也可以是Chebyshev-Gauss-Lobatto(CGL)配置点等. 这样, 数值方法既具有较好的稳定性又方便算法实施. 对非一致介质, 本文主要考虑规则介质交界面的情形, 并假设介电常数和磁导率是计算区域上的分片常数函数. 此方法自然地在介质交界处进行区域剖分, 基于合理的弱形式, 采用类似自然边界条件的方式处理交界面条件, 不需要使用附加的连接条件, 构建一致的数值格式, 并获得谱精度. 我们在时间方向采用LFCN格式离散, 该方法结合非线性部分的leap frog显隐格式和线性高阶导数部分的Crank-Nicolson隐式格式, 早期由文献[36]提出应用于Korteweg-de Vries 方程的Fourier谱方法, 并进行了线性化稳定性分析. 本文详细证明了非线性问题LFCN数值格式的稳定性和收敛性, 并得到$L^2$ -范数的最优误差估计. 数值算例中, 利用快速Legendre变换在Chebyshev点上计算非线性项, 提高算法效率. 数值结果显示了本文方法求解此类非线性问题的有效性和高精度, 同时证实了本文方法没有因为解的间断性而损失谱精度.

设$\Omega=I^x\times I^y=(0, 1)\times(0, 1),$ 考虑下面具有非线性电导率的二维Maxwell方程

其中${\textbf{E}}=(E_x(x, y, t), E_y(x, y, t))$, $H_z=H_z(x, y, t)$分别表示电场强度和磁场强度. $\epsilon$是介电常数, $\mu$是磁导率. 非线性项$\textbf{J}=(J_x,J_y)$表示电流密度, 其中$J_x=\sigma(|\textbf{E}|)E_x$, $J_y=\sigma(|\textbf{E}|)E_y$. $\sigma(\xi)$是光滑递增的一元函数. 考虑理想电导体PEC边界条件

和初始条件

非一致介质电磁场中的介电常数$\epsilon$和磁导率$\mu$在跨过不同的介质时产生了跳跃, 从而使解产生了间断或弱间断. 假设介质交界面在$x=\Gamma_1$和$y=\Gamma_2$处, 由Maxwell方程的积分形式, 可推导出如下交界面连接条件^[20]

其中$[\cdot]$表示在交界面处沿法向量的跳跃.

下面以上述交界面情况为例, 说明一下方程解是满足交界面条件(2.6)式的. 由(2.1)-(2.3)式知, 解$(E_x(t),E_y(t),H_z(t))$在整个区域$\Omega$上满足方程. 从而$E_x$在$y$方向是可导的, $E_y$在$x$方向是可导的, $H_z$在$x$和$y$方向都是可导的, 所以$[E_x]|_{y=\Gamma_2}=0,$$[E_y]|_{x=\Gamma_1}=0$以及$[H_z]|_{x=\Gamma_1}=[H_z]|_{y=\Gamma_2}=0$ 需要被满足. 进而由(2.1)和(2.2)式知, $[\epsilon E_x]|_{x=\Gamma_1}=0$和$[\epsilon E_y]|_{y=\Gamma_2}=0$需要被满足. 因此, 方程解是满足交界面条件(2.6)式的.

记$(\cdot,\cdot)_Q$和$\|\cdot\|_Q$是$L^2(Q)$空间的內积和范数, 并约定当$Q=\Omega$时省略下标$Q$. $H^{r}(\Omega)$$(r>0)$为通常的Sobolev空间, 其上的范数记为$\|\cdot\|_{H^{r}(\Omega)}$. 引入下列空间

二维非线性间断解Maxwell方程(2.1)-(2.4)的弱形式可写为: 求$(E_x(t),E_y(t),H_z(t))\in H^{0, 1}_{0}(\Omega)\times H^{1, 0}_{0}(\Omega)\times L^2(\Omega)$ 使得

这里假设介质被分成不同的矩形块, $\epsilon,\mu$在不同介质上是常数. 注意到在上面弱形式(2.7)和(2.8)式中, 使用了分部积分, 解$E_x$和$E_y$的正则性已经各向异性地被降低了. 因此(2.6)式中的连接条件$[\epsilon E_x]|_{x=\Gamma_1}=0$和$[\epsilon E_y]|_{y=\Gamma_2}=0$ 就类似有限元方法中自然边界条件一样以一种自然的方式进行了处理. 有关非线性Maxwell方程弱解的存在性可参见文献[9]. 我们考虑的解满足弱形式(2.7)-(2.9)式, 该弱形式有助于启发我们构造合理的逼近空间.

本节给出问题(2.1)-(2.5)的多区域Legendre-tau配置谱方法半离散格式, 并证明格式的稳定性和收敛性.

记

并设$\epsilon|_{\Omega_{lm}}$和$\mu|_{\Omega_{lm}}$是常数. 设$N^x_l,N^y_m$分别为$(a_{l-1},a_l),(b_{m-1},b_m)$上的积分节点数, 对$N^x_l, N^y_m~(1\le l\le L, 1\le m\le M)$, 定义

设${\Bbb P}_{\Bbb N}$是区间$I$上次数不超过$N$的多项式空间, 定义下列分片多项式空间

以及逼近空间

记参考区间$\hat I=(-1,1)$. 设$\hat x^l_i~(0\leq i\leq N^x_l)$, $\hat y^m_j~(0\leq j\leq N^y_m)$分别是参考区间$\hat I^x=(-1, 1)$, $\hat I^y=(-1, 1)$上LGL积分节点, $\hat{x}_i^{l,C}~(0\leq i\leq N^x_l)$, $\hat{y}_j^{m,C}~(0\leq j\leq N^y_m)$分别是参考区间$\hat I^x=(-1, 1)$, $\hat I^y=(-1, 1)$上CGL积分节点, 并定义

定义分片连续函数空间

建立函数在子区域$\Omega_{lm}$上与参考区域$\hat{\Omega}$上之间的对应关系

设插值算子$I_{{\mathbb N}}$可以是下列情形之一

i. LGL插值算子$I_{{\mathbb N}}^L:\tilde C(\bar{\Omega})\rightarrow {\Bbb P}_{\Bbb N}(\Omega)$, 满足

ii. CGL插值算子${I^C_N}: \tilde C(\bar{\Omega})\rightarrow {\Bbb P}_{\Bbb N}(\Omega)$, 满足

定义如下$L^2$ -正交投影算子

下面给出问题(2.1)-(2.5)的多区域Legendre-tau配置谱方法半离散格式: 求$(E_{x{\mathbb N}}(t),$$E_{y{\mathbb N}}(t),$$H_{z{\mathbb N}}(t))\in V_{\mathbb N}^{0, 1}\times V_{\mathbb N}^{1, 0}\times V_{\mathbb N}$使得

其中${\textbf{E}_{{\mathbb N}}(t)}=(E_{x{\mathbb N}}(t), E_{y{\mathbb N}}(t)).$ 初始值取为

在数值分析过程中, 考虑$I_{{\mathbb N}}$为LGL插值算子$I^L_{{\mathbb N}}$的情形. 为了方便起见, 仍采用记号$I_{{\mathbb N}}$. 对非线性项进行估计, 需要用到下面的引理.

${\bf 引理3.1}$ 如果$u\in {\Bbb P}_{N^x}(\hat{I}_x)\otimes{\Bbb P}_{N^y}(\hat{I}_y)$, 那么

对于多区域的情况, 如果$u\in {\Bbb P}_{{\mathbb N}}(\Omega),$有

${\bf证}$ 给出$u$的Legendre级数展开式

注意到Legendre多项式有性质$|L_l(x)|\leq 1(|x|\leq 1)$, $|L_m(y)|\leq 1(|y|\leq 1)$, 从而有

因此得到结论(3.5)式.

对于多区域情况, 注意到 $\|\hat{u}_{lm}\|_{L^2(\hat\Omega)}=\frac{2}{\sqrt{h^x_l h^y_m}}\big\|u|_{\Omega_{lm}}\big\|_{L^2({\Omega}_{lm})},$ 从而有

因此得到结论(3.6)式.

定义分裂Sobolev空间$\tilde H^r(\Omega)$,有

及其上范数

进一步可得到下面的逼近引理.

${\bf 引理3.2}$^[37-38] 对于任意$u\in \tilde{H}^r(\Omega)$,有

设$E_{x{\mathbb N}}$, $E_{y{\mathbb N}}$, $H_{z{\mathbb N}}$和右端分别有扰动$\tilde{E}_{x{\mathbb N}}$, $\tilde{E}_{y{\mathbb N}}$, $\tilde{H}_{z{\mathbb N}}$, $\tilde{f_1}$, $\tilde{f_2}$, $\tilde{f_3}$, 考虑如下扰动方程. 对$0< t \leq T$,有

其中, $\tilde{\textbf{E}}_{{\mathbb N}}=(\tilde{E}_{x{\mathbb N}}, \tilde{E}_{y{\mathbb N}}),$

在扰动方程(3.9)-(3.11)中, 令$v_x=\tilde{E}_{x{\mathbb N}}$, $v_y=\tilde{E}_{y{\mathbb N}}$, $v_z=\tilde{H}_{z{\mathbb N}}$, 可得

从而有

这里$I_{{\mathbb N}}$为LGL插值算子.

下面来估计非线性项. 定义如下离散范数

其中$w^l_i=\frac12 h^x_l\hat w^x_i,$$w^m_j=\frac12 h^y_m\hat w^y_j,$$\hat w^x_i$和$\hat w^y_j$分别为参考区间$\hat I^x,\hat I^y$上LGL求积公式的系数. 根据文献[39中(5.9.4)]式可知, 存在常数$C$使得

因此有

类似地, 有 $\big|(I_{{\mathbb N}}\tilde{J_y},\tilde{E}_{y{\mathbb N}})\big|\le \|\tilde{J_y}\|_{{\mathbb N}}\|\tilde{E}_{y{\mathbb N}}\|.$

现估计非线性项$\|\tilde{J}_x\|_{{\mathbb N}}$和$\|\tilde{J_y}\|_{{\mathbb N}}$. 假设逼近解是有界的, 并记

引入记号

令$B$是正常数, 对给定的$0<t_1\leq T$, 设有

为了方便表达, 在下面的证明过程中, 我们将会省略记号$s$. 当$0<s\le t_1$时, 对于非线性项$\tilde{J}_x$, 有

由微分中值定理得

其中$\zeta$介于$|\textbf{E}_{{\mathbb N}}+\tilde{\textbf{E}}_{{\mathbb N}}|$和$|\textbf{E}_{{\mathbb N}}|$ 之间. 利用(3.14)和(3.15)式, 由(3.16)和(3.17)式得到

其中$C$是依赖于$E_M, B, C^1_\sigma(E_M,B), C^0_\sigma(E_M), \frac {1}{\sqrt{\epsilon}}$的正常数. 类似地, 可得 $\|\tilde{J_y}\|_{{\mathbb N}} \leq C\|\sqrt{\epsilon}\tilde{\textbf{E}}_{{\mathbb N}}\|.$

因此, 对于$0<s\le t_1,$有

对$s$从$0$到$t(t\leq t_1)$积分得

令

综上可得, 在条件(3.15)式满足时

记$\underline{\epsilon}=\min \limits_{1\leq l \leq L,1\leq m \leq M}{\{\epsilon|_{\Omega_{lm}}\}}$, 有如下稳定性定理.

${\bf 定理3.1}$ 如果$\rho(T)\le \underline{\epsilon} B^2{\rm e}^{-CT}\big/\bar{N}^2$, 则有

${\bf证}$ 首先用反证法证明下面结论

否则存在$t_1$, $0<t_{1}<T$, 使得

由(3.21)式, 并根据Gronwall不等式, 得

根据引理3.1, 有

这与(3.24)式矛盾, 从而(3.23)式成立. 因此, $t_1=T$时, (3.21)式成立. 由Gronwall不等式, 即证得稳定性定理.

下面先给出误差估计时要用到的投影算子. 定义$H^1$ -投影算子$\hat P^{1, x}_N: H^1(\hat I^x)\rightarrow {\Bbb P}_{N}(\hat I^x)$, 满足

注意到$H^1$ -投影算子$\hat P^{1, x}_N$保持了边界条件, 并且满足下面的性质

类似地, 定义$\hat P^y_N$ 和 $\hat P^{1, y}_N$. 进而, 可导出下列投影算子

由算子$\hat P^{1, x}_{N^x_l}$的性质和(3.27)式, 对$u\in H^{1, 0}_{0}$, 有$P^{1, 0}_{{\mathbb N}} u \in V_{\mathbb N}^{1, 0}$, 并且对任意的$v \in V_{\mathbb N}$, 有

类似地

现给出下面的逼近引理.

${\bf 引理3.3}$^[27,40-41] 对于任意$u\in \tilde{H}^r(\Omega)$,有

取比较函数 $E_{x*}=P^{0, 1}_{{\mathbb N}} E_x \in V_{\mathbb N}^{0, 1}, E_{y*}=P^{1, 0}_{{\mathbb N}} E_y\in V_{\mathbb N}^{1, 0}$, $H_{z*}=P_{{\mathbb N}} H_z\in V_{\mathbb N}$, 可得如下截断误差方程: 对任意的$v_x \in V_{\mathbb N}^{0, 1}, v_y \in V_{\mathbb N}^{1, 0},v_z \in V_{\mathbb N},$有

其中, $\textbf{E}_{*}(t)=(E_{x*}(t), E_{y*}(t)),$

令$e_{x{\mathbb N}}(t)=E_{x*}(t)-E_{x{\mathbb N}}(t), e_{y{\mathbb N}}(t)=E_{y*}(t)-E_{y{\mathbb N}}(t), e_{z{\mathbb N}}(t)=H_{z*}(t)-H_{z{\mathbb N}}(t)$, 可得如下误差方程

其中, $\textbf{e}_{{\mathbb N}}(t)=(e_{x{\mathbb N}}(t), e_{y{\mathbb N}}(t)),$

类似前面半离散格式稳定性的讨论过程, 有如下收敛性定理.

${\bf 定理3.2}$ 设$(E_x, E_y, H_z)$是问题(2.1)-(2.5)的精确解, $(E_{x{\mathbb N}}, E_{y{\mathbb N}}, H_{z{\mathbb N}})$是半离散格式(3.1)-(3.4)的数值解. 设$r>2$, $E_x, E_y\in H^1(0, T; \tilde H^r(\Omega))$, $H_z\in L^2(0, T; \tilde H^r(\Omega))$, 则存在依赖解的正则性的常数$C$, 使得

${\bf证}$ 记

类似定理3.1的证明过程, 如果$\eta(T)\le \epsilon B^2{\rm e}^{-CT}/\bar{N}^2$, 则 $G(t)\le \eta(t){\rm e}^{Ct},~\forall~ 0<t\le T.$ 这里需要估计初始误差$G(0)$和截断误差$\|\sqrt{\epsilon}{ f_1}(s)\|^2 +\|\sqrt{\epsilon}{ f_2}(s)\|^2$, 验证误差$\eta(t)$满足定理3.1的条件. 对于初始误差, 由引理3.2得

接下来估计截断误差$\|\sqrt{\epsilon}{f_1}(s)\|^2 +\|\sqrt{\epsilon}{ f_2}(s)\|^2$. 简记$f_1(s)$为

对于$g_2$,有

记$E_{M}:=\|\textbf{E}\|_{L^{\infty}(\Omega)}$, $E_{*M}:=\|\textbf{E}_{*}\|_{L^{\infty}(\Omega)}$, 类似稳定性分析, 则有

进而根据引理3.3得

因此, $\int^t_0 \|\sqrt{\epsilon}{ f_1}(s)\|^2 {\rm d}s\leq C\hbar^{2r}.$ 类似地, 可得$\int^t_0 \|\sqrt{\epsilon}{f_2}(s)\|^2 {\rm d}s\leq C\hbar^{2r}.$ 因此, 当$r>2$且$\hbar$足够小时, 初始误差和截断误差$\eta(t)$满足定理3.1条件. 利用半离散的稳定性结论(3.22)式, 可得下面误差估计

再由引理3.2和引理3.3, 以及三角不等式, 就证得收敛性定理.

本节给出问题(2.1)-(2.5)的LFCN多区域Legendre-tau配置谱方法全离散格式, 并证明全离散格式的稳定性和收敛性.

设$\tau$是时间步长, $t_k=k\tau\,\,(k =0,1,\cdots,n_T; T = n_T \tau)$, 并定义

下面给出问题(2.1)-(2.5)的LFCN多区域Legendre-tau配置(LFCN-MLTC) 谱方法的三层全离散格式: 求$E_{x{\mathbb N}}^{k+1}\in V_{\mathbb N}^{0, 1}$, $E_{y{\mathbb N}}^{k+1}\in V_{\mathbb N}^{1, 0}$, $H_{z{\mathbb N}}^{k+1}\in V_{\mathbb N}$ 使得对于$1\leq k\leq n_T-1$有

初始值取为

假设解和右端分别有扰动$\tilde{E}^k_{x{\mathbb N}}$, $\tilde{E}^k_{y{\mathbb N}}$, $\tilde{H}^k_{z{\mathbb N}}$和$\tilde{f}_1^k$, $\tilde{f}_2^k$, $\tilde{f}_3^k$, 考虑LFCN-MLTC谱方法全离散格式的扰动方程: 对$1\leq k\leq n_{T}-1$,有

其中

在(4.5)-(4.7)式中, 分别取

整理可得

这里$I_{{\mathbb N}}$为LGL插值算子.

下面估计(4.8)式中的非线性项. 类似(3.13)式可得 $|(I_{{\mathbb N}}\tilde{J}_x^k,\tilde{E}_{x{\mathbb N}}^k)|\le \|\tilde{J}_x^k\|_{{\mathbb N}}\|\tilde{E}_{x{\mathbb N}}^k\|.$ 现估计非线性项$\|\tilde{J}_x^k\|_{{\mathbb N}}$. 假设逼近解是有界的, 记 $E_M=\max\limits_{0\leq k\leq n_{T}}\|{\textbf{E}}_{{\mathbb N}}^k\|_{L^{\infty}(\Omega)}.$ 令$B$是正常数, 对给定的$2\leq n\leq n_T$, 设有

类似半离散非线性项的估计(3.18)式, 当$1\leq k\le n-1$时, 对于非线性项$\tilde{J}_x^k$ 有如下估计

同理可得

对于右端误差$\tilde{f}^k_1, \tilde{f}^k_2, \tilde{f}^k_3$, 有

结合(4.10), (4.11)和(4.12)式, 由(4.8)式可得, 对于$1\leq k\le n-1,$有

上式关于$k$从$1$到$n-1(2\le n\le n_T)$求和, 并设$C\tau\le \frac12$, 整理可得

令

并且定义$\rho^1:=C(K^0+K^1)$.

综上可得, 对于$2\leq n\leq n_T$, 当条件(4.9)式满足时, 有

从而有如下稳定性定理.

${\bf定理4.1}$ 如果$\rho^{n_T}\le \underline{\epsilon} B^2{\rm e}^{-CT}/\bar{N}^2$, 则有

${\bf证}$ 采用归纳法证明定理结论.

首先对于$n=1$时, 有下面结论, $K^n\leq \rho^n{\rm e}^{C{n\tau}}.$ 对给定的$2\leq n \leq n_T$, 假设当$k\leq n-1$, 结论成立, 即

所以对$1\leq k \leq n-1$, 由定理条件以及引理3.1和(4.17)式, 类似(3.26)式, 有 $\|\tilde{{\textbf{E}}}_{{\mathbb N}}^k\|_{L^{\infty}(\Omega)}< B.$ 因此条件(4.9)式满足. 从而由(4.15)式和归纳假设(4.17)式得到

因此, $k=n$时结论成立. 从而证得定理.

本节讨论LFCN-MLTC谱方法全离散格式(4.1)-(4.3)式的收敛性, 并证明对于Legendre插值情况下的$L^2$ -范数最优阶误差估计.

引入合适的比较函数

记

给出下列截断误差方程

其中

由投影算子$P^{0, 1}_{{\mathbb N}},P^{1, 0}_{{\mathbb N}}$的性质(3.27)-(3.29)式以及$P_{\mathbb N}$的性质, 我们可以得到(4.18)式中的$f_1^k, f_2^k, f_3^k$.

令 $e^k_x=E^k_{x*}-E^k_{x{\mathbb N}}, e^k_y=E^k_{y*}-E^k_{y{\mathbb N}}, e^k_z=H^k_{z*}-H^k_{z{\mathbb N}}$, 得下列误差方程

其中, 记$\textbf{e}^k=(e_{x}^k, e_{y}^k),$

由稳定性结论, 可以得到下面收敛性定理.

${\bf定理4.2}$ 设$N, n, r$是正整数, 且$r> 2$, $(E_x, E_y, H_z)$是(2.1)-(2.5)的解, $(E^{n}_{x{\mathbb N}}, E^{n}_{y{\mathbb N}}, H^{n}_{z{\mathbb N}})$是全离散格式(4.1)-(4.4)的解, $E_x, E_y\in H^1(0, T; \tilde H^r(\Omega))\cap H^3(0, T;L^2(\Omega))$, $H_z\in C([T];$$\tilde H^r(\Omega))\cap H^3(0, T;L^2(\Omega))$. 当$\tau$和$\hbar$足够小, 则存在常数$C$使得

${\bf证}$ 对误差方程(4.21)-(4.23)应用稳定性定理4.1. 先验证误差方程中的初始和右端误差满足定理4.1中的条件.

对于$\|e_x^0\|^2+\|e_y^0\|^2+\|e_z^0\|^2$, 由(3.38)式给出如下估计

对于$\Vert{\sqrt{\epsilon} e^1_x}\Vert+\Vert{\sqrt{\epsilon} e^1_y}\Vert+\Vert{\sqrt{\mu} e^1_z}\Vert$, 由Taylor展式,

类似地, 可得

接着估计右端误差. 对$f^k_1$, 有下面估计

利用分部积分, 有

进一步, 可得

从而结合引理3.3有

类似(3.39)式及上式的推导, 可估计得到

因此

同理可得

因此, 对任意的$1\leq n \leq n_T$, 有

故当$r>2$且$\tau$和$\hbar$足够小时, 误差方程(4.21)-(4.23)的初始和右端误差满足定理4.1条件. 利用全离散稳定性结论(4.16)式, 并代入上述对初始和右端误差的估计结果, 有如下误差估计

再结合三角不等式, 并利用引理3.2和引理3.3, 即证得收敛性结论.

在本节中, 给出充分的数值算例来验证理论分析结果, 证实LFCN-MLTC方法求解非线性问题的有效性和高精度. 在算例中, 对非线性项采用Chebyshev插值, 选取CGL插值节点, 结合快速Legendre变换(FLT), 提高计算效率.

为了验证算法精度和有效性, 定义移位CGL点$(x^{l,C}_i,y^{m,C}_j)$上$t$时刻的最大模误差为

${\bf例5.1}$ 非线性电导率取为$\sigma(|\textbf{E}|)=|\textbf{E}|^2-|\textbf{E}|^4$, 则非线性项

取解$(E_x, E_y, H_z)$如下

其中$(x, y)\in[0,1]\times[0,1]$, $t\in(0, T]$, 参数$\omega, k_x, k_y$满足关系$\omega^2=\frac{k_x^2+k_y^2}{\mu\epsilon}$. 右端源项可由方程得到.

本算例中, 考虑非一致介质中具有非线性电导率的Maxwell方程间断解问题. 此时问题的解在交界面处会出现间断或弱间断. 这里假设介电常数$\epsilon$在$x=\frac12$产生跳跃, $\Omega$被剖分成$\Omega_1=(0, \frac12)\times(0, 1)$和$\Omega_2=(\frac12, 1)\times(0, 1)$. 问题的解在每个子区域上是分片光滑的.

设跳跃的介电常数$\epsilon$在$\Omega$上是分段常数

$k_x$在$\Omega$上也是分段常数

取$\mu=1$, $k_y=8$. $\Omega=\left(\frac{k_x^2+k_y^2}{\mu\epsilon}\right)^{\frac12}$在$\Omega$上也是分段常数.

$\bar{\Omega}$进行如下剖分

记${{\mathbb N}}^x=(N^x_{1}, N^x_{2}), {{\mathbb N}}^y=N^y_{1}$. 现检验本数值方法在时间和空间上的精度. 取$\tau=10^{-6}$时, 令${{\mathbb N}}^x$和${{\mathbb N}}^y$从$(12,24)$和$24$增加到$(16,32)$和$32$; 取${{\mathbb N}}^x=(16,32), {{\mathbb N}}^y=32$时, 让$\tau$从$10^{-4}$减少到$10^{-6}$. 表1和表2分别给出了$t=1$时刻相应的误差. 数值结果表明本方法在空间上具有高精度, 时间上是二阶精度, 验证了LFCN-MLTC法在求解此非线性间断解问题的有效性和高精度, 同时证实了本文方法没有因为解的间断性而损失谱精度.

${\bf例5.2}$$\Omega=(0,1)\times(0,1)$, 非线性电导率取为$\sigma(|\textbf{E}|)=|\textbf{E}|^2-|\textbf{E}|^4$. 如同文献[12], 解$(E_x, E_y, H_z)$取为

右端项可由方程得到.

现来验证LFCN-MLTC法的误差精度. 计算时, 取$\omega=1$, $\bar{\Omega}$被剖分成四个子单元

并记${{\mathbb N}}^x=(N^x_{1}, N^x_{2}), {{\mathbb N}}^y=(N^y_{1}, N^y_{2})$. 取$\tau=1.0\times10^{-5}$, 让${{\mathbb N}}^x={{\mathbb N}}^y$从$(6,6)$增加到$(12,12)$, 表3给出了LFCN-MLTC法在$t=1$时刻的最大模误差. 另外, 为了减少时间方向误差对空间误差的影响, 表中也给出了取较小的$\tau=1.0\times10^{-6}$, 在${{\mathbb N}}^x={{\mathbb N}}^y=(12,12)$时的误差结果. 数值结果表明了LFCN-MLTC法求解此类非线性问题的有效性和高精度.

${\bf例5.3}$ 设$\Omega=(0,1)\times(0,1)$, 非线性电导率取为$\sigma(|\textbf{E}|)=|\textbf{E}|^{1-\alpha}$, 可分别取$\alpha=0.3,$$0.5, 0.6, 0.8$. 如同文献[12], 解$(E_x, E_y, H_z)$取为

右端项可由方程得到.

同样地, $\bar{\Omega}$被剖分成四个子单元

并记${{\mathbb N}}^x=(N^x_{1}, N^x_{2}), {{\mathbb N}}^y=(N^y_{1}, N^y_{2})$. 取$\alpha=0.5$, $\tau=1.0\times10^{-5}$, 让${{\mathbb N}}^x={{\mathbb N}}^y$从$(6,6)$增加到$(10,10)$, 表4给出了$t=1$时刻的最大模误差. 表5给出了$\tau=1.0\times10^{-5}$, ${{\mathbb N}}^x={{\mathbb N}}^y=(10,10)$时, 分别取$\alpha=0.3, 0.5, 0.6, 0.8$, LFCN-MLTC法在$t=1$时刻的最大模误差. 数值结果表明本文方法求解此类非线性问题的可靠性.

${\bf例5.4}$ 类似文献[12], 考虑$\textbf{L}$ -型区域$\bar{\Omega}=[0,1]^2\setminus[0.5,1]^2$, 非线性电导率取为 $\sigma(|\textbf{E}|)=|\textbf{E}|^{\alpha}$. 解$(E_x, E_y, H_z)$取为

右端项为

现验证本文方法的误差精度. 计算时, $\textbf{L}-$型区域$\bar{\Omega}$被剖分成三个子单元

取$\epsilon=\mu=1, \omega=1$, 并记${{\mathbb N}}^x=(N^x_{1}, N^x_{2}), {{\mathbb N}}^y=(N^y_{1}, N^y_{2})$. 取$\alpha=0.3$, $\tau=1.0\times10^{-5}$, 让${{\mathbb N}}^x={{\mathbb N}}^y$从$(8,8)$增加到$(16,16)$, 表6给出了LFCN-MLTC法在$t=1$时刻的最大模误差. 取$\alpha=0.3, 0.5, 0.6, 0.8$, 表7给出了$\tau=1.0\times10^{-5}$, ${{\mathbb N}}^x={{\mathbb N}}^y=(16,16)$ 时, LFCN-MLTC法在$t=1$时刻的最大模误差. 数值结果表明本文方法求解$\textbf{L}$ -型区域下的非线性Maxwell方程问题的高精度和有效性.

本文讨论了具有非线性电导率的非一致介质的二维Maxwell方程的数值求解方法, 提出了LFCN多区域Legendre-tau配置谱方法. 此数值方法采用显隐数值格式, 对线性和非线性部分分别采用不同的处理方法. 线性项采用Legendre谱方法隐式处理, 非线性项采用谱配置法来显式处理. 此方法采用不同次的多项式空间逼近电磁场, 基于合理的弱形式, 逼近格式不需要使用额外附加的交界面连接条件, 构建一致的数值格式, 简化计算, 同时很好地保持了谱精度. 时间方向上采用分裂的LFCN三层格式. 非线性项放在中间已知层来处理. 给出半离散和全离散格式, 并证明了数值格式的稳定性和收敛性, 对于Legendre插值的情况得到了$L^2$ -范数意义下最优阶的误差估计. 利用快速Legendre变换在Chebyshev点上计算非线性项, 提高计算效率. 数值算例验证了理论分析结果, 显示了本文数值方法求解此类非线性非一致介质问题的有效性和高精度.