International Journal of Antennas and Propagation
Volume 2012 (2012), Article ID 956431, 8 pages
http://dx.doi.org/10.1155/2012/956431
Research Article

## Optimized Operator-Splitting Methods in Numerical Integration of Maxwell's Equations

1Key Lab of Intelligent Computing & Signal Processing, Anhui University, Hefei 230039, China
2Department of Electrical and Electronic Engineering, University of Hong Kong, Hong Kong

Received 11 March 2011; Accepted 27 April 2011

Academic Editor: Ning Yuan

Copyright © 2012 Z. X. Huang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

Optimized operator splitting methods for numerical integration of the time domain Maxwell's equations in computational electromagnetics (CEM) are proposed for the first time. The methods are based on splitting the time domain evolution operator of Maxwell's equations into suboperators, and corresponding time coefficients are obtained by reducing the norm of truncation terms to a minimum. The general high-order staggered finite difference is introduced for discretizing the three-dimensional curl operator in the spatial domain. The detail of the schemes and explicit iterated formulas are also included. Furthermore, new high-order Padé approximations are adopted to improve the efficiency of the proposed methods. Theoretical proof of the stability is also included. Numerical results are presented to demonstrate the effectiveness and efficiency of the schemes. It is found that the optimized schemes with coarse discretized grid and large Courant-Friedrichs-Lewy (CFL) number can obtain satisfactory numerical results, which in turn proves to be a promising method, with advantages of high accuracy, low computational resources and facility of large domain and long-time simulation. In addition, due to the generality, our optimized schemes can be extended to other science and engineering areas directly.

#### 1. Introduction

The finite-difference time-domain (FDTD(2, 2)) method [1, 2] has been widely used to simulate the transient solutions of electromagnetic problems involving the analysis and design of microwave structures, many other engineering applications, and the electromagnetic wave propagation in various media. Despite its simplicity and modeling versatility; however, FDTD(2, 2) is very computationally intensive due to its two inherent physical constraints, one being the numerical dispersion and another being the numerical stability. These limitations have always made it a matter of great interest to improve the efficiency of FDTD(2, 2) scheme and have led researchers to the development of various new schemes.

To improve the numerical dispersion, some high-order space strategies have been put forward. For example, Fang proposed the high-order FDTD(4, 4) method [3]. Yet, the method is difficult to handle material interface for modeling the complex three-dimensional objects. Another approach is the staggered FDTD(2, 4) method [46]. However, the method must set lower Courant-Friedrichs-Levy (CFL) number to obtain high-order numerical precision. In order to further explore efficient methods for optimum electromagnetic simulation, new improved time strategies referred as the high-order Runge-Kutta (R-K) approach was introduced in [7, 8]. However, the approach is dissipative and needs large amount of memory. Other alternative method is the alternating direction implicit FDTD (ADI-FDTD) algorithm [9, 10]. Although it saves CPU time owing to unconditional stability, undesirable numerical precision and dispersion will happen once the high CFL number is adopted. Another systematic approach to solve the time-dependent Maxwell equations with unconditionally stable numerical schemes was proposed and developed [11]. The basic idea of the methods is to employ a Lie-Trotter-Suzuki product formula to approximate the time evolution operator. As shown in [12], most of these methods can be seen as special cases of the time-operator-splitting methods where the original time evolution operator is split into a number of suboperators. There are two good reasons for this. Firstly, many such suboperators are simple and easy to implement. Secondly, splitting methods can preserve important mathematical and physical properties of the original system. Now methods of operator-splitting have been widely used and considered in various applications in science and engineering—from the evolution of techniques for solving linear equations that arise in reservoir simulation, astrophysical, and bioengineering applications to tsunami modeling and furthermore.

In this paper, particularly, we consider optimized operator-splitting methods for numerical solution of the time-dependent Maxwell’s equations. In Section 2, the time-domain Maxwell’s equations are rewritten as a time evolution matrix operator form. The novel approach is proposed to improve the efficiency of the time evolution operator based on splitting it into suboperator, and optimal time coefficients are obtained by reducing the norm of truncation terms to a minimum. The general high-order staggered finite difference is introduced for discretizing the three-dimension curl operator in the spatial domain. The explicit discretized formulas are presented in Section 3. New high-order Padé approximations and the stability analysis are also included. Section 4 presents numerical examples, and the conclusions are made in Section 5.

#### 2. General Formulations of Splitting Schemes

##### 2.1. The Unique Solution of Time Domain Maxwell Equations in Matrix Form

