# High-order unconditionally-stable four-step ADI-FDTD methods and numerical analysis.

1. INTRODUCTIONThe finite-difference time-domain (FDTD) method [1] has been proven to be an established numerical technique that provides accurate predictions of field behaviors for electromagnetic interaction problems [2-4]. Some enhanced FDTD methods have been proposed [5-9]. However, the applications of the FDTD method had been restricted by the well-known Courant-Friedrichs-Lewy (CFL) condition [10] on the times step and the numerical dispersion associated with space discretization.

Recently, to overcome the CFL condition on the time step size of the FDTD method, an unconditionally-stable alternating direction implicit (ADI)-FDTD method was developed [11,12]. The ADI-FDTD method has second-order accuracy both in time and space. Moreover, the numerical dispersion of two and three dimensions was analyzed in [13, 14]. Nevertheless, it presents large numerical dispersion error with large time steps. To improve the dispersion performance, several methods were proposed, such as parameter-optimized [15], error-reduced [16], artificial-anisotropy methods [17,18] and the simplified sampling biorthogonal ADI method [19].

Along the same line, other unconditionally-stable methods such as Crank-Nicolson (CN) [20, 21], split-step (SS) [22-25] and locally-one dimensional (LOD) [26, 27] FDTD methods were developed. Specially, the CN-FDTD method has higher accuracy than that of the ADI-FDTD method. The LOD-FDTD method can be considered as the split-step approach (SS1) with first-order accuracy in time, which consumes less CPU time than that of the ADI-FDTD method. On the other hand, using high-order schemes is a usual method to reduce the numerical dispersion error. Concretely, a higher-order ADI-FDTD method in 2-D domains was presented in [28]. Then, a comprehensive study of the stability and dispersion characteristics for a set of higher order 3-D ADI-FDTD methods was presented in [29]. Subsequently, an arbitrary-order 3-D LOD-FDTD method was proposed in [30]. Furthermore, high-order split-step FDTD methods in 3-D domains were developed in [31]. An arbitrary-order leapfrog ADI-FDTD method was presented in [32].

High-order four-step ADI-FDTD methods in 3-D domains are presented in this paper. Firstly, based on the exponential evolution operator (EEO), the Maxwell's equations can be split into four sub-procedures. Accordingly, the time step is divided into four sub-steps, and the proposed methods are denoted by four-step ADI-FDTD herein. In addition, high-order central finite-difference operators based on the Taylor central finite-difference method in [33] are used to approximate the spatial differential operators. Moreover, the uniform formulation of the proposed high-order methods is generalized. Secondly, all the high-order schemes are proven to be unconditionally stable by using the Fourier method. Furthermore, the numerical dispersion relation is derived analytically. Finally, numerical experiments are presented to demonstrate the unconditional stability and the efficiency of the proposed methods. Moreover, some discussions are provided to show the effects of the order of schemes, the propagation angle, the time step, and the mesh size on the numerical dispersion. It can be concluded that the proposed methods achieve better accuracy even with the coarser mesh, and such an improvement actually leads to other advantages such as higher computational efficiency and lower memory requirements.

2. FORMULATIONS OF THE PROPOSED HIGH-ORDER METHODS

In linear, lossless and isotropic medium, [epsilon] and [mu] are the electric permittivity and magnetic permeability, respectively. Then, the 3-D Maxwell's equations can be written in a matrix form as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (1)

The matrix [M] is split into four parts [A]/2, [B]/2, [A]/2, and [B]/2. Then, (1) can be written as

[partial derivative][??]/[partial derivative]t = [A]/2 x [??] + [B]/2 x [??] + [A]/2 x [??] + [B]/2 x [??]. (2)

Suppose that a numerical solution u(t) at a given time [t.sup.n] = n[DELTA]t is transported to that at the next time [t.sup.n+1] = (n + 1)[DELTA]t. Now, from the forward Taylor series development

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (3)

Therefore, combination with (3), the solution to (2) can be easily found as

[[??].sup.n+1] = exp([DELTA]t/2 x [A] + [DELTA]t/2 x [B] + [DELTA]t/2 x [A] + [DELTA]t/2 x [B])[[??].sup.n]. (4)

The exponential evolution operator (EEO) in (4) can be reformulated as follows

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (5)

Note that in the above equation [A] + [B] = [B] + [A] but [A][B] [not equal to] [B][A].

If the sequential splitting is used to split these EEOs, we can obtain

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (6)

By using the following Taylor series approximation,

exp ([delta][A]) = [[infinity].summation over (k=0)][([delta][A]).sup.k]/k! [approximately equal to] 1 + [delta][A] for [absolute value of [delta]] [much less than] 1. (7)

(6) can be approximated in the following manner for small [DELTA]t

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (8)

