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𝜕𝐇𝐄𝐇𝐄,𝜕𝑡=(𝐀+𝐁)𝐀=𝜇1𝜎𝐈3𝜇1𝐑𝟎3×3𝟎3×3𝟎,𝐁=3×3𝟎3×3𝜀1𝐑𝜀1𝜎𝐈3,(1) where 𝜇 and 𝜀 are the permeability and permittivity, 𝜎 and 𝜎 are electric and magnetic conductivities, 𝟎3×3 and 𝐈3 are 3×3 zeros matrix and identity matrix, 𝐑 is 3×3 matrix representing three-dimensional curl operator. Equation (1) can be cast in the following compact form𝑑𝐙𝑑𝑡=𝐋𝐙(𝑡).(2) 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𝐙(Δ𝑡)=𝑒𝐋Δ𝑡𝐙(0)𝑒(𝐀+𝐁)Δ𝑡𝐙(0),(3) 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𝑒Δ𝑡(𝐀+𝐁)=𝑚𝑙=1𝑒𝐁𝑑𝑙Δ𝑡𝑒𝐀𝑐𝑙Δ𝑡+𝑂Δ𝑡𝑛+1.(4) 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 Ψ1(𝑡)=Ψ(𝑡) of time evolution operator Ψ(𝑡)=𝑒𝐋𝑡) can also be reproduced by imposing additional time reversible constraints on the coefficients, namely, 𝑐𝑙=𝑐𝑚𝑙+1, and 𝑑𝑙=𝑑𝑚𝑙 with 𝑑𝑚=0. 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 𝑚1, of expensive force every time step. In this paper, particularly, we choose 𝑚=5 and minimize the truncation errors to 𝑂(Δ𝑡5) significantly with a little additional computational cost.

For 𝑚=5, the extended splitting result can be represented in the form𝑒Δ𝑡(𝐀+𝐁)=𝑒𝐀𝑐1Δ𝑡𝑒𝐁𝑑1Δ𝑡𝑒𝐀𝑐2Δ𝑡𝑒𝐁𝑑2Δ𝑡𝑒𝐀𝑐3Δ𝑡𝑒𝐁𝑑2Δ𝑡𝑒𝐀𝑐2Δ𝑡𝑒𝐁𝑑1Δ𝑡𝑒𝐀𝑐1Δ𝑡+Υ1Δ𝑡3+Υ2Δ𝑡5+𝑂Δ𝑡7.(5) Note also, the propagators can be calculated analytically as follows:𝑒𝐀𝑐𝑙Δ𝑡=expΔ𝑡𝑐𝑙𝜎𝜇𝐈31expΔ𝑡𝑐𝑙𝜎/𝜇𝜎𝐑𝟎3×3𝐈3,𝑒𝐁𝑑𝑙Δ𝑡=𝐈3𝟎3×31expΔ𝑡𝑑𝑙𝜎/𝜀𝜎𝐑expΔ𝑡𝑑𝑙𝜎𝜀𝐈3.(6) Here, the time reversible coefficients and the condition 5𝑙=1𝑐𝑙=5𝑙=1𝑑𝑙=1 have already been taken into account. With these assumptions, we actually only have three coefficients, namely, {𝑐1,𝑐2,𝑑1} to be determined. Using the Baker-Campbell-Hausdorff (BCH) formula, the explicit expression for Υ1 and Υ2 can be expressed asΥ1=𝑓1[[𝐀,𝐀,𝐁]]+𝑓2[[,Υ𝐀,𝐀,𝐁]]2=𝑓3[[[[𝐀,𝐀,𝐀,𝐀,𝐁]]]]+𝑓4[[[[𝐀,𝐀,𝐁,𝐀,𝐁]]]]+𝑓5[[[[𝐁,𝐀,𝐀,𝐀,𝐁]]]]+𝑓6[[[[𝐁,𝐁,𝐁,𝐀,𝐁]]]]+𝑓7[[[[𝐁,𝐁,𝐀,𝐀,𝐁]]]]+𝑓8[[[[.𝐀,𝐁,𝐁,𝐀,𝐁]]]](7) Here, [𝐀,𝐁]=𝐀𝐁𝐁𝐀 and 𝑓𝑖=𝑓𝑖(𝑐1,𝑐2,𝑑1),𝑖{1,2,,8}. The detailed expression of 𝑓𝑖 can be founded in [13]. The formula (5) represents a fourth-order scheme at Υ1=0, that is, 𝑓1 and 𝑓2 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 Υ2 to a minimum, and the coefficients {𝑐1,𝑐2,𝑑1} can be obtained by solving the system of equationsmin8𝑖=3𝑓2𝑖𝑐1,𝑐2,𝑑1𝑓subjectto:1𝑐1,𝑐2,𝑑1𝑓=02𝑐1,𝑐2,𝑑1=0.(8) 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𝑖𝑓𝑛𝑖,𝑗,𝑘𝜕𝑓|||𝜕𝑥𝑛𝑖,𝑗,𝑘=1Δ𝑥𝑀/2𝑠=1𝜉𝑠𝑓||𝑛𝑖+(2𝑠1)/2,𝑗,𝑘||𝑓𝑛𝑖(2𝑠1)/2,𝑗,𝑘+𝑂(Δ𝑥)𝑀.(9) Here, parameters 𝜉𝑠=(1)𝑠+1((𝑀1)!!)2/2𝑀2(𝑀/2+𝑠1)!(𝑀/2𝑠)!(2𝑠1)2  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 (4,𝑀) scheme is acquired.

