Investigation of Numerical Dispersion with Time Step of The FDTD Methods: Avoiding Erroneous Conclusions

Without loss of generality, a homogenous, lossless, isotropic medium and uniform mesh is considered in our study. Therefore, the analytical numerical dispersion of the FDTD(2,2) method [30] can be expressed as

where ${\tilde{k}_{x}}=\tilde{k}\sin\theta\cos\varphi$, $\tilde{k}_{y}=\tilde{k}\sin\theta\sin\varphi$, ${\tilde{k}_{z}}=\tilde{k}\cos\theta$ are numerical wavenumbers in the $x$, $y$ and $z$ direction, respectively, $\theta$ and $\phi$ are the azimuth and zenith angles, $\tilde{k}$ is the numerical wavenumber of electromagnetic waves in the Yee’s grid, $\Delta x$, $\Delta y$ and $\Delta z$ are mesh sizes in the $x$, $y$ and $z$ direction, respectively. $\omega$ denotes the angular frequency and $\Delta t$ is the time step.

To obtain stable numerical solutions for the explicit FDTD (2,2) method, time steps must satisfy the CFL condition [30], expressed as

For convenience in our following derivation, $S_{fdtd(2,2)}$ is defined as

where $\Delta t_{max\_{fdtd(2,2)}}$ is the maximum time step defined by the CFL condition in (3).

The analytical numerical dispersion of the FDTD(2,4) method [30] is written as

Then, the CFL condition for the FDTD(2,4) method [30] is expressed as

We further define $S_{fdtd(2,4)}$ as

where $\Delta t_{max\_{fdtd(2,4)}}$ is the maximum time step constrained by the CFL condition in (6).

Remark 1: $\tilde{k}\geq k$, $\Delta t\in(0,t_{max\_fdtd(2,2)}]$.

Since the explicit expression of $\tilde{k}$ can not be directly obtained from (2), the parameter scanning method [36, 37] is used to study the relationship between $\tilde{k}$ and $k$. We scanned $\theta\in[0,\pi]$, $\phi\in[0,2\pi]$ with other possible parameters in (2). This above statement always holds true. Although it is not mathematically rigorous, it should be enough for practical applications and our analysis. Here, we plot the maximum $\tilde{k}$ with $f=5$ GHz, $\Delta x$ = $\Delta y$ = $\Delta z$ = $6\times{10^{-3}}$ m. Both $\theta\in[0,\pi]$, $\phi\in[0,2\pi]$ are scanned and the maximum $\tilde{k}$ is obtained as shown in Fig. 1. It can be easy to find that when $\Delta t$ satisfies the CFL condition, namely $S\in(0,1]$, $\tilde{k}$ is always larger than $k$.

Lemma 1: For the FDTD(2,2) method,

Proof: Since the explicit expression of $\tilde{k}$ is not available, the implicit differentiation [38] is resorted to investigating the derivative of $\tilde{k}$ with respect to $S$ defined in (4). For convenience of derivations, the following symbols are defined as

and

Therefore, (2) is separated into LHS and RHS two parts. We first consider LHS. By assuming uniform mesh used in the simulations, $\Delta x=\Delta y=\Delta z=\Delta$, and substituting $\omega={{2\pi}}/{{\lambda\sqrt{\varepsilon\mu}}}$ and (4) into LHS, we obtain

where $N$ is the count of sampling cells per wavelength, defined as $N=\lambda/\Delta$. For practical simulations, $N\geq 10$ should be satisfied to obtain acceptable accurate results [30].

With some mathematical manipulations, the derivative of LHS with respect to $S$ is obtained as

where $Q=\left[{\pi S\cos\left({{{\pi S}}/{{\sqrt{3}N}}}\right)-\sqrt{3}N\sin\left({{{\pi S}}/{{\sqrt{3}N}}}\right)}\right]$. Since ${{\pi S}}/({{\sqrt{3}N}})\in(0,{\pi}/{2})$,

Then, let’s check the sign of $Q$ in (12). By taking its derivative, we get

Since $Q$ is a continous function, ${-{\pi^{2}}S\sin\left({\pi S/\sqrt{3}N}\right)/\left({\sqrt{3}N}\right)<0}$, and ${\left.{\left[{\pi\cos\left({{{\pi S}}/{{\sqrt{3}N}}}\right)S-\sqrt{3}N\sin\left({{{\pi S}}/{{\sqrt{3}N}}}\right)}\right]}\right|_{S=0}}=0$, we get