Now intermediate variables [[??].sup.n+1/4], [[??].sup.n+2/4], and [[??].sup.n+3/4] are introduced in between [[??].sup.n] and [[??].sup.n+1]. Then (8) can be computed in the following four sub-steps

([I] - [DELTA]t/4 x [A]) [[??].sup.n+1/4] = ([I] + [DELTA]t/4 x [B]) [[??].sup.n] (9a)

([I] - [DELTA]t/4 x [B]) [[??].sup.n+2/4] = ([I] + [DELTA]t/4 x [A]) [[??].sup.n+1/4] (9b)

([I] - [DELTA]t/4 x [A]) [[??].sup.n+3/4] = ([I] + [DELTA]t/4 x [B]) [[??].sup.n+2/4] (9c)

([I] - [DELTA]t/4 x [B]) [[??].sup.n+1] = ([I] + [DELTA]t/4 x [A]) [[??].sup.n+3/4]. (9d)

where [I] is a 6 x 6 identity matrix. [partial derivative]/[partial derivative][alpha] ([alpha] = x, y, z) is the finite-difference of order N that is used to approximate and replace the spatial differential operators in matrix [A] and [B]. By using the Taylor central finite-difference method in [33], the following approximation of the spatial derivative of the field component V in the [alpha] ([alpha] = x, y, z) direction is defined as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (10)

and V = [E.sub.[alpha]], [H.sub.[alpha]], [alpha] = x, y, z, [d.sub.(2l-1)/2] = -[d.sub.-(2l-1)/2] = [[(-1).sup.l-1][(2L-1)!!.sup.2]/[2.sup.2L-2](L + l - 1)!(L-l)![(2l-1).sup.2]], where L is equals to different integers (L = 1,2,3, ...). N = 2L is the order of the central finite difference approximation, coefficients d depend on the order of the approximation N, which are listed in Table 1 for N = 2, 4, 6, 8, and 10.

For instance, for the second-order approximation (N = 2), for components of [E.sub.x] and [H.sub.z] in sub-step 1, after a series of manipulation, can be expressed as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (11a)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (11b)

(11a) is a linear system of equations with a diagonally banded coefficient matrix of three non-zero elements; it can be solved efficiently with special numerical packages and it has a computational complexity of O (M). M is the total number of unknowns of the system. Note that when an order of N is applied, the diagonally banded coefficient matrix will have N +1 non-zero elements. (11b) is an explicit equation that can be computed directly. For other field components, equations and computational processes similar to (11a)-(11b) can be obtained. The formulations for the other high-order approximation and other field components can be found in a similar manner by using Table 1.

3. NUMERICAL STABILITY ANALYSIS

By using the Fourier method, assuming [k.sub.x], [k.sub.y], and [k.sub.z] to be the spatial frequencies along the x, y, and z directions, the field components in spectral domain at the nth time step can be denoted as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (12)

By substituting (12) into (9a)-(9d), the following equations can be generated

[U.sup.n+1] = [[[LAMBDA]].sub.4][[[LAMBDA]].sub.3][[[LAMBDA]].sub.2][[[LAMBDA]].sub.1][U.sup.n] = [[LAMBDA]][U.sup.n]. (13)

where [[LAMBDA]] is the growth matrix, [[[LAMBDA]].sub.1], [[[LAMBDA]].sub.2], [[[LAMBDA]].sub.3], and [[[LAMBDA]].sub.4] are the growth matrices of the sub-steps. Moreover, [[[LAMBDA]].sub.3] = [[[LAMBDA]].sub.1], [[[LAMBDA]].sub.4] = [[[LAMBDA]].sub.2], and b = [DELTA]t/(2[epsilon]), d = [DELTA]t/(2/[mu]), [P.sub.[alpha]] = -2 sin([k.sub.[alpha]][DELTA][alpha])/[DELTA][alpha], [A.sub.[alpha]] = 4 + bd[P.sup.2.sub.[alpha]] ([alpha] = x, y, z).

By using Maple 9.0, the eigenvalues of [[LAMBDA]] can be found, as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (14)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (15)

where [xi] = R/S, and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (16a)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (16b)

It is obvious that the values of [P.sub.[alpha]] corresponding to the second-, the fourth-, the sixth-, the eighth-, and the tenth-order central finite-difference schemes are all real numbers, thus the eigenvalues associated with these schemes can be represented as (15). Since [absolute value of [[lambda].sub.1]] = [absolute value of [[lambda].sub.2]] = [absolute value of [[lambda].sub.3]] = [absolute value of [[lambda].sub.4]] = [absolute value of [[lambda].sub.5]] = [absolute value of [[lambda].sub.6]] = 1, we can conclude that all the high-order schemes are unconditionally stable.

