A new paralleling-in-order based unconditionally stable FDTD method using Legendre polynomials

Abstract

The orthogonal expansion based unconditionally stable FDTD method has been developed and attracted much attention recently. It can be implemented by two schemes: marching-on-in-order and paralleling-in-order solution, which respectively use weighted Laguerre polynomials (WLP) and associated Hermite (AH) functions as temporal expansions and testing functions. This paper proposed to use Legendre polynomials as another kind of orthogonal basis to make a new paralleling-in-order based unconditionally stable FDTD method. Numerical result shows it has much better performance in long time simulation case when compared with AH FDTD method.

Keywords

Finite difference time domain (FDTD)Legendre polynomials paralleling-in-order unconditionally stability

1. Introduction

To overcome the numerical stability constraints of conventional finite-difference time-domain (FDTD) method [1, 2, 3, 4], many unconditionally stable methods to reduce or eliminate requirements of the stability condition has been proposed and developed, such as alternating-direction implicit method [5] and locally one-dimensional schemes [6], explicit and unconditionally stable FDTD method [7], and orthogonal expansions in time domain [8, 9, 10, 12]. For the orthogonal expansions schemes, field-versus-time variations in the FDTD space lattice are expanded using an appropriate set of orthogonal temporal basis and testing functions, such as WLP and AH functions, which lead to two different solution schemes: marching-on-in-order and paralleling-in-order, respectively. Both of them appear to be promising according to the reported work where the computational time can be reduced to at least ten percent of the conventional FDTD scheme [1]. In this paper, we explore the Legendre (LD) polynomials as another possible orthogonal expansion to incorporate with FDTD to form a paralleling-in-order based unconditionally stable FDTD method. Numerical example validates the effectiveness and the advantages of this new method.

2. Formulation

The 2-D time-domain Maxwell’s equations with the TEz wave case in lossy medium are considered

$\displaystyle\varepsilon\frac{\partial E_{x}\left({r,t}\right)}{\partial t}+% \sigma_{e}E_{x}\left({r,t}\right)=\frac{\partial H_{z}\left({r,t}\right)}{% \partial y}-J_{x}\left({r,t}\right)$ (1)

$\displaystyle\varepsilon\frac{\partial E_{y}\left({r,t}\right)}{\partial t}+% \sigma_{e}E_{y}\left({r,t}\right)=-\frac{\partial H_{z}\left({r,t}\right)}{% \partial x}-J_{y}\left({r,t}\right)$ (2)

$\displaystyle\mu\frac{\partial H_{z}\left({r,t}\right)}{\partial t}+\sigma_{m}% H_{z}\left({r,t}\right)=\frac{\partial E_{x}\left({r,t}\right)}{\partial y}-% \frac{\partial E_{y}\left({r,t}\right)}{\partial x}-M_{z}\left({r,t}\right)$ (3)

where $\varepsilon$ , $\mu$ , $\sigma_{e}$ and $\sigma_{m}$ are the permittivity, the permeability, the electric conductivity, and the magnetic loss of the medium respectively. $E_{\xi}\left({r,t}\right)$ and $J_{\xi}\left({r,t}\right)$ ( $\xi=x,y$ ) are the electric field component and the electric current densities. $H_{z}\left({r,t}\right)$ and $M_{z}\left({r,t}\right)$ are the magnetic field component and magnetic current densities.

We expand all the temporal quantities in terms of the associated Legendre polynomial given by [11]

$\displaystyle P_{q}\left(t\right)=\sqrt{\frac{2q+1}{2}}L_{q}\left({2\frac{t}{l% }-1}\right),\mbox{ }t\in\left[{0,l}\right]$ (4)

where $l$ is the time support for analyzing a causal responses, and $L_{q}$ is the Legendre polynomial with order $q$ [11], which are orthogonal in the interval [–1, 1] satisfying the following recurrence relation

$\displaystyle L_{q+1}\left(t\right)=\frac{2q+1}{q+1}tL_{q}\left(t\right)-\frac% {q}{q+1}L_{q-1}\left(t\right),q\geqslant 1$ (5)

and $L_{0}\left(t\right)=0$ , $L_{1}\left(t\right)=t$ . Given an time-support field function $u\left({r,t}\right)$ , it can be expanded by Eq. (4) as

$\displaystyle u\left({r,t}\right)=\sum\limits_{q=0}^{\infty}{u_{q}\left(r% \right)}P_{q}\left(t\right)$ (6)

where $u_{q}\left(r\right)$ is the $q$ -th expanding coefficients, and it can be calculated by [11]