For example, using Yee grid [1], the detailed expressions of 𝐻𝑥 and 𝐸𝑥 components in the (4,𝑀) scheme at the 𝑙th stage calculation after the 𝑛th time step are as follows:𝐻𝑥𝑛+𝑙/51𝑖,𝑗+21,𝑘+2=exp𝑤1×𝐻𝑥𝑛+(𝑙1)/51𝑖,𝑗+21,𝑘+2+1exp𝑤1𝑤1𝑘𝐸𝑦𝑛+(𝑙1)/51𝑖,𝑗+21,𝑘+2𝑗𝐸𝑧𝑛+(𝑙1)/51𝑖,𝑗+21,𝑘+2,𝐸𝑥𝑛+𝑙/51𝑖+2,𝑗,𝑘=exp𝑤2×𝐸𝑥𝑛+(𝑙1)/51𝑖+2+,𝑗,𝑘1exp𝑤2𝑤2𝑗𝐻𝑧𝑛+𝑙/51𝑖+2,𝑗,𝑘𝑘𝐻𝑦𝑛+𝑙/51𝑖+2,𝑗,𝑘(10) with 𝑤1=𝑐𝑙Δ𝑡𝜎(𝑖,𝑗+1/2,𝑘+1/2),𝑤𝜇(𝑖,𝑗+1/2,𝑘+1/2)2=𝑑𝑙Δ𝑡𝜎(𝑖+1/2,𝑗,𝑘).𝜀(𝑖+1/2,𝑗,𝑘)(11)

Here (i)we use Padé(0,3) and Padé(1,2) to approximate the expressions of exp(𝑤𝑖) and (1exp(𝑤𝑖))/𝑤𝑖:exp𝑤𝑖11+𝑤𝑖+𝑤2𝑖/2+𝑤3𝑖,/61exp𝑤𝑖𝑤𝑖1+𝑤𝑖/21+𝑤𝑖+𝑤2𝑖./3(12) 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𝛼𝜀=𝑠1𝑠1𝜀1𝑑𝑠1+1𝛼𝑠2𝑠2𝜀2𝑑𝑠2,(13) where 𝑠1 and 𝑠2 are the surfaces enclosed by curve-labeled single arrow (𝜀=𝜀1) and double arrow (𝜀=𝜀2) as indicated in Figure 2, α is free parameter, and we use 𝛼=9/8 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.

3.2. The Stability of the Proposed (4, M) Scheme

The conventional Fourier mode method is used to analyze the stability of the proposed (4,𝑀) scheme. For clarity, the discussion begins with a 1-D z-directed, x-polarized TEM wave. The equations can be rewritten as𝜕𝐻𝜕𝑡𝑦𝐸𝑥=10𝜇𝜕1𝜕𝑧𝜀𝜕0𝐻𝜕𝑧𝑦𝐸𝑥.(14) The field components in the 𝑛th time-step are denoted as𝐻𝑦𝐸𝑥|||||𝑛𝑧=𝑘=𝐻𝑛0𝑒𝑗0𝑘𝑧𝑘Δ𝑧𝐸𝑛0𝑒𝑗0𝑘𝑧𝑘Δ𝑧,(15) here, 𝑘𝑧 is the spatial frequency along the z-direction. Substituted (15) into (9), we obtain𝜕𝐹|||𝜕𝑧𝑛𝑘=𝑀/2𝑠=1𝜉𝑠𝑒𝑗0(𝑠1/2)𝑘𝑧Δ𝑧𝑒𝑗0(𝑠1/2)𝑘𝑧Δ𝑧Δ𝑧𝐹,𝐹=𝐻𝑦or𝐸𝑥.(16) Then (14) can be rewritten as𝜕𝐻𝜕𝑡𝑦𝐸𝑥=10𝜇𝜂𝑧1𝜀𝜂𝑧0𝐻𝑦𝐸𝑥(17) with 𝜂𝑧=𝑀/2𝑠=1𝜉𝑠𝑒𝑗0(𝑠1/2)𝑘𝑧Δ𝑧𝑒𝑗0(𝑠1/2)𝑘𝑧Δ𝑧Δ𝑧.(18)