4. NUMERICAL DISPERSION ANALYSIS

The numerical dispersion of the proposed methods can be found by following the procedure described in [14].

Assume the field to be a monochromatic wave with angular frequency [omega]

[U.sup.n+1] = [e.sup.j[omega][DELTA]t][U.sup.n]. (17)

Then, (13) can be expressed as

([e.sup.j[omega][DELTA]t][I] - [[LAMBDA]]) [U.sup.n] = 0. (18)

where [U.sup.n] is related to the initial field vector [U.sup.0] and defined by

[U.sup.n] = [U.sup.0][e.sup.j[omega][DELTA]tn]. (19)

For a nontrivial solution of (18), the determinant of the coefficient matrix should be zero as follows

det ([e.sup.j[omega][DELTA]t][I] - [[LAMBDA]]) = 0. (20)

With reference to the eigenvalues of [[LAMBDA]] above, the dispersion relationship of the high-order schemes can be deduced in (21).

[tan.sup.2]([omega][DELTA]t/2) = (S - R)/(S + R). (21)

Based on the previous arguments of the generalized eigenvalues, the dispersion relation in (21) can also be generalized to all the high-order schemes. This is achieved simply by replacing [P.sub.[alpha]] with those corresponding to the high-order difference schemes.

Assume that a wave propagating at angle [theta] and [phi] is in the spherical coordinate system. Then, [k.sub.x] = k sin [theta] sin [phi], [k.sub.y] = k sin [theta] cos [phi], [k.sub.z] = k cos [theta]. By substituting them into the dispersion relation (21), the numerical phase velocity [[??].sub.p] = [omega]/[??] can be solved numerically, where [??] is the numerical wave number. Moreover, several notations are introduced for clarity.

4.1. The Normalized Numerical Phase Velocity Error (NNPVE)

With the above definition of the numerical phase velocity, NNPVE at a propagation angle [theta] and [phi] can be defined as

NNPVE([theta], [phi]) = [absolute value of [[??].sub.p]([theta], [phi])/[upsilon] - 1] x 100%. (22)

4.2. The Maximum NNPVE

Here, in the entire range of [theta] and [phi], the maximum value of the NNPVE is denoted as the maximum NNPVE, i.e.,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (23)

4.3. The Normalized Numerical Phase Velocity Anisotropic Error (NNPVAE)

In the entire range of [phi], the normalized numerical phase velocity anisotropic error (NNPVAE) can be defined as,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (24)

4.4. The Maximum NNPVAE

In the entire range of [phi], the maximum value of the NNPVAE is denoted as the maximum NNPVAE, i.e.,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (25)

4.5. The Cell Per Wavelength (CPW)

For simplicity, uniform cells are considered here ([DELTA]x = [DELTA]y = [DELTA]z), CPW is introduced as the number of the cell per wavelength, i.e.,

CPW = [lambda]/[DELTA]x. (26)

where [lambda] is the wavelength with no numerical anisotropy.

4.6. The CFL Number (CFLN)

CFLN is defined as the ratio between the time step taken and the maximum CFL limit of the FDTD method, which can be expressed as

CFLN = [DELTA]t/[DELTA][t.sub.CFL_FDTD]. (27)

5. NUMERICAL RESULTS AND DISCUSSION

In this section, in order to demonstrate the efficiency of the proposed methods, numerical experiments are presented first, and then numerical dispersion characteristics of the second-, the fourth-, the sixth-, the eighth-, and the tenth-order four-step ADI-FDTD methods are presented, which based on the equations derived in the last sections. This study will be useful for the selection and evaluation of various higher order four-step ADI-FDTD methods.

5.1. Numerical Verifications of the Unconditional Stability

Although the unconditional stability has been theoretically proved for the proposed methods, numerical experiments are needed to confirm the theoretical result.

The FDTD method, the ADI-FDTD method, and the four-step ADI-FDTD method are utilized to simulate a cavity of 9 mm x 6 mm x 15 mm in size, all of which have the second-order accuracy in space. In addition, the cavity is filled with air and terminated with perfect electric conducting (PEC) boundaries. Moreover, a sinusoidal modulated Gaussian pulse of exp[-[(t - [t.sub.0]).sup.2]/[T.sup.2]] x sin[2[pi][f.sub.0](t - [t.sub.0])] is used as the excitation source, where T = 30 ps, [t.sub.0] = 3xT, [f.sub.0] = 20 GHz. The mesh size is chosen as [DELTA]x = [DELTA]y = [DELTA]z = 0.30 mm with 30 samples per wavelength at the highest frequency of the excitation source, leading to a mesh number of 30 x 20 x 50. [DELTA][t.sub.CFL_FDTD] = 0.57735 ps is the maximum time step size to satisfy the limitation of the 3D CFL condition in the conventional FDTD method. For the conventional FDTD method, CFLN = 1 and the step number is 6000, and the total simulation time is selected to be 3.4641ns. The simulations are performed on a computer of Pentium IV with 4GB RAM, and the computer program is developed with C++.