Maxwell’s equations in an isotropic medium can be rewritten in a matrix form as where and are the permeability and permittivity, and are electric and magnetic conductivities, and are zeros matrix and identity matrix, is matrix representing three-dimensional curl operator. Equation (1) can be cast in the following compact form Here, is the full electromagnetic field variable. Although only the time dependence is written explicitly, all these quantities additionally depend on space, but for simplicity of notation, we will omit the spatial dependence. If an initial configuration Z(0) is provided, the unique solution to (2) can be presented as where denotes the time step, and the operator has been split into two suboperators. The significance of such splitting will be understood below.

##### 2.2. The Splitting Method for Exponential Propagator in the Time Domain

Note, solution (3) is quite formal because the exponential propagator does not allow to be evaluated exactly at any given . However, at small enough values of , the total propagator can be split using the formula The coefficients   and   in this formula should be chosen in such a way to provide the highest possible value for   at a given integer number . Here is stage number needed in every integer time step, and is the order of the approximation. The main advantage of the above splitting is that the time-reversibility of solutions (following from the property of time evolution operator ) can also be reproduced by imposing additional time reversible constraints on the coefficients, namely, , and with . Note also that the splitting method is quite general to build numerical integrators of arbitrary orders   with stages. But we cannot choose the stage   too large, because this results in a too large number, namely , of expensive force every time step. In this paper, particularly, we choose and minimize the truncation errors to significantly with a little additional computational cost.

For , the extended splitting result can be represented in the form Note also, the propagators can be calculated analytically as follows: Here, the time reversible coefficients and the condition have already been taken into account. With these assumptions, we actually only have three coefficients, namely, to be determined. Using the Baker-Campbell-Hausdorff (BCH) formula, the explicit expression for and can be expressed as Here, and . The detailed expression of can be founded in [13]. The formula (5) represents a fourth-order scheme at , that is, and equal to zero. Now, we have three coefficients and only two equations, we are free to set one more equation. In this paper, we choose to reduce the fifth-order truncation error terms to a minimum, and the coefficients can be obtained by solving the system of equations One possible set of the coefficients are c1 = 0.1786, c2 = −0.0066, d1 = 0.7123, consists global minimum truncation error 0.00093. Now, we transfer to the discretization of spatial domain.

##### 2.3. The Mth-Order Difference Approximation for the First-Order Spatial Partial Derivative Operators

Let approximates the exact solution at the point in the nth time step. The following staggered th-order space difference operators are used to approximate the first-order spatial partial derivative in the -direction, that is, in three-dimensional curl operator . So Here, parameters   for minimum truncation error in (9). Similarly, and are constructed in a similar manner and used to approximate and , respectively.

#### 3. Practical Implementation

##### 3.1. Explicit Discretization Formulas

When one uses coefficients and of order four and substitutes the space difference operators with th-order accuracy for the first-order partial difference operators in , the scheme is acquired.

For example, using Yee grid [1], the detailed expressions of and components in the scheme at the th stage calculation after the th time step are as follows: with

Here (i)we use Padé and Padé to approximate the expressions of and : Our new approximate acquired super stability and efficiency than the Padé(2, 2) method in [14] when increases, as indicated in Figure 1. In addition, our approximate can be not only used to treat interior dielectric medium and conductor, but also directly applied in Berenger’s perfectly matched layer (PML) absorbing boundary conditions [15].(ii)The averaged permittivity over the patch can be expressed as where and are the surfaces enclosed by curve-labeled single arrow and double arrow as indicated in Figure 2, α is free parameter, and we use in our following numerical examples in order tobeconsistent with the high-order difference in spatial direction. In addition, the averaged conductivity can be treated in a similar way, which consumes little CPU time at the initial process by refined subcell modeling.

Figure 1: Comparison between two kinds of approximation. (a) For and (b) .
Figure 2: Subcell technique specifies for the Ex field calculation in scheme.
##### 3.2. The Stability of the Proposed (4, M) Scheme

The conventional Fourier mode method is used to analyze the stability of the proposed scheme. For clarity, the discussion begins with a 1-D z-directed, x-polarized TEM wave. The equations can be rewritten as The field components in the th time-step are denoted as here, is the spatial frequency along the z-direction. Substituted (15) into (9), we obtain Then (14) can be rewritten as with

Applying the scheme to (17), the time-marching relation can be expressed as where Solving for the eigenvalues, we can obtain where We may conclude that if , and the scheme is stable. Moreover, under this condition, the scheme is non-dissipative.