By considering (13), (15) and other variable signs in (12), we can easily get

With a similar manner, the derivative of RHS can be expressed as

where $\tilde{k}=\omega/{\tilde{v}_{p}}$, where $\tilde{v}_{p}$ is the numerical phase velocity in the Yee’s grid, and

with $M={{{\pi a\sin\theta\cos\varphi}}/{N}}$, $T={{{\pi a\cos\theta}}/{N}}$, and $a={c_{0}}/\tilde{v}_{p}$.

To further investigate signs of (18), we divide its domain into several parts based on signs of their function values of each term. Note that the function domain is $\varphi\in[0,2\pi]$ and $\theta\in[0,\pi]$. Therefore, we divide them into several parts to check their signs. The detailed results are summarized in Table I. Table I (a)-(c) correspond to three terms in (18). To better illustrate it, we take term1 in (18) as an example. When $\theta\in\left[{0,\pi/2}\right]$ and $\varphi\in\left[{0,\pi/2}\right]$, it is evidently that $\sin\theta\in\left[{0,1}\right],\cos\varphi\in\left[{0,1}\right]$. Note that $(\pi a/N)\in\left({0,\pi/2}\right)$, then $M\in[{0,\pi/2}]$. Therefore, we have $\sin\left(M\right)\geq 0$ and $\cos\left(M\right)\geq 0$, which implies that term1 is positive when $\theta\in\left[{0,\pi/2}\right]$ and $\varphi\in\left[{0,\pi/2}\right]$. The remaining situations can be obtained with similar manners. To sum them up, we can obtain the following inequality

By considering (16) and (19), we get

Evidently, Lemma 1 is analytically proved and we can further analyze effect of time steps on the NDE of the FDTD(2,2) method.

From Remark 1 and Lemma 1, we can easily find that the numerical wavenumber $\tilde{k}$ monotonically decreases as $S$ in the FDTD(2,2) method increases. When time steps grow larger and are bounded by (3), the numerical phase velocity $\tilde{v}_{p}$ of electromagnetic waves is closer to its physical counterpart $c_{0}$, and then NDE caused by the numerical dispersion would be relatively smaller. Therefore, the NDE reaches its minimum value when $S$ gets its maximum value defined by (3).

With similar procedures, we can investigate effects of time steps upon numerical dispersion for the FDTD(2,4) method.

Remark 2: $\exists t\_\in(0,t_{max\_fdtd(2,4)}]$, we have $\tilde{k}\geq k$, $\Delta t\in(0,t\_],\tilde{k}\leq k$, $\Delta t\in(t\_,t_{max\_fdtd(2,4)}]$.

Lemma 2: For the FDTD(2,4) method,

From Remark 2 and Lemma 2, the NDE $\tilde{k}$ firstly decreases as the $S$ of the FDTD(2,4) method get larger. When $\Delta t\in(t\_,t_{max\_fdtd(2,4)}]$, the NDE becomes larger. Therefore, the NDE reaches its minimum value at $t=t\_$.

To comprehensively investigate the numerical dispersion of the FDTD (2,2) method, the numerical phase velocity $\tilde{v}_{p}$ is calculated from (2). In our numerical studies, the frequency is 5 GHz and $\Delta x$ = $\Delta y$ = $\Delta z$ = $6\times{10^{-3}}$ m, which corresponds to $\lambda/10$. In all other simulations, the same parameters are used, otherwise stated.

As shown in Fig. 2, $\tilde{v}_{p}$ changes periodically over $\varphi$ for a fixed time step $S_{fdtd(2,2)}$ when $\theta=90^{o}$. However, $\tilde{v}_{p}$ increases as $S_{fdtd(2,2)}$ becomes larger with fixed $\phi$ values, which implies that $\tilde{v}_{p}$ approaches $c_{0}$ as time steps grow larger. Therefore, the NDE would become smaller as time steps get larger, which agrees with our analytical analysis in the previous section. Fig. 3 shows $\tilde{v}_{p}$ with respect to $S_{fdtd(2,2)}$ when $\phi=90^{o}$. Similar statements can be also obtained.