Figure 1 shows the [E.sub.z]-field at the observation point for the ADI-FDTD and four-step ADI-FDTD methods with CFLN = 6. From Fig. 1, the result of the four-step ADI-FDTD method is in better agreement with the FDTD method than that of the ADI-FDTD method. It can be seen that the proposed method did provide stable solutions even when the time step exceeds the CFL by six times. Further experiments (not shown here due to limitation of space) found that the solutions were bounded with even larger time steps and with higher orders. This verifies the unconditional stability of the proposed methods.

Table 2 shows the comparisons of results of three methods. As can be seen from Table 2, the four-step ADI-FDTD method with 0.3 mm mesh size and CFLN = 6 provides the same level of accuracy as the ADI-FDTD method with 0.3 mm mesh size and CFLN = 3. In addition, the CPU time and the memory requirement for two methods are the same. Therefore, the four-step ADI-FDTD method can use larger CFLN with the same level of the accuracy. Furthermore, the four-step ADI-FDTD method with CFLN = 3 is more accurate than that of the ADI-FDTD method with CFLN = 3. In addition, the four-step ADI-FDTD method with 0.3 mm mesh size provides the same accuracy as the ADI-FDTD method with 0.15 mm mesh size. In addition, the four-step ADI-FDTD scheme requires the CPU time of 52 s and the memory requirement of 3.9136 MB. However, the ADI-FDTD method increases the CPU time to 307 s and the memory requirement to 31.2995 MB, respectively. Consequently, with the same level of accuracy, the saving in the CPU time and the memory requirement of the four-step ADI-FDTD method can be more than 83.06% and 87.5% in comparisons with the ADI-FDTD method.

On the other hand, when the proposed fourth-order scheme works, the material interface and boundary conditions are treated using the fourth-order method in order to maintain the overall accuracy. To improve the computational efficiency, optimized FDTD methods based on the (2, 4) stencil were proposed in [34]. In the optimized (2, 4) FDTD methods, the spatial differential operators are approximated by the cell-centered finite difference scheme with four stencils, which features second-order accuracy in general. Subsequently, the optimized method in [34] was extended into the ADI-FDTD method and further reduced the dispersion errors by using the (2, 4) stencil approach [35]. Recently, to further reduce the dispersion error of the LOD-FDTD method, the parameter optimized method to the fourth-order LOD-FDTD was proposed based on the (2, 4) stencil [36]. Therefore, the similar approach will be used to the fourth-order four-step ADI-FDTD method. It is meaningful to propose the parameter optimized method for the four-step ADI-FDTD method based on the (2, 4) stencil, and it will be included in our future research plan.

5.2. Comparisons of Results with Second-order FDTD Methods

Figure 2 shows the NNPVE versus [phi] with CFLN = 6 and CPW = 30. From Fig. 2, it is apparent that the NNPVE of the proposed method is greatly reduced compared with the ADI-FDTD method. For instance, with [theta] = 45[degrees] and [phi] = 45[degrees], the NNPVE of the proposed method is reduced by more than 67% in comparison with the ADI-FDTD method. Moreover, the NNPVE reaches minimum at [theta] = 45[degrees], [phi] = 45[degrees].

Figures 3 and 4 present the NNPVE and the NNPVAE versus [phi] with CFLN = 3, 6, CPW = 30, respectively. As can be seen from Fig. 3, the NNPVE increases when CFLN increases. However, the increase of the NNPVE of the proposed method is much less pronounced than that of the ADI-FDTD method. In addition, the proposed method with CFLN = 6 has the same NNPVE value compared with the ADI-FDTD method with CFLN = 3. On the other hand, From Fig. 4, the NNPVAE reaches maximum at [phi] = 45[degrees]. Other characteristics are similar to the NNPVE in Fig. 3, which are not shown here again.

Figures 5 and 6 show the maximum NNPVE and the maximum NNPVAE versus CFLN with CPW = 60, 120, respectively. From Fig. 5, with CPW = 60 and CFLN = 20, the maximum NNPVE of the ADI-FDTD method is 10%. However, the maximum NNPVE of the proposed method is 3%, which is lower than that of the ADI-FDTD method. Furthermore, the maximum NNPVE increases as CFLN increases. The proposed method with CPW = 60 has the same error compared with the ADI-FDTD method with CPW = 120. On the other hand, from Fig. 6, it can be seen that Fig. 6 is similar to Fig. 5, except that the value of the maximum NNPVAE is different from the value of the maximum NNPVE. Such an improvement of the accuracy, which is realized by the proposed method with the coarse mesh, leads to other advantages, such as higher computational efficiency and lower memory requirements.