$\displaystyle u_{q}\left(r\right)=\int\limits_{-\infty}^{+\infty}{u\left({r,t}% \right)}P_{q}\left(t\right)dt$ (7)

From the intrinsic features of Legendre function [11], the differential relationship can be described as

$\displaystyle P_{q}\left(t\right)=\frac{1}{2\sqrt{\left({2q+3}\right)\left({2q% +1}\right)}}P_{q+1}^{\prime}\left(t\right)-\frac{1}{2\sqrt{\left({2q+1}\right)% \left({2q-1}\right)}}P_{q-1}^{\prime}\left(t\right)$ (8)

If the field derivative of $u\left({r,t}\right)$ to $t$ is expanded as

$\displaystyle u^{\prime}\left({r,t}\right)=\sum\limits_{q=0}^{\infty}{u_{q}^{% \left(1\right)}\left(r\right)P_{q}\left(t\right)}$ (9)

where $u_{q}^{\left(1\right)}\left(r\right)$ is $q$ -th expanding coefficients for $u^{\prime}\left({r,t}\right)$ , then incorporated with Eq. (8), it can be deduced as

$\displaystyle u^{\prime}(r,t)=\left(\sum\limits_{q=0}^{\infty}\left(\frac{l}{2% \sqrt{(2q+1)(2q-1)}}u_{q-1}^{(1)}(r)-\frac{l}{2\sqrt{(2q+3)(2q+1)}}u_{q+1}^{(1% )}(r)\right)P_{q}(t)\right)^{\prime}$ (10)

Connected Eqs (10) and (6), we can get

$\displaystyle u_{q}\left(r\right)=\frac{l}{2\sqrt{\left({2q+1}\right)\left({2q% -1}\right)}}u_{q-1}^{(1)}\left(r\right)-\frac{l}{2\sqrt{\left({2q+3}\right)% \left({2q+1}\right)}}u_{q+1}^{(1)}\left(r\right)$ (11)

When assemble $\left\{{u_{q}\left(r\right)}\right\}_{q=0,1\cdots Q-1}$ as a $Q$ -tuple $U$ , and $\{{u_{q}^{(1)}\left(r\right)}\}_{q=0,1\cdots Q-1}$ as $U^{\left(1\right)}$ , a matrix-multiply relationship can be obtained from Eq. (11) as following

$\displaystyle U=\alpha_{L}U^{\left(1\right)}$ (12)

where $\alpha_{L}$ is integral matrix

$\displaystyle\alpha_{L}=\frac{l}{2}\left[{\begin{array}[]{l}\ \ \ \ \ \ \ \ \ % {-1/}\sqrt{1\cdot 3}\\ \ \ \ \ {1/}\sqrt{1\cdot 3}\ \ \ \ \ {-1/}\sqrt{3\cdot 5}\\ {1/}\sqrt{3\cdot 5}\ \ \ \ \ \ \ \ \ \ \qquad\ddots\\ \ \ \ \ \ \ \ \ \ \ \ddots\qquad\qquad\ \ \ \ {-1/}\sqrt{\left({2Q-3}\right)% \left({2Q-1}\right)}\\ \qquad\ \ \ \ {1/}\sqrt{\left({2Q-3}\right)\left({2Q-1}\right)}\\ \end{array}}\right]_{Q\times Q}$ (13)

Alternatively Eq. (12) can be rewritten as

$\displaystyle U^{(1)}=\alpha_{L}^{-1}U$ (14)

Similar to the paralleling-in-order based AH FDTD method [12], we can apply a $Q$ -tuple-domain transformation to Eqs (1)–(3), and discretize them as following

$\displaystyle\alpha_{e\left({i,j}\right)}\left.{E_{x}}\right|_{i,j}=({\left.{H% _{z}}\right|_{i,j}-\left.{H_{z}}\right|_{i,j-1}})/\Delta\overline{y}_{j}-\left% .{J_{x}}\right|_{i,j}$ (15)

$\displaystyle\alpha_{e\left({i,j}\right)}\left.{E_{y}}\right|_{i,j}=-({\left.{% H_{z}}\right|_{i,j}-\left.{H_{z}}\right|_{i-1,j}})/\Delta\overline{x}_{i}-% \left.{J_{y}}\right|_{i,j}$ (16)

$\displaystyle\left.{\alpha_{m\left({i,j}\right)}H_{z}}\right|_{i,j}=({\left.{E% _{x}}\right|_{i,j+1}-\left.{E_{x}}\right|_{i,j}})/\Delta y_{j}-({\left.{E_{y}}% \right|_{i+1,j}-\left.{E_{y}}\right|_{i,j}})/\Delta x_{i}-\left.{M_{z}}\right|% _{i,j}$ (17)

