HighPerformance Computing and Engineering Applications in Electromagnetics
View this Special IssueResearch Article  Open Access
Z. X. Huang, X. L. Wu, W. E. I. Sha, B. Wu, "Optimized OperatorSplitting Methods in Numerical Integration of Maxwell's Equations", International Journal of Antennas and Propagation, vol. 2012, Article ID 956431, 8 pages, 2012. https://doi.org/10.1155/2012/956431
Optimized OperatorSplitting Methods in Numerical Integration of Maxwell's Equations
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 highorder staggered finite difference is introduced for discretizing the threedimensional curl operator in the spatial domain. The detail of the schemes and explicit iterated formulas are also included. Furthermore, new highorder 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 CourantFriedrichsLewy (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 longtime simulation. In addition, due to the generality, our optimized schemes can be extended to other science and engineering areas directly.
1. Introduction
The finitedifference timedomain (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 highorder space strategies have been put forward. For example, Fang proposed the highorder FDTD(4, 4) method [3]. Yet, the method is difficult to handle material interface for modeling the complex threedimensional objects. Another approach is the staggered FDTD(2, 4) method [4–6]. However, the method must set lower CourantFriedrichsLevy (CFL) number to obtain highorder numerical precision. In order to further explore efficient methods for optimum electromagnetic simulation, new improved time strategies referred as the highorder RungeKutta (RK) 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 (ADIFDTD) 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 timedependent Maxwell equations with unconditionally stable numerical schemes was proposed and developed [11]. The basic idea of the methods is to employ a LieTrotterSuzuki product formula to approximate the time evolution operator. As shown in [12], most of these methods can be seen as special cases of the timeoperatorsplitting 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 operatorsplitting 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 operatorsplitting methods for numerical solution of the timedependent Maxwell’s equations. In Section 2, the timedomain 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 highorder staggered finite difference is introduced for discretizing the threedimension curl operator in the spatial domain. The explicit discretized formulas are presented in Section 3. New highorder 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 threedimensional 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 timereversibility 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 BakerCampbellHausdorff (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 fourthorder 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 fifthorder 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 c_{1} = 0.1786, c_{2} = −0.0066, d_{1} = 0.7123, consists global minimum truncation error 0.00093. Now, we transfer to the discretization of spatial domain.
2.3. The MthOrder Difference Approximation for the FirstOrder Spatial Partial Derivative Operators
Let approximates the exact solution at the point in the nth time step. The following staggered thorder space difference operators are used to approximate the firstorder spatial partial derivative in the direction, that is, in threedimensional 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 thorder accuracy for the firstorder 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 curvelabeled 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 highorder 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.
(a)
(b)
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 1D zdirected, xpolarized TEM wave. The equations can be rewritten as The field components in the th timestep are denoted as here, is the spatial frequency along the zdirection. Substituted (15) into (9), we obtain Then (14) can be rewritten as with
Applying the scheme to (17), the timemarching 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 nondissipative.
In the threedimensional (3D) case, the continuoustime discretespace 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 CourantFriedrichsLevy 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.

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 tenpoint 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.

(a)
(b)
4.2. Scattering of the Dielectric Sphere
Next, consider a dielectric sphere illuminated by a plane wave propagating in the zdirection and E polarized in the xdirection. 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 tengrid 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 fourthorder accuracy in time domain and fourthorder, secondorder accuracy in spatial domain, respectively.
(a)
(b)
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.
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 largedomain and longtime simulation.
5. Conclusions
We present optimized operatorsplitting methods for numerical solution of the timedependent Maxwell equations in the time domain. The general highorder staggered finite difference is introduced for approximating the threedimensional 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
 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. View at: Google Scholar
 A. Taflove and S. C. Hagness, Computational Electrodynamics: The FiniteDifference TimeDomain Method, Artech House, London, UK, 3rd edition, 2005.
 J. Fang, Time domain finite difference computation for Maxwell's equations, Ph.D. thesis, University of California, Berkeley, Calif, USA, 1989.
 A. Yefet and P. G. Petropoulos, “A staggered fourthorder accurate explicit finite difference scheme for the timedomain Maxwell's equations,” Journal of Computational Physics, vol. 168, no. 2, pp. 286–315, 2001. View at: Publisher Site  Google Scholar
 S. V. Georgakopoulos, C. R. Birtcher, C. A. Balanis, and R. A. Renaut, “Higherorder finitedifference schemes for electromagnetic radiation, scattering, and penetration, part I: theory,” IEEE Antennas and Propagation Magazine, vol. 44, no. 1, pp. 134–142, 2002. View at: Publisher Site  Google Scholar
 S. V. Georgakopoulos, C. R. Birtcher, C. A. Balanis, and R. A. Renaut, “Higherorder finitedifference schemes for electromagnetic radiation, scattering, and penetration, part 2: applications,” IEEE Antennas and Propagation Magazine, vol. 44, no. 2, pp. 92–101, 2002. View at: Publisher Site  Google Scholar
 J. L. Young, D. Gaitonde, and J. S. Shang, “Toward the construction of a fourthorder difference scheme for transient EM wave simulation: staggered grid approach,” IEEE Transactions on Antennas and Propagation, vol. 45, no. 11, pp. 1573–1580, 1997. View at: Google Scholar
 J. S. Shang, “Highorder compactdifference schemes for timedependent Maxwell equations,” Journal of Computational Physics, vol. 153, no. 2, pp. 312–333, 1999. View at: Google Scholar
 T. Namiki, “A new FDTD algorithm based on alternatingdirection implicit method,” IEEE Transactions on Microwave Theory and Techniques, vol. 47, no. 10, pp. 2003–2007, 1999. View at: Google Scholar
 F. Zheng and Z. Chen, “Toward the development of a threedimensional unconditionally stable finitedifference timedomain method,” IEEE Transactions on Microwave Theory and Techniques, vol. 48, no. 9, pp. 1550–1558, 2000. View at: Publisher Site  Google Scholar
 J. S. Kole, M. T. Figge, and H. De Raedt, “Unconditionally stable algorithms to solve the timedependent Maxwell equations,” Physical Review E, vol. 64, no. 6, Article ID 066705, pp. 1–14, 2001. View at: Google Scholar
 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. 34, pp. 522–541, 2009. View at: Publisher Site  Google Scholar
 I. P. Omelyan, I. M. Mryglod, and R. Folk, “Optimized ForestRuth and Suzukilike algorithms for integration of motion in manybody systems,” Computer Physics Communications, vol. 146, no. 2, pp. 188–202, 2002. View at: Publisher Site  Google Scholar
 T. Hirono, W. Lui, S. Seki, and Y. Yoshikuni, “A threedimensional fourthorder finitedifference timedomain scheme using a symplectic integrator propagator,” IEEE Transactions on Microwave Theory and Techniques, vol. 49, no. 9, pp. 1640–1648, 2001. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar
 W. C. Chew, Waves and Fields in Inhomogenous Media, Van Nostrand Reinhold, New York, NY, USA, 1990.
Copyright
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.