Figures 7 and 8 show the maximum NNPVE and the maximum NNPVAE versus CPW with CFLN = 3, 6, respectively. As can be seen from Figs. 7 and 8, the maximum NNPVE and the maximum NNPVAE have the similar characteristics, except that the value of the maximum NNPVE is different from the maximum NNPVAE. Concretely, from Fig. 7, the maximum NNPVE decreases as CPW increases. For the same CFLN value, the maximum NNPVE of the proposed method is lower than that of the ADI-FDTD method. Moreover, the proposed method with CFLN = 6 has the same Maximum NNPVE value compared with the ADI-FDTD method with CFLN = 3.

5.3. Effect of Higher Order Schemes on Numerical Dispersion

Figures 9 and 10 present the NNPVE versus [phi] with [theta] = 45[degrees] and the NNPVAE versus [phi] for the proposed method of different orders, respectively. It is apparent that the NNPVE is reduced when a high-order approximation is used. However, the amount of the reduction is related to the CPW value. For instance, with the fine mesh (CPW = 30), the NNPVE of the fourth-order scheme is reduced by about 75% of that of the second-order scheme. However, when the order is increased beyond four, the NNPVE decreases little (see Figs. 9(a) and (b)). On the other hand, with the coarse mesh (CPW = 5), the NNPVE can be decreased with the higher order schemes, which can be seen from Fig. 9(c). Nevertheless, the reduction in the NNPVE is not linearly proportional to the orders. The reduction starts to level off beyond a certain order. On the other hand, from Fig. 10, the NNPVAE has the similar phenomenon to the NNPVE in Fig. 9, except that the NNPVAE in Fig. 10 reaches maximum at [phi] = 45[degrees] and minimum at [phi] = 0[degrees], 90[degrees]. However, the NNPVE in Fig. 9 reaches maximum at [phi] = 0[degrees], 90[degrees] and minimum at [phi] = 45[degrees].

5.4. Effect of Propagation Angle ([phi]) on Numerical Dispersion

Figure 11 shows the NNPVE versus [phi] with CFLN = 10 and CPW = 60 for the proposed fourth-order scheme. It is seen that the NNPVE reaches minimum at the diagonal direction ([theta] = 45[degrees], [phi] = 45[degrees]) and maximum at three axial directions, ([theta] = 90[degrees], [phi] = 0[degrees]), ([theta] = 90[degrees], [phi] = 90[degrees]) and [theta] = 0[degrees]. This represents a numerical anisotropy that is similar to that of the conventional FDTD method. In addition, even with CFLN = 10 (large time step), the maximum NNPVE is less than 0.75%.

5.5. Effect of the CFLN Value on Numerical Dispersion

Figure 12 shows how the CFLN value affects the NNPVE. From Fig. 12, the NNPVE is less than 0.2% with CFLN = 3 while is less than 0.7% with CFLN = 6. In other words, the NNPVE increases when CFLN increases. This is consistent with the fact the Taylor's series approximation to obtain (8) is accurate under the assumption of small time step [DELTA]t. The smaller the time step [DELTA]t, the better accuracy of (8). It is concluded that even with a high order, the time step [DELTA]t should not be too large.

5.6. Effect of the CPW Value on Numerical Dispersion

Figures 13 and 14 present how the CPW value affects the NNPVE and the NNPVAE for both the second- and the fourth-order schemes, respectively. It can be seen from Figs. 13 and 14, the NNPVE and the NNPVAE have the similar characteristics, except that the NNPVE reaches minimum at [phi] = 45[degrees] and the NNPVAE reaches maximum at [phi] = 45[degrees]. Moreover, from Fig. 13, for the second-order scheme, when the CPW value increases from 30 to 60, the NNPVE is reduced by four time, and when the CPW value increases from 30 to 120, the NNPVE is reduced by twelve times. It can be seen that the NNPVE is decreased by refining the mesh for both the second- and the fourth-order schemes. This is explained that as CPW increases, the mesh sizes [DELTA]x, [DELTA]y and [DELTA]z become smaller, the approximation of (10) is more accurate.

Figures 15 and 16 present the maximum NNPVE and the maximum NNPVAE versus CPW for the proposed high-order schemes with CFLN = 1, 6, respectively. From Figs. 15 and 16, the maximum NNPVE and the maximum NNPVAE have the similar characteristics. In detail, when CFLN is smaller, it is apparent that the increase of the order is more effective than the increase of CPW. In other words, the higher order scheme with the coarsest mesh may result in higher computation efficiency than the lower order scheme with the finest mesh for the same required accuracy. At the same time, the memory requirement by the updating procedures may also be reduced. On the other hand, when CFLN is larger, it presents an extreme case that the increase of the order reduces a little in the maximum NNPVE and the maximum NNPVAE, whereas the increase of CPW does reduce the maximum NNPVE and the maximum NNPVAE a lot, up to four times. Nevertheless, increasing CPW leads to such disadvantage as higher memory requirements.