where

$\displaystyle\alpha_{e\left({i,j}\right)}=\left.\varepsilon\right|_{i,j}\alpha% _{L}^{-1}+\left.{\sigma_{e}}\right|_{i,j}I$ (18)

$\displaystyle\alpha_{m\left({i,j}\right)}=\left.{\mu_{m}}\right|_{i,j}\alpha_{% L}^{-1}+\left.{\sigma_{m}}\right|_{i,j}I$ (19)

where $\left.{E_{x}}\right|_{i,j}$ , $\left.{E_{y}}\right|_{i,j}$ , $\left.{H_{z}}\right|_{i,j}$ , $\left.{J_{x}}\right|_{i,j}$ , $\left.{J_{y}}\right|_{i,j}$ and $\left.{M_{z}}\right|_{i,j}$ are $Q$ -tuple representations of fields and sources respectively. $I$ is the $Q$ -dimensional identity matrix. By assembling Eqs (15)–(17) and eliminating the electric field components, a five-diagonal banded matrix equation for Hz component can be obtained

$\displaystyle a_{l(i,j)}\left.{H_{z}}\right|_{i-1,j}+a_{r(i+1,j)}\left.{H_{z}}% \right|_{i+1,j}+a_{m(i,j)}\left.{H_{z}}\right|_{i,j}$ (20) $\displaystyle\ \ \ \ \ +a_{d(i,j)}\left.{H_{z}}\right|_{i,j-1}+a_{u(i,j+1)}% \left.{H_{z}}\right|_{i,j+1}=b_{i,j}$

where

$\displaystyle a_{u\left({i,j+1}\right)}=-\alpha_{e\left({i,j+1}\right)}^{-1}/% \Delta\overline{y}_{j+1}/\Delta y_{j}$ (21)

$\displaystyle a_{d\left({i,j}\right)}=-\alpha_{e\left({i,j}\right)}^{-1}/% \Delta\overline{y}_{j}/\Delta y_{j}$ (22)

$\displaystyle a_{l\left({i,j}\right)}=-\alpha_{e\left({i,j}\right)}^{-1}/% \Delta\overline{x}_{i}/\Delta x_{i}$ (23)

$\displaystyle a_{r\left({i+1,j}\right)}=-\alpha_{e\left({i+1,j}\right)}^{-1}/% \Delta\overline{x}_{i+1}/\Delta x_{i}$ (24)

$\displaystyle a_{m(i,j)}=-\left({a_{r(i+1,j)}+a_{l(i,j)}+a_{u(i,j+1)}+a_{d(i,j% )}+\alpha_{m\left({i,j}\right)}}\right)$ (25)

$\displaystyle b_{i,j}=-\left.\left({\left.{\alpha_{m(i,j+1)}^{-1}J_{x}}\right|% _{i,j+1}-\left.{\alpha_{m(i,j)}^{-1}J_{x}}\right|_{i,j}}\right)\right/\Delta y% +\left.\left({\left.{\alpha_{m(i+1,j)}^{-1}J_{y}}\right|_{i+1,j}-\left.{\alpha% _{m(i,j)}^{-1}J_{y}}\right|_{i,j}}\right)\right/\Delta x\left.{-M_{z}}\right|_% {i,j}$ (26)

By using eigenvalue transformation from $\alpha_{L}X=XV$ , where $X$ and $V$ are the eigenvector matrix and diagonal matrix composed of eigenvalues $\left\{{\lambda_{q}}\right\}$ respectively, Eq. (2) can be changed to the paralleling-in-order solution. For the $q$ -th decoupled equation, we have

$\displaystyle A\left({1/\lambda_{q}}\right)\left.{H_{z}^{\ast}}\right|^{q}=% \left.{b^{\ast}}\right|^{q}$ (27)

where $A\left(\cdot\right)$ is a banded sparse matrix, with the similar form as from AH FDTD method [12], and $\left.{b^{\ast}}\right|^{q}$ is the transformed variables from $\left.b\right|_{i,j}=X\left.{b^{\ast}}\right|_{i,j}$ . Finally, we can obtain a paralleling-in-order scheme to calculate all of the expanding coefficients of electromagnetic fields, and then the time-domain responses can be reconstructed from Eq. (6).

Figure 1.

The geometry configuration for a 2-D parallel plate waveguide with a PEC slot and a partly filled dielectric medium.

3. Numerical results