In the three-dimensional (3D) case, the continuous-time discrete-space Maxwell’s equations can be written as Considering , (23) can be rewritten in tensor form as where is the tensor matrix defined by the spherical angles [16]. Using the similar analysis, the same formula as (21) is obtained in the case except that   Generally speaking, the stability of the scheme is assessed by the Courant-Friedrichs-Levy limit, that is, CFLmax, which can be expressed as where is the temporal stability factor, and the numerical time evolution operator will not blow up for all , which can be determined according to the following inequality here, is the spatial stability factor, which is defined as , and is the dimensional number. The numerical results of CFLmax with scheme, and FDTD(2, 2) are listed in Table 1. For comparison, we also plot the dispersion curves for (4, 4) and FDTD(2, 2) method in Figure 3. As we can see from Figure 3, the relative phase velocity error decreases with the increasing of sampled points per wavelength (PPW). As expected, (4, 4) scheme acquires better numerical dispersion than FDTD(2, 2) method.

Table 1: CFLmax and accuracy for different methods.
Figure 3: Numerical dispersion curves for (4, 4) scheme and FDTD(2, 2) method.

#### 4. Numerical Examples

Remember that we cannot choose the discretized order of to be too large, because this results in a too large number of expensive forces in real application. In our following numerical examples, we mainly concentrate on the performers of .

##### 4.1. Radiation of a Dipole

We considered a computational domain of 46 by 46 by 46 cells surrounded by a ten-point PML. A vertical dipole P was located at point (23,23,23), in the center of the domain. To test the efficiency of our new Padé approximate applied in PML, Figure 4 shows the Ez  field emanating along the plane of after 89 time steps. Notice that the part of the field not in the PML radiates concentrically from the source, as it should. In addition, the results for the vertical field Ez at point , two cells from the PML, with FDTD(2, 2) and (4, 4) scheme are summarized in Table 2. Note that for roughly the same computational cost, the (4, 4) scheme gives results that are more accurate than the FDTD(2, 2) method, which in turn proves to be a promising method, with advantages of high accuracy, low computational resources, and facility of large domain and long time simulation.

Table 2: Comparison of results for FDTD(2, 2) and (4, 4) scheme.
Figure 4: Radiation from a dipole with a ten-point PML. (a) The distribution of the Ez field. (b) The contour of the Ez amplitude field.
##### 4.2. Scattering of the Dielectric Sphere

Next, consider a dielectric sphere illuminated by a plane wave propagating in the z-direction and E polarized in the x-direction. The frequency of the incident wave is 300 MHz. The sphere has a diameter of 1.0 m, relative permittivity , and conductivity of zero. We use uniform grid Δx = Δy = Δz = Δ. The total computational domain is 80 by 80 by 80 cells, total field occupies 32 by 32 by 32 cells, and a ten-grid PMLs are used. We denote the relative radar cross section (RCS) error as where RCS* is the analytical solution, RCS is the solution with FDTD(2, 2), (4, 2), or (4, 4) scheme. Figure 5 shows the relative error of RCS computed with fourth-order accuracy in time domain and fourth-order, second-order accuracy in spatial domain, respectively.

Figure 5: The relative error of computed RCS with different scheme at CFL = 0.5. (a) E-plane (b) H-plane.

As we can see, the (4, 4) scheme is more accurate than the (4, 2) scheme under the same discretized grid and CFL number. When the grid enlarges to , the error of (4, 2) scheme becomes unacceptable. Figure 6 shows the comparisons between (4, 4) and (4, 2) schemes at . In this case, whatever we choose or , the results computed by (2, 2) scheme are divergent. The reason may be that CFL number exceeds the stability of (2, 2) scheme. The results for (4, 4) scheme are still acceptable except at some particular angles. Figure 7 shows the comparisons between (4, 4) and (2, 2) scheme with different discretized grid and CFL number. It is clear that with high CFL number and coarse grid, the results of the (4, 4) scheme are still acceptable to some extent. The time spent in (2, 2) method is longer, about 16 minutes, and it is about 12 minutes for (4, 4) scheme. The memory consumed is around 30 M for (2, 2) method, and it is about 20 M for (4, 4) scheme.

Figure 6: The relative error of computed E-plane RCS with different scheme at CFL = 0.6.
Figure 7: The relative error of computed E-plane RCS with (4, 4) and (2, 2) schemes.