Fig. 4 illustrates $\tilde{v}_{p}$ with different $S_{fdtd(2,2)}$ with respect to $\varphi$ and $\theta$. It is easy to find that each $\tilde{v}_{p}$ surface does not intersect each other. Therefore, we can conclude that $\tilde{v}_{p}$ approaches $c_{0}$ as $S_{fdtd(2,2)}$ becomes larger. The numerical dispersion results show that the optimum time step of the FDTD(2,2) method is the largest time step defined by the CFL condition, which also agrees with our theoretical analysis.

$\tilde{v}_{p}/c_{0}$ of the FDTD (2,4) method is also calculated to investigate its numerical dispersion with different $S_{fdtd(2,4)}$. As shown in Fig. 5, $\tilde{v}_{p}/c_{0}$ changes periodically with respect to $\varphi$ when $\theta=90^{\circ}$. $\tilde{v}_{p}$ approaches to $c_{0}$ at the beginning when $S_{fdtd(2,4)}$ becomes large, and then becomes larger than $c_{0}$ as $S_{fdtd(2,4)}$ further increases. Therefore, the NDE of the FDTD (2,4) method firstly decreases and then increases over $S_{fdtd(2,4)}$. The statements obtained in Fig. 6 are similar to those of Fig. 5. They agree with our theoretical analysis.

Here, we provide a numerical method to calculate the optimum time step of the FDTD(2,4) method. The optimum time step is defined as

which can make difference between the numerical wavenumber and its analytical counterpart with $\theta\in[0,\pi]$ and $\varphi\in[0,2\pi]$ minimum.

A wave propagation example is used to verify the accuracy of the two FDTD methods with different time steps in one dimensional case. A Gaussian function with the waveform $\phi\left(t\right)={e^{-16{{\left({0.7-ct}\right)}^{2}}}}$ is selected as the excitation, and uniform mesh size $\Delta=5\times{10^{-2}}$ m, and the simulation time is $t=3.6685\times{10^{-8}}$ s. Note that reflected electromagnetic waves do not reach probe in the simulations.

A two dimensional cavity with dimension of $1$ m $\times$ $2$ m and perfectly electric conductor (PEC) boundaries is considered and its resonant frequencies in both transverse magnetic (TM) and transverse electric (TE) modes are calculated through the two FDTD methods. The mesh size is $\Delta=4\times{10^{-2}}$ m.

Case A: a cubic cavity with side length of $1$ m in the three dimensional case is considered. The mesh size is $\Delta=4\times{10^{-2}}$ m. The resonant frequencies are calculated from the two FDTD methods.

Fig. 14 shows the relative error of resonant frequencies obtained from the FDTD(2,2) method. It is similar to one- and two-dimensional cases. Although there are some numerical fluctuations, the relative error becomes smaller for all the three modes as $S_{fdtd(2,2)}$ increases. It agrees well with our investigations.

Fig. 15 shows the resonant frequencies obtained from the FDTD(2,4) method. It is easy to find that relative errors of resonant frequencies show similar behaviors. The relative error gradually decreases and then slightly becomes larger as $S_{fdtd(2,4)}$ become larger. The relative error with small $S_{fdtd(2,4)}$ is larger than that with large $S_{fdtd(2,4)}$. It may account for other numerical errors.

Case B: we use the FDTD(2,2) method to further verify our findings for a practical application in this case. A missile is illustrated by a plane wave and the surface current density is calculated from the FDTD(2,2) method with different time steps. The geometry configurations of the missile in our simulations are shown in Fig. 16(a). It is with dimension of 2.3114$\lambda$$\times$1.4647$\lambda$$\times$0.5334$\lambda$. The plane wave incidents from the $x$-axis with 700 MHz. Fig. 16(b) shows the reference surface current density calculated from the FEKO. Fig. 17(a) and (b) show the surface current density obtained in the FDTD(2,2) method with $S_{fdtd(2,2)}=0.1$ and $1$. It can be found that the surface current distribution obtained from the FDTD(2,2) method agree well with the reference solution. The accuracy of results in our simulations seem unchanged with $S_{fdtd(2,2)}=0.1$ and $1$. There are various factors that can account for those results in this case, like complex geometry, staircase error, possible reflection from the total field/scattered field boundary. It turns out that time steps of the FDTD(2,2) method almost have no significant effects on the accuracy in the practical simulation. Therefore, large time step is preferred in the practical simulations since short simulation time is required and almost the same level of accuracy can also be achieved.