Energy analysis of the ADI-FDTD method with fourth-order accuracy in time for Maxwell's equations

1.   Introduction

Maxwell's equations are very important partial differential equations and they play a significant role in electromagnetic theory. Numerical solutions for Maxwell's equations are extensive in science and engineering, for example, for radio-frequencies, antennas, microwaves, wireless engineering and the design of CPUs in microelectronics. There are many efficient numerical methods for solving Maxwell's equations, such as the finite-difference time-domain (FDTD) method ^[1,2], the finite element method ^[3], the weak Galerkin finite element method ^[4] and so on.

The FDTD method (also called Yee$'$s scheme) was first introduced by Yee ^[1] in 1966, and later applied to many problems in computational electromagnetics ^{[5,6,7,8,9,10,11,12,13,14,15,16,17,18,19]}. However, the traditional FDTD method is conditionally stable and must satisfy the Courant-Friedrichs-Lewy (CFL) stability condition; see ^[5,20,21,22]. For the 2-D Maxwell's problem, the CFL condition is

where $c$ is the maximum of the wave velocity, i.e., the speed of light in vacuum, $\Delta t$ is the time step size and $\Delta x$ and $\Delta y$ are the spatial step sizes. This implies that the maximum time step size is limited by the minimum spatial size, and that computation of the FDTD method needs much CPU time when the spatial step sizes are small.

Many numerical methods have been proposed for Maxwell's equations to get rid of this restriction. One of these methods is the alternating direction implicit FDTD method, which was proposed in 1999 by Namiki ^[2] and which zheng et al. ^[12] proved to be unconditionally stable. The ADI-FDTD scheme consists of two stages and can be solved directly by using the Thomas algorithm, since each stage includes tri-diagonal systems of linear equations. By truncation error analysis, it was found that this method has second-order accuracy in both time and space. The ADI-FDTD method is quickly applied to various electromagnetic computations and the ADI-FDTD methods with high-order accuracy were developed. Based on Yoshida's work ^[9], Tan and Ding ^[23] proposed a new ADI-FDTD method with fourth-order accuracy in time for 2D-wave propagation in a lossless, isotropic medium. The new ADI-FDTD method was proved to be unconditionally stable by the Fourier method, i.e., this method is free from the constraints of the CFL conditions. However, energy analysis of the new method in terms of energy preservation and convergence are not available.

Analysis of the energy of the methods and construction of energy-preserving methods are popular and valuable. As Li and Vu-Quoc said in ^[10], to some extent, the ability to preserve some invariant properties of the original differential equation is a criterion to judge the success of a numerical simulation. Therefore, it is important to develop energy-preserving algorithms and analyze the energy of existing algorithms for solving Maxwell's equations. It is well known that the FDTD method is a popular numerical method in computational electromagnetics, but it is not unconditionally stable and fails to satisfy the energy conservation property. In order to make up for these shortcomings, various energy-preserving methods have been applied to solve Maxwell's equations based on the ADI and splitting techniques. For example, a new splitting FDTD scheme for Maxwell's equations was proposed in ^[24]. It proved that the new scheme was unconditionally stable and energy-preserving. Furthermore, Chen et al. ^[25] proposed the symmetric energy-conserved splitting FDTD scheme (i.e., symmetric EC-S-FDTD). It was proved that this symmetric EC-S-FDTD method is unconditionally stable and has second-order convergence in both time and space. In ^[26], a symmetric EC-S-FDTD method for Maxwell's equations in negative index metamaterials was proposed. We note that these works are limited to cases with homogeneous boundary conditions.

In addition, a proper summation by parts (SBP) formula was found for the approximate derivative in ^[27]. A "simultaneous approximation term" was used to treat the inhomogeneous boundary conditions. Appel$\ddot{o}$ and Bokil ^[28] presented the construction of novel SBP-FDTD methods for the numerical discretization of the Maxwell-Duffing models and derived energy estimates for the semi-discrete methods that are analogous to the continuous energy estimates. Boundary and interface conditions were handled by the simultaneous approximation technique. In ^[29], a new framework to construct time-stable finite-difference schemes was proposed for hyperbolic systems base on the application of strong boundary conditions. Sufficient conditions for strong time stability and conservation were derived for the linear advection equation and coupled system of hyperbolic equations using the energy method. For other energy-preserving methods of Maxwell's equations, please see the relevant references ^{[30,31,32,33]}. In quantum mechanics, the conservative methods are also favored by researchers; see ^{[34,35,36,37,38]}.

In this paper, we consider the following 2-D transverse electric (TE) problem in a lossless medium and without sources and charges:

with $(x, y)\in\Omega = [0, a]\times[0, b]$ and $t\in (0, T]$, where $\varepsilon$ is the electric permittivity of the medium, $\mu$ is the magnetic permeability and $\textbf{E} = (E_x(x, y, t), E_y(x, y, t))$ and $H_z = H_z(x, y, t)$ are the electric field and the magnetic field, respectively.

Motivated by Tan and Ding ^[23], we present a ADI-FDTD scheme with fourth-order accuracy in time for the TE problems (1.1)–(1.3), and an analysis of the energy of the scheme. We use the energy method to derive the numerical energy identity of this method. It is found that the method is approximately energy-preserving. In addition, numerical experiments are provided to verify the numerical performance of the ADI-FDTD scheme with fourth-order accuracy in time.

The remaining part of the paper is organized as follows. In Section 2, we introduce some preliminaries and notations. In Section 3, we give the detailed energy analysis of the scheme. In Section 4, we present some numerical results to illustrate the accuracy, convergence and energy identity of the scheme.

2.   Preliminaries and notations

We consider the perfectly electric conductor (PEC) boundary condition as follows:

where $\partial\Omega$ denotes the boundary of $\Omega$ and $\vec{\bf{n}}$ is the outward normal vector on $\partial\Omega$. The initial conditions are

It is well known that for suitably smooth data, the problem has a unique solution for all time; see^[8]. Here, we assume that the solution of the Maxwell's equations (1.1)–(1.3) has the following regularity property:

To simplify the notations we only consider the case of constant coefficients which are independent of $x$ and $y$. But the proposed method is valid for the case of variable coefficients with $\varepsilon = \varepsilon(x, y)$ and $\mu = \mu(x, y)$.

For the problems in a lossless medium, Poynting's theorem states that electromagnetic energy stays constant for all time. The following lemma is an integral form of the Poynting's theorem.

Lemma 2.1. ^[24]If E and H are the solutions of the Maxwell's equations (1.1)–((1.3) in a lossless medium and satisfy the PECboundary condition (2.1), then for any $t\geq 0$, it holds that

Consider the space domain $\Omega = [0, a]\times[0, b]$ and time interval [0, T]. Let $\Delta x > 0$ and $\Delta y > 0$ be space steps and $\Delta t > 0$ be the time step; define

where $0\leq i \leq I, 0\leq j \leq J$ and $0\leq n \leq N$; $I$, $J$ and $N$ are integers. Define $(x_\alpha, y_\beta, t^\theta) = (\alpha\Delta x, \beta\Delta y, \theta\Delta t)$, where $\alpha$ is either $i$ or $i+\frac{1}{2}$, $\beta$ is either $j$ or $j+\frac{1}{2}$ and $\theta$ is either $n$ or $n+\frac{1}{2}$. For a function $F(x, y, t)$, define

and $\delta_u \delta_v F_{\alpha, \beta}^m = \delta_u(\delta_v F_{\alpha, \beta}^m)$ where $u, v\in\{x, y\}.$

Denote by $E_{x_{i+\frac{1}{2}, j}}^m$, $E_{y_{i, j+\frac{1}{2}}}^m$ and $H_{z_{i+\frac{1}{2}, j+\frac{1}{2}}}^m$ the approximations of the electric field $E_x(i+\frac{1}{2}, j, t^m)$, $E_y(i, j+\frac{1}{2}, t^m)$ and $H_z(i+\frac{1}{2}, j+\frac{1}{2}, t^m)$, respectively. For a grid function $F_{\alpha, \beta}$, where $\alpha = i, i+\frac{1}{2}$ and $\beta = j, j+\frac{1}{2}$, define the discrete $L^2$ norms:

For a grid function $\textbf{F} = (U, V)$ over the mesh for the electric field, the $L^2$ norms of the vector-valued function is defined as follows:

Now we give the ADI-FDTD scheme with fourth-order accuracy in time for the TE problems (1.1)–(1.3). Let

and

and

The ADI-FDTD method consists of two steps:

Step 1:

Step 2:

where $I$ is a $3\times 3$ identity matrix and $u^{n}$ and $u^{n+1}$ are the fields at integer time steps. We know that this method is second-order accurate in time. Motivated by the method in ^[23], we design the new numerical scheme for the TE problems (1.1)–(1.3) as

where the coefficients $\alpha_{0}$ and $\alpha_{1}$ are given by

From ^[9], we know that the coefficients are determined by

It is easy to prove that the schemes described by (2.16)–(2.21) are unconditionally stable and have fourth-order accuracy in time. In order to analyze the numerical energy of the schemes given by (2.16)–(2.21), we rewrite them as follows:

Stage 1:

Stage 2:

Stage 3:

Stage 4:

Stage 5:

Stage 6:

The boundary and initial conditions are given by

and

3.   Energy analysis

The goal of this section is to derive the numerical energy identity of the schemes described by (2.22)–(2.27) and prove that this scheme is unconditionally stable. To this end, we fist prove the following lemma.

Lemma 3.1. Let ${E_{y_{i, j+\frac{1}{2}}}^{n}}$ and ${H_{z_{i+\frac{1}{2}, j+\frac{1}{2}}}^{n}}$ be the grid function in the ADI-FDTD schemes (2.22)–(2.27) that satisfies the boundary condition (2.1). Then, it holds that

Proof. Using the definition of the operator $\delta_x$, we have

The right side term of the above equation becomes

This completes the proof of Lemma 3.1.

Next, we use Lemma 3.1 to derive the numerical energy identity for the schemes (2.22)–(2.27).

Theorem 3.2. Let $\mathit{\boldsymbol{E}}^n = {(E_{x_{i+\frac{1}{2}, j}}^{n}, E_{y_{i, j+\frac{1}{2}}}^{n})}$ and $H_{z_{i+\frac{1}{2}, j+\frac{1}{2}}}^{n}$ be the solutionsof the schemes (2.22)–(2.27); then, it holds that

where

Proof. Square both sides of the equations in Stage 1 (2.22), so it is easy to get

We sum over all the terms over all $i$ and $j$, and then add these equations together, yielding

and

and

For the above equations given by (3.3)–(3.5), using Lemma 3.1 and the PEC boundary condition, we have

In a similar way, from the equations in (2.23), we can obtain that

and

and

We take the sum of (3.7)–(3.9), and get that

Similar to the derivation of (3.6), from (2.23) we get

Similarly, from (2.24)–(2.27), we also obtain the following identities:

and

and

Take the sum of (3.6)–(3.14), and then we have

It shows that the ADI-FDTD schemes (2.22)–(2.27) are approximately energy-preserving. This proof is completed.

4.   Numerical experiments

In this section, we provide some numerical results to verify the numerical performance of the ADI-FDTD schemes (2.16)–(2.21) for the TE models (1.1)–(1.3).

4.1.   Example 1

Set $\Omega = [0, 1]\times[0, 1]$ surrounded by a perfect conductor and consider the TE models (1.1)–(1.3) in a lossless medium with normalized electric permittivity and magnetic permeability, i.e., $\varepsilon = 1$ and $\mu = 1$. The exact solutions of the problems (1.1)–(1.3) are

The drive routines were written in Matlab, and the computation was run using a 2.20 GHz PC with 8 GB RAM and a Windows 10 operating system.

By calculation and simulation, we obtained some numerical results, as shown in Figures 1–9. The parameters were $\Delta t = 1/40$, $\Delta x$ = $\Delta y$ = $\Delta t^2$ and $T = 1$.

In Figures 1 and 2, we provide the surfaces between the exact solution $E_x$ and the numerical solution $E^n_x$ with $T = 1$. From the figures we can see clearly that the numerical solution behavior can approximate well the exact solution.

In Figure 3, the contour of the error $E_x-E^n_x$ is given. We can see clearly that the maximums absolute values of errors were close to $8\times 10^{-6}$, which means that the new scheme is effective for the problems (1.1)–(1.3).

In addition, we provide the surfaces between the exact solution $E_y$ and the numerical solution $E^n_y$ in Figures 4 and 5, as well as the contour of the error $E_y-E^n_y$ in Figure 6. One can see that the maximums of absolute values of errors were close to $2.3\times 10^{-5}$ from Figure 6. From these figures, we can see clearly that the numerical solution $E^n_y$ behavior can approximate well the exact solution $E_y$.

Similarly, Figures 7 and 8 show the surfaces between the exact solution $H_z$ and the numerical solution $H^n_z$. The contour of the error $H_z-H^n_z$ is presented in Figure 9; we can see that the maximums of absolute values of errors were close to $1.25\times 10^{-5}$. From these results, it is clear that the numerical solution $H^n_z$ behavior can approximate well the exact solution $H_z$.

In the following, we continue to describe some experiments to test the stability, convergence and energy conservation ability of the ADI-FDTD schemes (2.16)–(2.21). Let $ErrE$ and $ErrH$ be defined by

where $E_x(i+\frac{1}{2}, j, t_n), E_y(i, j+\frac{1}{2}, t_n)$ and $H_z(i+\frac{1}{2}, j+\frac{1}{2}, t_n)$ are the exact solutions of the problems (1.1)–(1.3) and $E^n_{x_{i+\frac{1}{2}, j}}, E^n_{y_{i, j+\frac{1}{2}}}, H^n_{z_{i+\frac{1}{2}, j+\frac{1}{2}}}$ are the solutions of the ADI-FDTD methods (2.16)–(2.21) for $n\geq 0$. $I_{n+1}^2$ is the left side of the identity (3.1) and $I_{n}^2$ is the right side of the identity (3.1). $E_1$ is the energy difference between the $n+1$ and $n$ levels. $E_2$ is the difference between the two sides of the Eq (3.1).

In calculation, we take $\Delta x = \Delta y = (\Delta t)^2$, $T = 1$ and $T = 2$. Tables 1 and 2 give the errors of the numerical solution of (1.1)–(1.3) as computed by the ADI-FDTD methods described by (2.16)–(2.21) in the discrete $L^2$ norms and given convergence rates in different time step sizes $\Delta t = T/N$.

Tables 3 and 4 give the energy difference between the level $n$ and level $n+1$ in the discrete $L^2$ norm and the difference between the two sides of the identity given by (3.1).

Table 5 shows the results for $I_n-I_0$, $E_1$ and $E_2$ for the spatial steps $\Delta x = \Delta y = 0.01$ and time step of $\Delta t = 0.001$.

From the numerical results in Tables 1–5, we get the following observations:

● Tables 1 and 2 show the errors and convergence rates with $T = 1$ and $T = 2$. We can see that the ADI-FDTD schemes (2.22)–(2.27) are efficient and have fourth-order accuracy in time.

● Tables 3 and 4 show the energy error between layer $n$ and $n+1$ and the difference between the left and right sides of the equation (3.1) with $T = 1$, $T = 2$. From the tables, we can see that $E_1$ and $E_2$ tend to zero, which verifies the theoretical results in Theorem 3.2.

● Table 5 shows the values of $I_n-I_0$, $E_1$ and $E_2$ at $n = 100,200,400,800, 1000, 2000$ and 4000. From the second line we can see that the ADI-FDTD schemes (2.22)–(2.27) are stable and approximately energy-preserving. The third line shows the energy error between layer $n$ and $n+1$; we can see that the values of $E_1$ were close to $5.5\times 10^{-14}$. The fourth line shows that the difference between the left and right sides of the equation defined by (3.1) was close to $10^{-15}$, which verifies Theorem 3.2.

4.2.   Example 2

Set $\Omega = [0, 1]\times[0, 1]$, $k_x = 1$, $k_y = 2$, $\varepsilon = \mu = 1$; the exact solutions of the problems (1.1)–(1.3) are

This example is presented for studying error estimates of the schemes (2.16)–(2.21) with $k_x\neq k_y$. The parameters were $\Delta t = 1/40$, $\Delta x$ = $\Delta y$ = $\Delta t^2$ and $T = 1$. All calculation results are shown in Figures 10–18.

In Figures 10 and 11, we provide the surfaces between the exact solution $E_x$ and the numerical solution $E^n_x$ with $t = 1$. From the figures we can see clearly that the numerical solution behavior can approximate well the exact solution.

The contour of the error $E_x-E^n_x$ is presented in Figure 12. One can see that the maximums of absolute values of errors was close to $1.64\times 10^{-4}$, i.e., the new schemes (2.16)–(2.21) are effective for the problems (1.1)–(1.3).

In addition, we provide the surfaces between the exact solution $E_y$ in Figure 13 and the numerical solution $E^n_y$ in Figure 14, as well as the contour of the error $E_y-E^n_y$ in Figure 15. One can see that the maximums of absolute values of errors was close to $6.92\times 10^{-5}$ from this Figure 15. From these figures we can see clearly that the numerical solution $E^n_y$ behavior can approximate well the exact solution $E_y$.

Similarly, Figures 16 and 17 show the surfaces between the exact solution $H_z$ and the numerical solution $H^n_z$. Figure 18 shows the contour of the error $H_z-H^n_z$, we can see that the maximums of absolute values of errors was close to $2.07\times 10^{-4}$. From these figures we can see clearly that the numerical solution $H^n_z$ behavior can approximate well the exact solution $H_z$, i.e., the the numerical schemes (2.16)–(2.21) are efficient.

Furthermore, taking $\Delta x = \Delta y = (\Delta t)^2$ and $T = 1$, some calculation results are shown in Tables 6 and 7. Table 6 shows the error of the numerical solution of (1.1)–(1.3) as computed by using (2.16)–(2.21) in the discrete $L^2$ norms and convergence rates in different time step sizes $\Delta t = T/N$. The results in Table 6, clearly show that the method is efficient and has fourth-order accuracy in time. $E_{1}$ which is the difference in energy between level $n$ and level $n+1$ and $E_{2}$ which is the difference between the two sides of the identity (3.1) are presented in Table 7. The results of $E_1$ and $E_2$ in Table 7 verify the theoretical results in Theorem 3.2.

Let $E_3$ be defined by

we show the errors $E_1$, $E_2$ and $E_3$ in Table 8 with $\Delta x = \Delta y = 0.01$ and $\Delta t = 0.001$. It is clearly shown that for a long time, the numerical solution of the schemes (2.16)–(2.21) keep the approximate energy-preserving nature in the discrete energy norms and are consistent with the theoretical results obtained in Theorem 3.2.

5.   Conclusions

In this paper, an ADI-FDTD method with fourth-order accuracy in time for the 2-D TE problem without sources and charges has been proposed and the numerical identity of the method has been derived. It is strictly proved that the proposed method is approximately energy-preserving, and that the two perturbation terms in the energy identity will affect the energy conservation in theory. Numerical experiments to compute the energies and convergence orders in time were carried out, and the computed results confirmed the theoretical analysis. The perturbation terms derived can be used to improve the application of the ADI-FDTD scheme.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant Nos. 11871138, 11801389, 11771314, 12171343), Sichuan Science and Technology Program (Grant Nos. 2020YJ0110, 2022JDTD0019), National-Local Joint Engineering Laboratory of System Credibility Automatic Verification (Grant No. ZD20220105) and funding from the V. C. and V. R. Key Lab of Sichuan Province.

Conflict of interest

The authors declare that they have no conflict of interest.