As indicated in figures, we can come to a conclusion as follows.(i)The smaller apace discretized grid we fix, the higher numerical precision we obtain, no matter what scheme we adopt.(ii)With the same spatial discretizated scheme, the higher order of time domain discretized, the higher CFL number we get.(iii)The (4, 4) scheme with coarse discretized grid and high CFL number can reach satisfactory numerical results, which in turn proves to be a promising method, with advantages of high accuracy, low computational resources, and facility of large-domain and long-time simulation.

#### 5. Conclusions

We present optimized operator-splitting methods for numerical solution of the time-dependent Maxwell equations in the time domain. The general high-order staggered finite difference is introduced for approximating the three-dimensional curl operator in the spatial domain. The efficiency of the scheme, especially the (4, 4) scheme, has been verified by some numerical examples. The major shortcoming of the scheme is that it consumes more CPU time than the FDTD(2, 2) method when the same grid size is used. Effective parallel algorithm is an open question for further study.

#### Acknowledgments

This work is supported by the Key National Natural Science Foundation of China (no. 60931002) and Universities of Natural Science Foundation of Anhui Province (no. KJ2011A002).

#### References

1. K. S. Yee, “Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media,” IEEE Transactions on Antennas and Propagation, vol. 14, pp. 302–307, 1966.
2. A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method, Artech House, London, UK, 3rd edition, 2005.
3. J. Fang, Time domain finite difference computation for Maxwell's equations, Ph.D. thesis, University of California, Berkeley, Calif, USA, 1989.
4. A. Yefet and P. G. Petropoulos, “A staggered fourth-order accurate explicit finite difference scheme for the time-domain Maxwell's equations,” Journal of Computational Physics, vol. 168, no. 2, pp. 286–315, 2001.
5. S. V. Georgakopoulos, C. R. Birtcher, C. A. Balanis, and R. A. Renaut, “Higher-order finite-difference schemes for electromagnetic radiation, scattering, and penetration, part I: theory,” IEEE Antennas and Propagation Magazine, vol. 44, no. 1, pp. 134–142, 2002.
6. S. V. Georgakopoulos, C. R. Birtcher, C. A. Balanis, and R. A. Renaut, “Higher-order finite-difference schemes for electromagnetic radiation, scattering, and penetration, part 2: applications,” IEEE Antennas and Propagation Magazine, vol. 44, no. 2, pp. 92–101, 2002.
7. J. L. Young, D. Gaitonde, and J. S. Shang, “Toward the construction of a fourth-order difference scheme for transient EM wave simulation: staggered grid approach,” IEEE Transactions on Antennas and Propagation, vol. 45, no. 11, pp. 1573–1580, 1997.
8. J. S. Shang, “High-order compact-difference schemes for time-dependent Maxwell equations,” Journal of Computational Physics, vol. 153, no. 2, pp. 312–333, 1999.
9. T. Namiki, “A new FDTD algorithm based on alternating-direction implicit method,” IEEE Transactions on Microwave Theory and Techniques, vol. 47, no. 10, pp. 2003–2007, 1999.
10. F. Zheng and Z. Chen, “Toward the development of a three-dimensional unconditionally stable finite-difference time-domain method,” IEEE Transactions on Microwave Theory and Techniques, vol. 48, no. 9, pp. 1550–1558, 2000.
11. J. S. Kole, M. T. Figge, and H. De Raedt, “Unconditionally stable algorithms to solve the time-dependent Maxwell equations,” Physical Review E, vol. 64, no. 6, Article ID 066705, pp. 1–14, 2001.
12. M. A. Botchev, I. Faragó, and R. Horváth, “Application of operator splitting to the Maxwell equations including a source term,” Applied Numerical Mathematics, vol. 59, no. 3-4, pp. 522–541, 2009.
13. I. P. Omelyan, I. M. Mryglod, and R. Folk, “Optimized Forest-Ruth- and Suzuki-like algorithms for integration of motion in many-body systems,” Computer Physics Communications, vol. 146, no. 2, pp. 188–202, 2002.
14. T. Hirono, W. Lui, S. Seki, and Y. Yoshikuni, “A three-dimensional fourth-order finite-difference time-domain scheme using a symplectic integrator propagator,” IEEE Transactions on Microwave Theory and Techniques, vol. 49, no. 9, pp. 1640–1648, 2001.
15. J. P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” Journal of Computational Physics, vol. 114, no. 2, pp. 185–200, 1994.
16. W. C. Chew, Waves and Fields in Inhomogenous Media, Van Nostrand Reinhold, New York, NY, USA, 1990.