To verify the proposed new method, a numerical example is analyzed by simulating the TEz wave propagation in a parallel plate waveguide, as shown in Fig. 1. It is with a PEC slot of the thickness 0.2 mm and the distance 0.2 mm, and a partly filled dielectric material of the thickness 0.8 mm. The dielectric medium parameters are given as two cases: Case I, $\varepsilon=$ 11 $\varepsilon_{0}$ , $\mu=\mu_{0}$ , $\sigma_{e}=$ 0.003 S/m and $\sigma_{m}=$ 0 $\Omega$ /m; Case II, $\varepsilon=$ 2 $\varepsilon_{0}$ , $\mu=$ $\mu_{0}$ , $\sigma_{e}=$ 30000 S/m and $\sigma_{m}=$ 0 $\Omega$ /m. There are 140 $\times$ 8 uniform cells ( $\Delta_{x}=\Delta_{y}=$ 0.1 mm) in the computational domain. A Gaussian pulse sinusoidally modulated is used as the electric current source profile

$\displaystyle J_{y}\left(t\right)=\exp\left({-\left({\left({t-t_{c}}\right)/t_% {d}}\right)^{2}}\right)\sin\left({2\pi f_{c}\left({t-t_{c}}\right)}\right)$ (28)

where $t_{d}=1/(2f_{c})$ , $t_{c}=$ 4 $t_{d}$ , $f_{c}=$ 12 GHz. And the total simulation time is set as $l=$ 1.28 ns for case I and $l=$ 12.8 ns for case II, then it leads to the marching-in-on-time steps for $N=$ 6000 and $N=$ 60000, respectively, to meet the Courant-Friedrich-Levy (CFL) stability condition [1]. And the number of orders for LD functions is chosen as 80 and 300 respectively to obtain a good approximation of field components.

Table 1

The comparison of computational resources for the case of I

	${}_{\Delta t}$ (ps)	Memory (MB)	CPU time (s)
FDTD ( $N=$ 6000)	0.21	1.8	2.97
AH FDTD ( $Q=$ 80)	21	2.9	1.32
LD FDTD ( $Q=$ 80)	21	2.9	1.32

Figure 2.

The calculated results of transient electric field $E_{y}$ for the case of I. (a) the calculated waves (b) the relative error.

The $E_{y}$ electric field responses at measurement point $p_{1}$ and $p_{2}$ , located at the center of the slot and behind the medium respectively, are calculated, which are both in agreement with the conventional FDTD method as shown in Figs 2 and 3. For comparison, the AH FDTD method is also used in these two cases and both make the results from FDTD method as the reference value to calculate their relative error, which is defined by $R=20\log_{10}[{\left|{E_{y}^{\textit{AH/LD FDTD}}(t)-E_{y}^{\textit{FDTD}}(t)}% \right|/\max[{E_{y}^{\textit{FDTD}}(t)}]}]$ . For the case of I, one can find that the good results for these two methods in Fig. 2a, the relative errors are almost below 35 dB, and there is almost the same maximum value of relative error for these two methods, shown in Fig. 2b. But for the case of II, a difference come out as shown in Fig. 3a, where the result for LD FDTD method is still agreeable when number of orders is chosen as $Q=$ 300 while the result of AH FDTD method is not satisfying, unless it’s Q is chosen as 800. Especially, we can see this difference clearly from their relative errors, as shown in Fig. 3b and 3c, where the former are relevant to global time support and the latter shows its initial stage for detail representation. On the contrary, the accuracy seems to be insufficient for AH FDTD method when Q is selected as a relative small one, such as $Q=$ 300. Instead, it means that when Q reaches 800, the results from AH FDTD method can achieve a comparable accuracy with the ones from LD FDTD method. One should note that for case II the waveform at point $p_{2}$ has a larger amplitude attenuation and a longer delay than the result at point $p_{1}$ , due to the high dielectric medium located between them.

Table 2

The comparison of computational resources for the case of II

	${\Delta t}$ (ps)	Memory (MB)	CPU time (s)
FDTD ( $N=$ 60000)	0.21	1.8	30.8
AH FDTD ( $Q=$ 300)	21	11.8	1.55
AH FDTD ( $Q=$ 800)	21	28.9	1.95
LD FDTD ( $Q=$ 300)	21	11.8	1.55

Figure 3.

The calculated results of transient electric field $E_{y}$ for the case of II. (a) the calculated waves (b) the relative errors (c) the initial stage of relative errors.

