1. Introduction
The finite-difference time-domain (FDTD) method [1], [2] is a widely used numerical method for the analysis of electromagnetic problems. Although the explicit FDTD method has the advantage of a simple algorithm, it has the drawback that the Courant-Friedrichs-Lewy (CFL) condition restricts the choice of the time step size. To remove the CFL condition, we often resort to the implicit FDTD methods such as the alternating-direction implicit (ADI) [3], [4], locally one-dimensional (LOD) [5]-[7], and Crank-Nicolson (CN) [8], [9] FDTD methods.
The ADI- and LOD-FDTD methods can be efficiently implemented due to the operator-splitting technique. Note that the use of these methods yields the operator splitting error, particularly when the time step size is chosen to be large. On the other hand, the original CN-FDTD method without operator splitting is highly accurate, but requires solving a large banded matrix with a time-consuming solver. To alleviate this problem, we have introduced the iterated CN (ICN) scheme into the FDTD method [10], in which an implicit CN procedure can be replaced with an explicit iteration process [11]. Note that the ICN-FDTD method is constrained by the CFL condition as the explicit FDTD method [10]. Recently, the ICN-FDTD method has been extended with the domain decomposition [12], body-of-revolution technique [13] and frequency-dependent formulation [14]. Although the validity of the ICN-FDTD method has been demonstrated, its numerical dispersion property has not yet been discussed.
In this letter, the numerical dispersion analysis of the ICN-FDTD method is performed. We first review the formulation of the ICN-FDTD method with two iterations. We next derive the numerical dispersion relation of the ICN-FDTD method from the amplification matrix. The numerical dispersion property is discussed with attention to the eigenvalue of the amplification matrix.
2. Formulation
We consider linear, isotropic, and lossless materials in a one-dimensional problem. Maxwell’s equations in a matrix form are expressed as
| \[\begin{align} \frac{{\partial}\boldsymbol{\phi}}{{\partial}t}=\boldsymbol{A\phi} \tag{1} \end{align}\] |
where \(\boldsymbol{\phi}=[H_x, E_y]^T\) and
| \[\begin{aligned} \boldsymbol{A}= \begin{bmatrix} 0&\frac{1}{\mu_0}{\frac{\partial}{\partial z}}\\ \frac{1}{\varepsilon_0}{\frac{\partial}{\partial z}}&0\\ \end{bmatrix} \end{aligned}\] |
in which \(T\) is the transpose, \(\varepsilon_0\) and \(\mu_0\) represent the permittivity and permeability of a vacuum, respectively.
We apply the two-iteration technique to Eq. (1). First, Eq. (1) is calculated as
| \[\begin{align} \boldsymbol{\phi}^{n+1}_{\rm 1st}=\boldsymbol{\phi}^n+{{\Delta}t}\boldsymbol{A\phi}^n \tag{2} \end{align}\] |
where \(\Delta \it{t}\) is the time step. Using the result of Eq. (2), the intermediate field is estimated as
| \[\begin{align} \boldsymbol{\phi}^{n+\frac{1}{2}}_{\rm 1st}=\frac{1}{2}(\boldsymbol{\phi}^{n+1}_{\rm 1st}+\boldsymbol{\phi}^n). \tag{3} \end{align}\] |
Using the field of Eq. (3), the second iteration is carried out as
| \[\begin{align} \boldsymbol{\phi}^{n+1}_{\rm 2nd}=\boldsymbol{\phi}^n+{{\Delta}t}\boldsymbol{A\phi}^{n+\frac{1}{2}}_{\rm 1st}. \tag{4} \end{align}\] |
The improved intermediate field is
| \[\begin{align} \boldsymbol{\phi}^{n+\frac{1}{2}}_{\rm 2nd}=\frac{1}{2}(\boldsymbol{\phi}^{n+1}_{\rm 2nd}+\boldsymbol{\phi}^n). \tag{5} \end{align}\] |
Finally, using the result of Eq. (5), the solution is obtained as
| \[\begin{align} \boldsymbol{\phi}^{n+1}=\boldsymbol{\phi}^n+{\Delta}t {{\boldsymbol{A\phi}^{n+\frac{1}{2}}_{\rm 2nd}}}. \tag{6} \end{align}\] |
The FDTD equations are obtained with substituting the above \(\boldsymbol{\phi}\) and \(\boldsymbol{A}\) into Eqs. (2)-(6).
3. Numerical Dispersion Analysis
It is important to study the wavelength characteristics of the numerical dispersion, since they offer useful information on the accuracy of the method by which the wideband characteristics are treated [15]. We here derive the numerical dispersion relation, paying attention to the eigenvalue of the amplification matrix.
A direct solution without the intermediate field is obtained from Eqs. (2)-(6), leading to [10]
| \[\begin{align} \boldsymbol{\phi}^{n+1}=\Big(\boldsymbol{I}+\Delta t\boldsymbol{A}+\frac{\Delta t^2}{2}\boldsymbol{A}^2+\frac{\Delta t^3}{4}\boldsymbol{A}^3\Big)\boldsymbol{\phi}^{n} \tag{7} \end{align}\] |
where \(\boldsymbol{I}\) is the identity matrix. Here, we define the plane wave as
| \[\begin{align} & H_x=H_{x0}e^{j(\omega n \Delta t-k_zl\Delta z)} \tag{8} \\ & E_y=E_{y0}e^{j(\omega n \Delta t-k_zl\Delta z)} \tag{9} \end{align}\] |
in which \(n\) and \(l\) denote the indexes in the \(t\) and \(z\) axes, respectively, \(\omega\) is the angular frequency, \(k_z=2\pi/\lambda\) is the wavenumber, \(\lambda\) is the wavelength, and \(\Delta z\) is the spatial sampling width. Substituting Eqs. (8) and (9) into Eq. (7), we obtain the following matrix form:
| \[\begin{align} \boldsymbol{\phi}^{n+1}=\boldsymbol{F\phi}^{n} \tag{10} \end{align}\] |
where
| \[\begin{aligned} & \boldsymbol{F}= \begin{bmatrix} 1-2\alpha^2&-j\frac{2}{c\mu_0}\alpha(1-\alpha^2)\\ -j\frac{2}{c\varepsilon_0}\alpha(1-\alpha^2)&1-2\alpha^2\\ \end{bmatrix} \\ & \alpha=\frac{c\Delta t}{\Delta z}\sin{\frac{k_z \Delta z}{2}} \end{aligned}\] |
in which \(\boldsymbol{F}\) is the total amplification matrix, and \(c\) is the speed of light in a vacuum. Equation (10) can be rewritten as
| \[\begin{align} (e^{j\omega \Delta t}\boldsymbol{I}-\boldsymbol{F})\boldsymbol{\phi}^n=0. \tag{11} \end{align}\] |
For Eq. (11) to have a solution, the following relation must be satisfied [16]:
| \[\begin{align} \det(e^{j\omega \Delta t}\boldsymbol{I}-\boldsymbol{F})=0. \tag{12} \end{align}\] |
After some manipulations, we finally derive the following numerical dispersion relation:
| \[\begin{align} e^{j\omega \Delta t}=1-2\alpha^2\pm j2\alpha(1-\alpha^2). \tag{13} \end{align}\] |
We next derive the eigenvalues of \(\boldsymbol{F}\). The absolute value of the eigenvalue is obtained as
| \[\begin{align} |\lambda_{\rm eig}|=\sqrt{1-4\alpha^4+4\alpha^6}. \tag{14} \end{align}\] |
For Eq. (10) to be stable, \(|\lambda_{\rm eig}| \leq 1\) must be satisfied. Note that Eq. (14) is less than one, since \(|\alpha|\) is less than 1 as long as the CFL condition is satisfied. That is, the time step of the ICN-FDTD method is restricted by the same CFL condition as in the traditional explicit one [10]. For this case, the field is dissipative. To account for the dissipative property, a complex angular frequency \((\omega=\omega_{\rm R}+j\omega_{\rm I})\) should be introduced into Eq. (13) as
| \[\begin{align} e^{-\omega_{\rm I}\Delta t}(\cos{\omega_{\rm R} \Delta t} {+} j\sin{\omega_{\rm R}\Delta t}) = 1 {-} 2\alpha^2\pm j2\alpha(1 {-} \alpha^2). \tag{15} \end{align}\] |
Comparing the real and imaginary terms in the above equation, we can rewrite the numerical dispersion relation as follows:
| \[\begin{align} \begin{split} \omega=\frac{1}{\Delta t}\arcsin{\frac{2\alpha-2\alpha^3}{\sqrt{1-4\alpha^4+4\alpha^6}}}\\-j\frac{1}{\Delta t}\log{\sqrt{1-4\alpha^4+4\alpha^6}}. \end{split} \tag{16} \end{align}\] |
It is found that the numerical dispersion relation of the ICN-FDTD method becomes complex. The imaginary part represents the dissipative property.
Figure 1 (a) and (b) show the normalized phase velocity and normalized attenuation factor, respectively, as a function of sampling width normalized to the wavelength (\(\Delta z/ \lambda\)). Here, a width of \(\Delta z/ \lambda=0.02\) and \(0.05\) corresponds to the case where one wavelength is sampled by fifty and twenty points, respectively. The numerical dispersion property is calculated from \(\omega/(k_z c)\) using \(\omega\) obtained from Eq. (16), that is its real and imaginary parts correspond to the normalized phase velocity and normalized attenuation factor, respectively. In this calculation, \(\Delta t\) is the maximum time step that satisfies the CFL condition of the explicit FDTD method. For reference, the results of the explicit FDTD and CN-FDTD methods are also included. It is seen that the ICN-FDTD method has a small phase error relative to the explicit method, but exhibits almost the same one as the CN counterpart. Unfortunately, only the ICN-FDTD method is slightly dissipative resulting from the attenuation factor, since the iterative calculation terminates the series expansion of the CN procedure. Nevertheless, the ICN-FDTD method yields almost the same numerical results as those from the explicit FDTD method when analyzing optical devices [10], [14].
 |