Applying the (4,𝑀) scheme to (17), the time-marching relation can be expressed as𝐻0𝑛+1𝐸0𝑛+1𝐻=𝑆𝑛0𝐸𝑛0,(19) where𝑆=5𝑙=1𝜂10𝑧𝑑𝑙Δ𝑡𝜀1𝜂1𝑧𝑐𝑙Δ𝑡𝜇01.(20) Solving for the eigenvalues, we can obtain 𝜆1,2=tr(𝑆)±𝑗04tr(𝑆)22,(21) wheretr(𝑆)=2+5𝑙=1𝑔𝑙(𝜀𝜇)1Δ2𝑡𝜂2𝑧𝑙,𝑔𝑙=1𝑖1𝑗1<𝑖2𝑗2<<𝑖𝑙𝑗𝑙5𝑐𝑖1𝑑𝑗1𝑐𝑖2𝑑𝑗2𝑐𝑖𝑙𝑑𝑗𝑙+1𝑖1<𝑗1𝑖2<𝑗2𝑖𝑙<𝑗𝑙5𝑑𝑖1𝑐𝑗1𝑑𝑖2𝑐𝑗2𝑑𝑖𝑙𝑐𝑗𝑙.(22) We may conclude that |𝜆1,2|=1 if |tr(𝑆)|2, 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𝜕𝐇𝐄=𝜂𝜕𝑡0𝑥𝑒𝑥+𝜂𝑦𝑒𝑦+𝜂𝑧𝑒𝑧𝜇×𝜂𝑥𝑒𝑥+𝜂𝑦𝑒𝑦+𝜂𝑧𝑒𝑧𝜀𝐇𝐄.×0(23) Considering (𝜂2𝑥+𝜂2𝑦+𝜂2𝑧)<0, (23) can be rewritten in tensor form as𝜕𝐇𝐄=0𝜕𝑡𝜂2𝑥+𝜂2𝑦+𝜂2𝑧𝜇𝐊𝜂2𝑥+𝜂2𝑦+𝜂2𝑧𝜀𝐇𝐄,𝐊0(24) 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 tr(𝑆)=2+𝑚𝑙=1𝑔𝑙(𝜀𝜇)1Δ2𝑡𝜂2𝑥+𝜂2𝑦+𝜂2𝑧𝑙.(25)   Generally speaking, the stability of the scheme is assessed by the Courant-Friedrichs-Levy limit, that is, CFLmax, which can be expressed asCFLmax=Δ𝑡max𝜆𝑆,(26) where Δ𝑡max 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|||||2+5𝑙=1𝑔𝑙Δ𝑡2𝑙|||||||||2,Δ𝑡Δ𝑡max,𝑔𝑙=1𝑖1𝑗1<𝑖2𝑗2<<𝑖𝑙𝑗𝑙5𝑐𝑖1𝑑𝑗1𝑐𝑖2𝑑𝑗2𝑐𝑖𝑙𝑑𝑗𝑙+1𝑖1<𝑗1𝑖2<𝑗2𝑖𝑙<𝑗𝑙5𝑑𝑖1𝑐𝑗1𝑑𝑖2𝑐𝑗2𝑑𝑖𝑙𝑐𝑗𝑙,(27) here, 𝜆𝑆 is the spatial stability factor, which is defined as 𝜆𝑆=2𝐷×𝑀/2𝑠=1|𝜉𝑠|, and 𝐷 is the dimensional number. The numerical results of CFLmax with (4,𝑀) 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.

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 𝑘=23 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 (12,22,12+1/2), 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.

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 𝜀𝑟=4.0, 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||Error=RCSRCS||||RCS||,(28) 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.

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 Δ=0.1, the error of (4, 2) scheme becomes unacceptable. Figure 6 shows the comparisons between (4, 4) and (4, 2) schemes at CFL=0.6. In this case, whatever we choose Δ=0.05or Δ=0.1, 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 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 (4,𝑀) 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.


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