From the above investigations, we can conclude that higher order schemes can achieve lower NNPVE and NNPVAE; however, the improvement is related to the propagation angle, the time step and the mesh size. When the time step is increased, the NNPVE and the NNPVAE of each higher order scheme becomes larger, and thus the scheme to decrease the effect of temp oral discretization on numerical dispersion will also b e needed.

6. CONCLUSION

High-order unconditionally-stable four-step ADI-FDTD methods in 3-D domains have been presented. Based on the exponential evolution operator (EEO), the Maxwell's equations have been split into four sub-procedures, and accordingly, the time step is divided into four sub-steps. Moreover, the high-order central finite-difference schemes based on the Taylor central finite-difference method have been used to approximate the spatial differential operators. The unconditional stability of the proposed methods has been theoretically proven and numerical verification by the experiments. Furthermore, the generalized form of the dispersion relations of the proposed methods has also been provided. To trade off between the computational efficiency and accuracy, the effects of the order of schemes, the propagation angle, the time step, and the mesh size on the dispersion characteristics have been illustrated through the numerical results.

The improved accuracy realized by the higher order four-step ADI-FDTD schemes actually leads to other advantages of them, which are higher computation efficiency and lower memory requirement, since they may obtain better accuracy even with coarser mesh. Therefore, the proposed schemes will be of interest and usefulness in electromagnetics problems and can be applied into waveguide, antenna or EMC problems. Generalizing these extensions will be our future work.

ACKNOWLEDGMENT

This work was supported by the National Natural Science Foundation of China (61171029 and 61201036), the Natural Science Foundation of Guangdong Province (S2012040006893) and the Program of State Key Laboratory of Millimeter Waves (K201102), Southeast of University.

REFERENCES

[1.] Yee, K. S., "Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media," IEEE Trans. on Antennas on Propag., Vol. 14, No. 3, 302-307, May 1966.

[2.] Su, D., D. M. Fu, and Z. H. Chen, "Numerical modeling of active devices characterized by measured S-parameters in FDTD," Progress In Electromagnetics Research, Vol. 80, 381-392, 2008.

[3.] Li, J., L. X. Guo, and H. Zeng, "FDTD method investigation on the polarimetric scattering from 2-D rough surface," Progress In Electromagnetics Research, Vol. 101, 173-188, 2010.

[4.] Izadi, M., M. Z. A. A. Kadir, and C. Gomes, "Evaluation of electromagnetic fields associated with inclined lightning channel using second order FDTD-hybrid methods," Progress In Electromagnetics Research, Vol. 117, 209-236, 2011.

[5.] Koh, I. S., H. Kim, J. M. Lee, J. G. Yook, and C. S. Pil, "Novel explicit 2-D FDTD scheme with isotropic dispersion and enhanced stability," IEEE Trans. on Antennas on Propag., Vol. 54, No. 11, 3505-3510, Nov. 2006.

[6.] Wang, C. C. and C. W. Kuo, "An efficient scheme for processing arbitrary lumped multiport devices in the finite-difference time-domain method," IEEE Trans. on Microw. Theory and Tech., Vol. 55, No. 5, 958-965, May 2007.

[7.] Zygiridis, T. T. and T. D. Tsiboukis, "Improved finite-difference time-domain algorithm based on error control for lossy materials," IEEE Trans. on Microw. Theory and Tech., Vol. 56, No. 6, 1440-1445, Jun. 2008.

[8.] Tofighi, M. R., "FDTD modeling of biological tissues Cole-Cole dispersion for 0.5-30 GHz using relaxation time distribution samples-novel and improved implementations," IEEE Trans. on Microw. Theory and Tech., Vol. 57, No. 10, 2588-2596, Oct. 2009.

[9.] Kim, H., I. S. Koh, and J. G. Yook, "Enhanced total-field/scattered-field technique for isotropic-dispersion FDTD scheme," IEEE Trans. on Antennas on Propag., Vol. 58, No. 10, 3407-3411, Oct. 2010.

[10.] Taflove, A. and S. C. Hagness, Computational Electrodynamics: The Finite-difference Time-domain Method, 2nd Edition, Artech House, Boston, MA, 2000.

[11.] Namiki, T., "A new FDTD algorithm based on alternating-direction implicit method," IEEE Trans. on Microw. Theory and Tech., Vol. 47, No. 10, 2003-2007, Oct. 1999.