Fig. 1 Numerical dispersion properties. |
4. Conclusion
We have investigated the numerical dispersion property of the FDTD method based on the ICN scheme. The numerical dispersion relation is derived with attention to the eigenvalue of the amplification matrix. For the eigenvalue less than one, it is shown that the ICN-FDTD method has a small phase error and is slightly dissipative. Our future work is to show the effectiveness of the ICN-FDTD method compared to the explicit FDTD method.
References
[1] K. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propag., vol.14, no.3, pp.302-307, May 1966.
CrossRef
[2] A. Taflove and S.C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method, 3rd ed., Artech House, Norwood, MA, 2005.
[3] T. Namiki, “A new FDTD algorithm based on alternating-direction implicit method,” IEEE Trans. Microw. Theory Tech., vol.47, no.10, pp.2003-2007, Oct. 1999.
CrossRef
[4] F. Zheng, Z. Chen, and J. Zhang, “A finite-difference time-domain method without the Courant stability conditions,” IEEE Microw. Guided Wave Lett., vol.9, no.11, pp.441-443, Nov. 1999.
CrossRef
[5] J. Shibayama, M. Muraki, J. Yamauchi, and H. Nakano, “Efficient implicit FDTD algorithm based on locally one-dimensional scheme,” Electron. Lett., vol.41, no.19, pp.1046-1047, Sept. 2005.
CrossRef
[6] E.L. Tan, “Unconditionally stable LOD-FDTD method for 3-D Maxwell’s equations,” IEEE Microw. Wireless Compon. Lett., vol.17, no.2, pp.85-87, Feb. 2007.
CrossRef
[7] I. Ahmed, E.-K. Chua, E.-P. Li, and Z. Chen, “Development of the three-dimensional unconditionally stable LOD-FDTD method,” IEEE Trans. Antennas Propag., vol.56, no.11, pp.3596-3600, Nov. 2008.
CrossRef
[8] Y. Yang, R.S. Chen, and E.K.N. Yung, “The unconditionally stable Crank Nicolson FDTD method for three-dimensional Maxwell’s equations,” Microw. Opt. Tech. Lett., vol.48, no.8, pp.1619-1622, May 2006.
CrossRef
[9] H.K. Rouf, F. Costen, and S.G. Garcia, “3D Crank-Nicolson finite difference time domain method for dispersive media,” Electron. Lett., vol.45, no.19, pp.961-962, Sept. 2009.
CrossRef
[10] J. Shibayama, T. Nishio, J. Yamauchi, and H. Nakano, “Explicit FDTD method based on iterated Crank-Nicolson scheme,” Electron. Lett., vol.58, no.1, pp.16-18, Jan. 2022.
CrossRef
[11] S.A. Teukolsky, “Stability of the iterated Crank-Nicholson method in numerical relativity,” Phys. Rev. D, vol.61, no.8, pp.087501-1-2, March 2000.
CrossRef
[12] P. Wu, X. Wang, Y. Xie, H. Jiang, and T. Natsuki, “Iterated Crank-Nicolson procedure with enhanced absorption for nonuniform domains,” IEEE J. Multiscale Multiphys. Comput. Tech., vol.7, pp.61-68, March 2022.
CrossRef
[13] S. Wu, Y. Dong, L. Liu, F. Su, and X. Chen, “ICN-FDTD scheme with absorption boundary condition for nonuniform rotational symmetric geometrics,” Int. J. Numer. Model. Electron. Networks Devices Fields, vol.36, no.2, e3057, March/April 2023.
CrossRef
[14] J. Shibayama, A. Kawahara, J. Yamauchi, and H. Nakano, “Frequency-dependent finite-difference time-domain method based on iterated Crank-Nicolson scheme,” Electron. Lett., vol.59, no.1, e12695. Jan. 2023.
CrossRef
[15] J. Shibayama, M. Muraki, R. Takahashi, J. Yamauchi, and H. Nakano, “Performance evaluation of several implicit FDTD methods for optical waveguide analyses,” J. Lightw. Technol., vol.24, no.6, pp.2465-2472, June 2006.
CrossRef
[16] S. Ogurtsov and G. Pan, “An updated review of general dispersion relation for conditionally and unconditionally stable FDTD algorithms,” IEEE Trans. Antennas Propag., vol.56, no.8, pp.2572-2583, Aug. 2008.
CrossRef
