Optimized Operator-Splitting 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 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 [4–6]. 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𝜕𝐀=−𝜇₋₁𝜕𝑡𝐇𝐄=(𝐀+𝐁)𝐇𝐄,𝟎𝜎^∗_3×3𝐈₃−𝜇⁻¹𝐑𝟎_3×3,𝐁=𝟎^3×3𝟎_3×3𝜀⁻¹𝐑−𝜀⁻¹𝜎𝐈₃,(1) where 𝜇 and 𝜀 are the permeability and permittivity, 𝜎 and 𝜎^∗ are electric and magnetic conductivities, 𝟎_3×3 and 𝐈₃ 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 Ψ⁻¹(𝑡)=Ψ(−𝑡) 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 𝑂(Δ𝑡⁵) 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Δ𝑡}+Υ₁Δ𝑡³+Υ₂Δ𝑡⁵+𝑂Δ𝑡⁷.(5) Note also, the propagators can be calculated analytically as follows:𝑒^{𝐀𝑐𝑙Δ𝑡}=⎛⎜⎜⎝exp−Δ𝑡𝑐^𝑙𝜎^∗𝜇𝐈³−1−exp−Δ𝑡𝑐^𝑙𝜎^∗/𝜇𝜎^∗𝟎_3×3𝐑𝐈₃𝑒_{𝐁𝑑𝑙Δ𝑡}⎞⎟⎟⎠,=⎛⎜⎜⎜⎝𝐈³𝟎_3×31−exp−Δ𝑡𝑑_𝑙𝜎/𝜀𝜎𝐑exp−Δ𝑡𝑑^𝑙𝜎𝜀𝐈³⎞⎟⎟⎟⎠.(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, {𝑐₁,𝑐₂,𝑑₁} to be determined. Using the Baker-Campbell-Hausdorff (BCH) formula, the explicit expression for Υ₁ and Υ₂ can be expressed asΥΥ₁=𝑓₁[𝐀,[𝐀,𝐁]]+𝑓₂[𝐀,[𝐀,𝐁]]₂,=+𝑓₃𝑓[₅𝐀,[𝐀,[𝐀,[𝐀,𝐁]]]]+𝑓₄[𝐀,[𝐀,[𝐁,[𝐀,𝐁]]]]+𝑓₇[𝐁,[𝐀,[𝐀,[𝐀,𝐁]]]]+𝑓₆[𝐁,[𝐁,[𝐁,[𝐀,𝐁]]]][𝐁,[𝐁,[𝐀,[𝐀,𝐁]]]]+𝑓₈[𝐀,[𝐁,[𝐁,[𝐀,𝐁]]]].(7) Here, [𝐀,𝐁]=𝐀𝐁−𝐁𝐀 and 𝑓_𝑖=𝑓_𝑖(𝑐₁,𝑐₂,𝑑₁),𝑖∈{1,2,…,8}. The detailed expression of 𝑓_𝑖 can be founded in [13]. The formula (5) represents a fourth-order scheme at Υ₁=0, 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 equationsmin⎧⎪⎨⎪⎩⎷⁸_𝑖=3subjectto:𝑓₁𝑓^2𝑖𝑐₁,𝑐₂,𝑑₁⎫⎪⎬⎪⎭𝑓𝑐₁,𝑐₂,𝑑₁=0₂𝑐₁,𝑐₂,𝑑₁=0.(8) One possible set of the coefficients are c₁ = 0.1786, c₂ = −0.0066, d₁ = 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+𝑠−1)!(𝑀/2−𝑠)!(2𝑠−1)²  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:𝐻^{𝑥𝑛+𝑙/5}𝑖,𝑗+1=exp−𝑤₁2,𝑘+12+1−exp−×𝑤𝐻^{𝑥𝑛+(𝑙−1)/5}₁𝑖,𝑗+12,𝑘+12𝑤₁⋅ℜ_𝑘⋅𝐸^{𝑦𝑛+(𝑙−1)/5}𝑖,𝑗+1−ℜ_𝑗2,𝑘+12⋅𝐸^{𝑧𝑛+(𝑙−1)/5}𝑖,𝑗+1𝐸2,𝑘+12_{𝑥𝑛+𝑙/5},=exp𝑖−+𝑤12₂,𝑗,𝑘+1−exp−×𝑤𝐸^{𝑥𝑛+(𝑙−1)/5}₂𝑖+12,𝑗,𝑘𝑤₂⋅ℜ_𝑗⋅𝐻^{𝑧𝑛+𝑙/5}𝑖+1−ℜ_𝑘2,𝑗,𝑘⋅𝐻^{𝑦𝑛+𝑙/5}𝑖+12,𝑗,𝑘(10) with 𝑤₁=𝑐^𝑙Δ𝑡𝜎^∗(𝑖,𝑗+1/2,𝑘+1/2)𝑤𝜇(𝑖,𝑗+1/2,𝑘+1/2)₂,=𝑑^𝑙Δ𝑡𝜎(𝑖+1/2,𝑗,𝑘)𝜀(𝑖+1/2,𝑗,𝑘).(11)

Here (i)we use Padé(0,3) and Padé(1,2) to approximate the expressions of exp(−𝑤_𝑖) and (1−exp(−𝑤_𝑖))/𝑤_𝑖:exp−𝑤_𝑖≈11−exp−𝑤1_𝑖+𝑤_𝑖+𝑤^2𝑖/2+𝑤^3𝑖/6,𝑤_𝑖≈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−𝛼𝑠₂_𝑠2𝜀₂𝑑𝑠₂,(13) 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 𝛼=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.

(a)

(b)

Figure 1 

Comparison between two kinds of approximation. (a) For 𝑒−𝑤 and (b) (1−𝑒−𝑤)/𝑤.

Figure 2 

Subcell technique specifies for the Ex field calculation in (4,𝑀) 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 (4,𝑀) scheme. For clarity, the discussion begins with a 1-D z-directed, x-polarized TEM wave. The equations can be rewritten as𝜕𝜕𝑡𝐻𝐸^𝑦_𝑥=⎛⎜⎜⎜⎝0−1−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𝜕𝜕𝑡𝐻𝐸^𝑦_𝑥=⎛⎜⎜⎜⎝0−1𝜇𝜂^𝑧−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𝑛+10𝑛+1}=𝑆𝐻𝐸^𝑛0𝑛0,(19) where𝑆=⁵_𝑙=1⎛⎜⎜⎝−𝜂10_𝑧𝑑_𝑙Δ_𝑡𝜀1⎞⎟⎟⎠⎛⎜⎜⎝1−𝜂^𝑧𝑐_𝑙Δ_𝑡01𝜇⎞⎟⎟⎠.(20) Solving for the eigenvalues, we can obtain 𝜆_1,2=tr(𝑆)±𝑗⁰√4−tr(𝑆)²2,(21) wheretr(𝑆)=2+⁵_𝑙=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𝑔_𝑙(𝜀𝜇)⁻¹Δ^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+₅_𝑙=1𝑔_𝑙−Δ𝑡²^𝑙|||||≤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.

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.

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.

Table 2 

Comparison of results for FDTD(2, 2) and (4, 4) scheme.

(a)

(b)

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 𝜀_𝑟=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 asError=||RCS−RCS^∗||||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.

(a)

(b)

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 Δ=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.

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

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