[12.] Zheng, F., Z. Chen, and J. Zhang, "Toward the development of a three-dimensional unconditionally stable finite-difference time-domain method," IEEE Trans. on Microw. Theory and Tech., Vol. 48, No. 9, 1550-1558, Sep. 2000.

[13.] Sun, G. L. and C. W. Trueman, "Analysis and numerical experiments on the numerical dispersion of two-dimensional ADI-FDTD," IEEE Antennas Wireless Propag. Lett., Vol. 2, No. 1, 78-81, 2003.

[14.] Zheng, F. and Z. Chen, "Numerical dispersion analysis of the unconditionally stable 3-D ADI-FDTD method," IEEE Trans. on Microw. Theory and Tech., Vol. 49, No. 5, 1006-1009, May 2001.

[15.] Wang, M. H., Z. Wang, and J. Chen, "A parameter optimized ADI-FDTD method," IEEE Antennas Wireless Propag. Lett., Vol. 2, No. 1, 118-121, 2003.

[16.] Ahmed, I. and Z. Chen, "Error reduced ADI-FDTD methods," IEEE Antennas Wireless Propag. Lett., Vol. 4, 323-325, 2005.

[17.] Zheng, H. X. and K. W. Leung, "An efficient method to reduce the numerical dispersion in the ADI-FDTD," IEEE Trans. on Microw. Theory and Tech., Vol. 53, No. 7, 2295-2301, Jul. 2005.

[18.] Zhang, Y., S. W. Lu, and J. Zhang, "Reduction of numerical dispersion of 3-D higher order alternating-direction-implicit finite-difference time-domain method with artificial anisotropy," IEEE Trans. on Microw. Theory and Tech., Vol. 57, No. 10, 2416-2428, Oct. 2009.

[19.] Kong, K. B., S. O. Park, and J. S. Kim, "Stability and numerical dispersion of 3-D simplified sampling biorthogonal ADI method," Journal of Electromagnetic Waves and Application, Vol. 24, No. 1, 1-12, 2010.

[20.] Sun, G. and C. W. Trueman, "Approximate Crank-Nicolson schemes for the 2-D finite-difference time-domain method for TEz waves," IEEE Trans. on Antennas on Propag., Vol. 52, No. 11, 2963-2972, Nov. 2004.

[21.] Xu, K., Z. H. Fan, D. Z. Ding, and R. S. Chen, "GPU accelerated unconditionally stable Crank-Nicolson FDTD method for the analysis of three-dimensional microwave circuits," Progress In Electromagnetics Research, Vol. 102, 381-395, 2010.

[22.] Fu, W. and E. L. Tan, "Development of split-step FDTD method with higher-order spatial accuracy," Electron. Lett., Vol. 40, No. 20, 1252-1253, Sep. 2004.

[23.] Xiao, F., X. H. Tang, L. Guo, and T. Wu, "High-order accurate split-step FDTD method for solution of Maxwell's equations," Electron. Lett., Vol. 43, No. 2, 72-73, Jan. 2007.

[24.] Chu, Q. X. and Y. D. Kong, "Three new unconditionally-stable FDTD methods with high-order accuracy," IEEE Trans. on Antennas on Propag., Vol. 57, No. 9, 2675-2682, Sep. 2009.

[25.] Kong, Y. D. and Q. X. Chu, "Reduction of numerical dispersion of the six-stages split-step unconditionally-stable FDTD method with controlling parameters," Progress In Electromagnetics Research, Vol. 122, 175-196, 2012.

[26.] Shibayama, J., M. Muraki, J. Yamauchi, and H. Nakano, "Efficient implicit FDTD algorithm based on locally one-dimensional scheme," Electron. Lett., Vol. 41, No. 19, 1046-1047, Sep. 2005.

[27.] Ahmed, I., E. K. Chua, E. P. Li, and Z. Chen, "Development of the three-dimensional unconditionally-stable LOD-FDTD method," IEEE Trans. on Antennas on Propag., Vol. 56, No. 11, 3596-3600, Nov. 2008.

[28.] Wang, Z., J. Chen, and Y. Chen, "Development of a higher-order ADI-FDTD method," Microwave Optical Technol. Lett., Vol. 37, No. 1, 8-12, Apr. 2003.

[29.] Fu, W. and E. L. Tan, "Stability and dispersion analysis for higher order 3-D ADI-FDTD method," IEEE Trans. on Antennas on Propag., Vol. 53, No. 11, 3691-3696, Nov. 2005.

[30.] Liu, Q. F., Z. Chen, and W. Y. Yin, "An arbitrary order LOD-FDTD method and its stability and numerical dispersion," IEEE Trans. on Antennas on Propag., Vol. 57, No. 8, 2409-2417, Aug. 2009.