Tables 1 and 2 show the comparison of the computational resources. We can see that the simulation takes much more time for the FDTD method compared with proposed method, especially for the case of II. While the tradeoff for the proposed method is that it consumes more memory than conventional FDTD method, which is similar to the AH FDTD method. One should note that the advantage of simulation time above for LD or AH FDTD methods tends to be obtained in the case of fine structures, such as the examples in this paper, for the conventional FDTD method should deal with a large number of marching-on-in-time steps under this case, which is a time-consuming process. In addition, from Table 2, we can find the advantages compared with AH FDTD method that the proposed method can use relative smaller memory storage and slightly fewer CPU times to get a readily results.

4. Conclusion

A new paralleling-in-order based unconditionally stable FDTD method is proposed using Legendre polynomials in this paper. By using the intrinsic differential features of the Legendre function, an integral matrix is derived to guarantee a direct $Q$ -tuple-domain transformation for time domain Maxwell equation. A numerical example with 2-D TEz case validates its effectiveness and shows it better performance, when compared with AH FDTD method, in long time simulation applications such as shielding effectiveness analysis. Although the proposed formulation in this paper is for 2D cases, it has a potential to be applied in analyzing the 3D structures in future. The main challenge might be how to implement a high efficient 3-D method based on this work in next step, for there will be thirteen nonzero elements in each line of the coefficient matrix of the banded equation form the 2-D case of Eq. (27) where there are only five nonzero elements, which will generate the increasing the computational resources in term of the CPU time and memory size.

Footnotes

Acknowledgments

This work is supported by the National Natural Science Foundation of China under Grant 51477183 and 51407198, and Natural Science Foundation of Jiangsu Province under Grant SBK2017043113. This support is gratefully acknowledged.

References

Taflove

and Hagness

S.C.

, Finite-Difference Time-Domain Solution of Maxwell’s Equations, Wiley Encyclopedia of Electrical and Electronics Engineering (2016), 1–33.

Xiong

Chen

Cai

Z.Y.

and Chen

, A numerically efficient method for the FDTD analysis of the shielding effectiveness of large shielding enclosures with thin-slots, Int J Appl Electrom 40(4) (Dec. 2012), 251–258.

Zhou

Y.H.

Huang

Z.Y.

Shi

L.H.

and Fu

S.C.

, Analysis of frequency-dependent field-to-transmission line coupling with Associated Hermite FDTD method, Int J Appl Electrom 49(4) (2015), 443–451.

Sha

Huang

Z.X.

Chen

M.S.

and Wu

X.L.

, Survey on symplectic finite-difference time-domain schemes for Maxwell’s equations, IEEE Trans. Antennas Propag. 56(2) (2008), 493–500.

Namiki

, A new FDTD algorithm based on alternating-direction implicit method, IEEE Trans. Microwave Theory Tech. 47 (Oct. 1999), 2003–2007.

Ahmed

Chua

E.K.

E.P.

and Chen

, Development of the three-dimensional unconditionally stable LOD-FDTD method, IEEE Trans. Antennas Propag. 56(11) (Nov. 2008), 3596–3600.

Gaffar

and Jiao

, An explicit and unconditionally stable FDTD method for electromagnetic analysis, IEEE Trans. Microw. Theory Tech. 62(11) (Nov. 2014), 2538–2550.

Chung

Y.S.

Sarkar

T.K.

Jung

B.H.

and Salazar-Palma

, An unconditionally stable scheme for the finite-difference time-domain method, IEEE Trans. Microw. Theory Tech. 51(3) (2003), 697–704.

Zhao

J.Y.

Yin

W.Y.

Zhu

M.D.

and Luo

, Time domain EFIEPMCHW method combined with adaptive marching-on-in order procedure for studying on time-and frequency-domain responses of some composite structures, IEEE Trans. Electromagn. Compat. 55(6) (Dec. 2013), 1220–1230.

10.

Huang

Z.Y

Shi

L.H.

Chen

and Zhou

Y.H.

, A new unconditionally stable scheme for FDTD Method Using Associated Hermite orthogonal functions, IEEE Trans. Antennas Propagat. 62(9) (Sept. 2014), 4804–4809.

11.

Bayin

Ş.S.

, Legendre Equation and Polynomials, in Mathematical Methods in Science and Engineering, John Wiley & Sons, Inc., Hoboken, NJ, USA. 2006.

12.

Huang

Z.Y.

Shi

L.H.

Zhou

Y.H.

and Chen

, An Improved Paralleling-in-order Solving Scheme for AH-FDTD Method using Eigenvalue Transformation, IEEE Trans. Antennas Propagat. 63(5) (2015), 2135–2140.