[31.] Kong, Y. D. and Q. X. Chu, "High-order split-step unconditionally-stable FDTD methods and numerical analysis," IEEE Trans. on Antennas on Propag., Vol. 59, No. 9, 3280-3289, Sep. 2011.

[32.] Yang, S. C., Z. Chen, Y. Q. Yu, and W. Y. Yin, "An unconditionally stable one-step arbitrary-order leapfrog ADI-FDTD method and its numerical properties," IEEE Trans. on Antennas on Propag., Vol. 60, No. 4, 1995-2003, Apr. 2012.

[33.] Xiao, F., X. H. Tang, and H. Ma, "High-order US-FDTD based on the weighted finite-difference method," Microwave Optical Technol. Lett., Vol. 45, No. 2, 142-144, Apr. 2005.

[34.] Sun, G. and C. W. Trueman, "Optimized finite-difference time-domain methods on the (2, 4) stencil," IEEE Trans. on Microw. Theory and Tech., Vol. 53, No. 3, 832-842, Mar. 2005.

[35.] Fu, W. and E. L. Tan, "A parameter optimized ADI-FDTD method based on the (2, 4) stencil," IEEE Trans. on Antennas on Propag., Vol. 54, No. 6, 1836-1842, Jun. 2006.

[36.] Liu, Q. F., W. Y. Yin, Z. Chen, and P. G. Liu, "An efficient method to reduce the numerical dispersion in the LOD-FDTD method based on the (2, 4) stencil," IEEE Trans. on Antennas on Propag., Vol. 58, No. 7, 2384-2393, Jul. 2010.

Yongdan Kong (1), *, Qingxin Chu (1,2), and Ronglin Li (1)

(1) School of Electronic and Information Engineering, South China University of Technology, Guangzhou, Guangdong 510640, China

(2) The State Key Laboratory of Millimeter Waves, Nanjing, Jiangsu 210096, China

Received 22 October 2012, Accepted 21 December 2012, Scheduled 10 January 2013

* Corresponding author: Qingxin Chu (qxchu@scut.edu.cn).

Table 1. Coefficients of high-order central finite-difference schemes. N/d [d.sub.1/2] [d.sub.3/2] [d.sub.5/2] 2 1 4 9/8 -1/24 6 225/192 -25/384 3/640 8 1.196289 -0.079753 0.009570 10 1.211243 -0.089722 0.013843 N/d [d.sub.7/2] [d.sub.9/2] 2 4 6 8 -0.000698 10 -0.001766 0.000119 Table 2. Comparisons of results with three methods. Scheme Mesh CFLN Step Result (GHz) Relative size (mm) number [TE.sub.101] error (%) (19.437) FDTD 0.30 1 6000 19.43 0.0360 0.15 1 12000 19.43 0.0360 ADI-FDTD 0.30 3 2000 19.37 0.3447 0.30 6 1000 19.22 1.1164 0.15 6 2000 19.38 0.2933 Four-step 0.30 3 2000 19.41 0.1389 ADI-FDTD 0.30 6 1000 19.37 0.3447 0.15 6 2000 19.42 0.0875 Scheme Mesh Result (GHz) Relative Result (GHz) size (mm) [TE.sub.101] error (%) [TM.sub.111] (26.926) (30.046) FDTD 0.30 26.91 0.0594 30.03 0.15 26.92 0.0233 30.04 ADI-FDTD 0.30 26.74 0.6908 29.84 0.30 26.29 2.3620 29.32 0.15 26.76 0.6165 29.86 Four-step 0.30 26.86 0.2451 29.98 ADI-FDTD 0.30 26.74 0.6908 29.84 0.15 26.88 0.1708 29.99 Scheme Mesh Relative CPU Memory size (mm) error (%) time (s) (MB) FDTD 0.30 0.0533 19 2.0598 0.15 0.02 274 16.4784 ADI-FDTD 0.30 0.6856 26 3.9136 0.30 2.4163 13 3.9136 0.15 0.6191 307 31.2995 Four-step 0.30 0.2197 52 3.9136 ADI-FDTD 0.30 0.6856 26 3.9136 0.15 0.1864 617 31.2995

Printer friendly Cite/link Email Feedback | |

Author: | Kong, Yongdan; Chu, Qingxin; Li, Ronglin |
---|---|

Publication: | Progress In Electromagnetics Research |

Article Type: | Report |

Geographic Code: | 7ISRA |

Date: | Mar 1, 2013 |

Words: | 5800 |

Previous Article: | Scattering of an axial Gaussian beam by a conducting spheroid with non-confocal chiral coating. |

Next Article: | Plane wave beam produced by an exclusive medium. |

Topics: |