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

Authors: Z. X. Huang; X. L. Wu; W. E. I. Sha; B. Wu
Journal: International Journal of Antennas and Propagation
 (2012)

---

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


---

## Body

## 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 Section2, 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)∂∂t(HE)=(A+B)(HE),A=(-μ-1σ*I3-μ-1R03×303×3),B=(03×303×3ε-1R-ε-1σI3),
where μ and ε are the permeability and permittivity, σ and σ* are electric and magnetic conductivities, 03×3 and I3 are 3×3 zeros matrix and identity matrix, R is 3×3 matrix representing three-dimensional curl operator. Equation (1) can be cast in the following compact form(2)dZdt=LZ(t).
Here, Z(t)=[H(t),E(t)]T 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(3)Z(Δt)=eLΔtZ(0)≡e(A+B)ΔtZ(0),
where Δt denotes the time step, and the operator L=A+B 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 eLΔt does not allow to be evaluated exactly at any given Δt. However, at small enough values of Δt, the total propagator can be split using the formula(4)eΔt(A+B)=∏l=1meBdlΔteAclΔt+O(Δtn+1).
The coefficients cl  and dl  in this formula should be chosen in such a way to provide the highest possible value for n  at a given integer number m. Here m is stage number needed in every integer time step, and n is the order of the approximation. The main advantage of the above splitting is that the time-reversibility  Ψ(-t)Z(t)=Z(0) of solutions (following from the property Ψ-1(t)=Ψ(-t) of time evolution operator Ψ(t)=eLt) can also be reproduced by imposing additional time reversible constraints on the coefficients, namely, cl=cm-l+1, and dl=dm-l  with dm=0. Note also that the splitting method is quite general to build numerical integrators of arbitrary orders  n with m stages. But we cannot choose the stage m  too large, because this results in a too large number, namely m-1, of expensive force every time step. In this paper, particularly, we choose m=5 and minimize the truncation errors to O(Δt5) significantly with a little additional computational cost.Form=5, the extended splitting result can be represented in the form(5)eΔt(A+B)=eAc1ΔteBd1ΔteAc2ΔteBd2ΔteAc3ΔteBd2ΔteAc2ΔteBd1ΔteAc1Δt+Υ1Δt3+Υ2Δt5+O(Δt7).
Note also, the propagators can be calculated analytically as follows:(6)eAclΔt=(exp⁡(-Δtclσ*μ)I3-1-exp⁡(-Δtclσ*/μ)σ*R03×3I3),eBdlΔt=(I303×31-exp⁡(-Δtdlσ/ε)σRexp⁡(-Δtdlσε)I3).
Here, the time reversible coefficients and the condition ∑l=15cl=∑l=15dl=1 have already been taken into account. With these assumptions, we actually only have three coefficients, namely, {c1,c2,d1} to be determined. Using the Baker-Campbell-Hausdorff (BCH) formula, the explicit expression for Υ1 and Υ2 can be expressed as(7)Υ1=f1[A,[A,B]]+f2[A,[A,B]],Υ2=f3[A,[A,[A,[A,B]]]]+f4[A,[A,[B,[A,B]]]]+f5[B,[A,[A,[A,B]]]]+f6[B,[B,[B,[A,B]]]]+f7[B,[B,[A,[A,B]]]]+f8[A,[B,[B,[A,B]]]].
Here, [A,B]=AB-BA and fi=fi(c1,c2,d1),i∈{1,2,…,8}. The detailed expression of fi can be founded in [13]. The formula (5) represents a fourth-order scheme at Υ1=0, that is, f1 and f2 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 {c1,c2,d1} can be obtained by solving the system of equations(8)min⁡{∑i=38fi2(c1,c2,d1)}subject to:f1(c1,c2,d1)=0f2(c1,c2,d1)=0.
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

Letf|i,j,kn=f(iΔx,jΔy,kΔz;nΔt) approximates the exact solution f at the point (iΔx,jΔy,kΔz) in the nth time step. The following staggered Mth-order space difference operators are used to approximate the first-order spatial partial derivative in the ζ-direction, that is,  ∂ζ,ζ∈{x,y,z}  in three-dimensional curl operator R. So(9)Ri⋅fi,j,kn≡∂∂xf|i,j,kn=1Δx∑s=1M/2ξs[f|i+(2s-1)/2,j,kn-f|i-(2s-1)/2,j,kn]+O(Δx)M.
Here, parameters ξs=(-1)s+1((M-1)!!)2/2M-2(M/2+s-1)!(M/2-s)!(2s-1)2  for minimum truncation error in (9). Similarly, ℜj·fi,j,kn and ℜk·fi,j,kn are constructed in a similar manner and used to approximate ∂yfi,j,kn and ∂zfi,j,kn, respectively.


## 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)∂∂t(HE)=(A+B)(HE),A=(-μ-1σ*I3-μ-1R03×303×3),B=(03×303×3ε-1R-ε-1σI3),
where μ and ε are the permeability and permittivity, σ and σ* are electric and magnetic conductivities, 03×3 and I3 are 3×3 zeros matrix and identity matrix, R is 3×3 matrix representing three-dimensional curl operator. Equation (1) can be cast in the following compact form(2)dZdt=LZ(t).
Here, Z(t)=[H(t),E(t)]T 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(3)Z(Δt)=eLΔtZ(0)≡e(A+B)ΔtZ(0),
where Δt denotes the time step, and the operator L=A+B 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 eLΔt does not allow to be evaluated exactly at any given Δt. However, at small enough values of Δt, the total propagator can be split using the formula(4)eΔt(A+B)=∏l=1meBdlΔteAclΔt+O(Δtn+1).
The coefficients cl  and dl  in this formula should be chosen in such a way to provide the highest possible value for n  at a given integer number m. Here m is stage number needed in every integer time step, and n is the order of the approximation. The main advantage of the above splitting is that the time-reversibility  Ψ(-t)Z(t)=Z(0) of solutions (following from the property Ψ-1(t)=Ψ(-t) of time evolution operator Ψ(t)=eLt) can also be reproduced by imposing additional time reversible constraints on the coefficients, namely, cl=cm-l+1, and dl=dm-l  with dm=0. Note also that the splitting method is quite general to build numerical integrators of arbitrary orders  n with m stages. But we cannot choose the stage m  too large, because this results in a too large number, namely m-1, of expensive force every time step. In this paper, particularly, we choose m=5 and minimize the truncation errors to O(Δt5) significantly with a little additional computational cost.Form=5, the extended splitting result can be represented in the form(5)eΔt(A+B)=eAc1ΔteBd1ΔteAc2ΔteBd2ΔteAc3ΔteBd2ΔteAc2ΔteBd1ΔteAc1Δt+Υ1Δt3+Υ2Δt5+O(Δt7).
Note also, the propagators can be calculated analytically as follows:(6)eAclΔt=(exp⁡(-Δtclσ*μ)I3-1-exp⁡(-Δtclσ*/μ)σ*R03×3I3),eBdlΔt=(I303×31-exp⁡(-Δtdlσ/ε)σRexp⁡(-Δtdlσε)I3).
Here, the time reversible coefficients and the condition ∑l=15cl=∑l=15dl=1 have already been taken into account. With these assumptions, we actually only have three coefficients, namely, {c1,c2,d1} to be determined. Using the Baker-Campbell-Hausdorff (BCH) formula, the explicit expression for Υ1 and Υ2 can be expressed as(7)Υ1=f1[A,[A,B]]+f2[A,[A,B]],Υ2=f3[A,[A,[A,[A,B]]]]+f4[A,[A,[B,[A,B]]]]+f5[B,[A,[A,[A,B]]]]+f6[B,[B,[B,[A,B]]]]+f7[B,[B,[A,[A,B]]]]+f8[A,[B,[B,[A,B]]]].
Here, [A,B]=AB-BA and fi=fi(c1,c2,d1),i∈{1,2,…,8}. The detailed expression of fi can be founded in [13]. The formula (5) represents a fourth-order scheme at Υ1=0, that is, f1 and f2 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 {c1,c2,d1} can be obtained by solving the system of equations(8)min⁡{∑i=38fi2(c1,c2,d1)}subject to:f1(c1,c2,d1)=0f2(c1,c2,d1)=0.
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

Letf|i,j,kn=f(iΔx,jΔy,kΔz;nΔt) approximates the exact solution f at the point (iΔx,jΔy,kΔz) in the nth time step. The following staggered Mth-order space difference operators are used to approximate the first-order spatial partial derivative in the ζ-direction, that is,  ∂ζ,ζ∈{x,y,z}  in three-dimensional curl operator R. So(9)Ri⋅fi,j,kn≡∂∂xf|i,j,kn=1Δx∑s=1M/2ξs[f|i+(2s-1)/2,j,kn-f|i-(2s-1)/2,j,kn]+O(Δx)M.
Here, parameters ξs=(-1)s+1((M-1)!!)2/2M-2(M/2+s-1)!(M/2-s)!(2s-1)2  for minimum truncation error in (9). Similarly, ℜj·fi,j,kn and ℜk·fi,j,kn are constructed in a similar manner and used to approximate ∂yfi,j,kn and ∂zfi,j,kn, respectively.


## 3. Practical Implementation

### 3.1. Explicit Discretization Formulas

When one uses coefficientscl and dl of order four and substitutes the space difference operators with Mth-order accuracy for the first-order partial difference operators in R, the (4,M) scheme is acquired.For example, using Yee grid [1], the detailed expressions of Hx and Ex components in the (4,M) scheme at the lth stage calculation after the nth time step are as follows:(10)Hxn+l/5(i,j+12,k+12)=exp⁡(-w1)×Hxn+(l-1)/5(i,j+12,k+12)+1-exp⁡(-w1)w1⋅{Rk⋅Eyn+(l-1)/5(i,j+12,k+12)-Rj⋅Ezn+(l-1)/5(i,j+12,k+12)},Exn+l/5(i+12,j,k)=exp⁡(-w2)×Exn+(l-1)/5(i+12,j,k)+1-exp⁡(-w2)w2⋅{Rj⋅Hzn+l/5(i+12,j,k)-Rk⋅Hyn+l/5(i+12,j,k)}
with
(11)w1=clΔtσ¯*(i,j+1/2,k+1/2)μ(i,j+1/2,k+1/2),w2=dlΔtσ¯(i+1/2,j,k)ε¯(i+1/2,j,k).Here(i)
we use Padé(0,3) and Padé(1,2) to approximate the expressions of exp(-wi) and (1-exp(-wi))/wi:(12)exp⁡(-wi)≈11+wi+wi2/2+wi3/6,1-exp⁡(-wi)wi≈1+wi/21+wi+wi2/3.
Our new approximate acquired super stability and efficiency than the Padé(2, 2) method in [14] when wi 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 S can be expressed as(13)ε¯=αs1∬s1ε1ds1+1-αs2∬s2ε2ds2,
where s1 and s2 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.Comparison between two kinds of approximation. (a) Fore-w and  (b) (1-e-w)/w.
(a)(b)Figure 2
Subcell technique specifies for the Ex field calculation in(4,M) 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,M) scheme. For clarity, the discussion begins with a 1-D z-directed, x-polarized TEM wave. The equations can be rewritten as(14)∂∂t(HyEx)=(0-1μ∂∂z-1ε∂∂z0)(HyEx).
The field components in the nth time-step are denoted as(15)(HyEx)|z=kn=(H0n⋅e-j0⋅kzkΔzE0n⋅e-j0⋅kzkΔz),
here, kz is the spatial frequency along the z-direction. Substituted (15) into (9), we obtain(16)∂F∂z|kn=∑s=1M/2ξse-j0(s-1/2)kzΔz-ej0(s-1/2)kzΔzΔz⋅F,F=HyorEx.
Then (14) can be rewritten as(17)∂∂t(HyEx)=(0-1μηz-1εηz0)(HyEx)
with
(18)ηz=∑s=1M/2ξse-j0(s-1/2)kzΔz-ej0(s-1/2)kzΔzΔz.Applying the(4,M) scheme to (17), the time-marching relation can be expressed as(19)(H0n+1E0n+1)=S(H0nE0n),
where(20)S=∏l=15(10-ηzdlΔtε1)(1-ηzclΔtμ01).
Solving for the eigenvalues, we can obtain (21)λ1,2=tr⁡(S)±j04-tr⁡(S)22,
where(22)tr⁡(S)=2+∑l=15gl((εμ)-1Δt2ηz2)l,gl=∑1≤i1≤j1<i2≤j2<⋯<il≤jl≤5ci1dj1ci2dj2⋯cildjl+∑1≤i1<j1≤i2<j2≤⋯≤il<jl≤5di1cj1di2cj2⋯dilcjl.
We may conclude that |λ1,2|=1 if |tr(S)|≤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(23)∂∂t(HE)=(0-ηxex+ηyey+ηzezμ×ηxex+ηyey+ηzezε×0)(HE).
Considering (ηx2+ηy2+ηz2)<0, (23) can be rewritten in tensor form as(24)∂∂t(HE)=(0--(ηx2+ηy2+ηz2)μK̅⋅-(ηx2+ηy2+ηz2)εK̅⋅0)(HE),
where K̅ 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 (25)tr⁡(S)=2+∑l=1mgl((εμ)-1Δt2(ηx2+ηy2+ηz2))l.
Generally speaking, the stability of the scheme is assessed by the Courant-Friedrichs-Levy limit, that is, CFLmax, which can be expressed as(26)CFLmax=Δtmax⁡λS,
where Δtmax is the temporal stability factor, and the numerical time evolution operator will not blow up for all Δt, which can be determined according to the following inequality(27)|2+∑l=15gl(-Δt2)l|≤2,|Δt|≤Δtmax⁡,gl=∑1≤i1≤j1<i2≤j2<⋯<il≤jl≤5ci1dj1ci2dj2⋯cildjl+∑1≤i1<j1≤i2<j2≤⋯≤il<jl≤5di1cj1di2cj2⋯dilcjl,
here, λS is the spatial stability factor, which is defined as  λS=2D×∑s=1M/2|ξs|,  and D is the dimensional number. The numerical results of CFLmax with (4,M) 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.
MethodCFLmaxAccuracyD = 2D = 3(4, 2)1.050.86O(Δt4)+O(Δς2)(4, 4)0.900.73O(Δt4)+O(Δς4)(4, 6)0.850.69O(Δt4)+O(Δς6)FDTD(2, 2)0.710.58O(Δt2)+O(Δς2)Figure 3
Numerical  dispersion curves for (4, 4) scheme and FDTD(2, 2) method.


## 3.1. Explicit Discretization Formulas

When one uses coefficientscl and dl of order four and substitutes the space difference operators with Mth-order accuracy for the first-order partial difference operators in R, the (4,M) scheme is acquired.For example, using Yee grid [1], the detailed expressions of Hx and Ex components in the (4,M) scheme at the lth stage calculation after the nth time step are as follows:(10)Hxn+l/5(i,j+12,k+12)=exp⁡(-w1)×Hxn+(l-1)/5(i,j+12,k+12)+1-exp⁡(-w1)w1⋅{Rk⋅Eyn+(l-1)/5(i,j+12,k+12)-Rj⋅Ezn+(l-1)/5(i,j+12,k+12)},Exn+l/5(i+12,j,k)=exp⁡(-w2)×Exn+(l-1)/5(i+12,j,k)+1-exp⁡(-w2)w2⋅{Rj⋅Hzn+l/5(i+12,j,k)-Rk⋅Hyn+l/5(i+12,j,k)}
with
(11)w1=clΔtσ¯*(i,j+1/2,k+1/2)μ(i,j+1/2,k+1/2),w2=dlΔtσ¯(i+1/2,j,k)ε¯(i+1/2,j,k).Here(i)
we use Padé(0,3) and Padé(1,2) to approximate the expressions of exp(-wi) and (1-exp(-wi))/wi:(12)exp⁡(-wi)≈11+wi+wi2/2+wi3/6,1-exp⁡(-wi)wi≈1+wi/21+wi+wi2/3.
Our new approximate acquired super stability and efficiency than the Padé(2, 2) method in [14] when wi 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 S can be expressed as(13)ε¯=αs1∬s1ε1ds1+1-αs2∬s2ε2ds2,
where s1 and s2 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.Comparison between two kinds of approximation. (a) Fore-w and  (b) (1-e-w)/w.
(a)(b)Figure 2
Subcell technique specifies for the Ex field calculation in(4,M) 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,M) scheme. For clarity, the discussion begins with a 1-D z-directed, x-polarized TEM wave. The equations can be rewritten as(14)∂∂t(HyEx)=(0-1μ∂∂z-1ε∂∂z0)(HyEx).
The field components in the nth time-step are denoted as(15)(HyEx)|z=kn=(H0n⋅e-j0⋅kzkΔzE0n⋅e-j0⋅kzkΔz),
here, kz is the spatial frequency along the z-direction. Substituted (15) into (9), we obtain(16)∂F∂z|kn=∑s=1M/2ξse-j0(s-1/2)kzΔz-ej0(s-1/2)kzΔzΔz⋅F,F=HyorEx.
Then (14) can be rewritten as(17)∂∂t(HyEx)=(0-1μηz-1εηz0)(HyEx)
with
(18)ηz=∑s=1M/2ξse-j0(s-1/2)kzΔz-ej0(s-1/2)kzΔzΔz.Applying the(4,M) scheme to (17), the time-marching relation can be expressed as(19)(H0n+1E0n+1)=S(H0nE0n),
where(20)S=∏l=15(10-ηzdlΔtε1)(1-ηzclΔtμ01).
Solving for the eigenvalues, we can obtain (21)λ1,2=tr⁡(S)±j04-tr⁡(S)22,
where(22)tr⁡(S)=2+∑l=15gl((εμ)-1Δt2ηz2)l,gl=∑1≤i1≤j1<i2≤j2<⋯<il≤jl≤5ci1dj1ci2dj2⋯cildjl+∑1≤i1<j1≤i2<j2≤⋯≤il<jl≤5di1cj1di2cj2⋯dilcjl.
We may conclude that |λ1,2|=1 if |tr(S)|≤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(23)∂∂t(HE)=(0-ηxex+ηyey+ηzezμ×ηxex+ηyey+ηzezε×0)(HE).
Considering (ηx2+ηy2+ηz2)<0, (23) can be rewritten in tensor form as(24)∂∂t(HE)=(0--(ηx2+ηy2+ηz2)μK̅⋅-(ηx2+ηy2+ηz2)εK̅⋅0)(HE),
where K̅ 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 (25)tr⁡(S)=2+∑l=1mgl((εμ)-1Δt2(ηx2+ηy2+ηz2))l.
Generally speaking, the stability of the scheme is assessed by the Courant-Friedrichs-Levy limit, that is, CFLmax, which can be expressed as(26)CFLmax=Δtmax⁡λS,
where Δtmax is the temporal stability factor, and the numerical time evolution operator will not blow up for all Δt, which can be determined according to the following inequality(27)|2+∑l=15gl(-Δt2)l|≤2,|Δt|≤Δtmax⁡,gl=∑1≤i1≤j1<i2≤j2<⋯<il≤jl≤5ci1dj1ci2dj2⋯cildjl+∑1≤i1<j1≤i2<j2≤⋯≤il<jl≤5di1cj1di2cj2⋯dilcjl,
here, λS is the spatial stability factor, which is defined as  λS=2D×∑s=1M/2|ξs|,  and D is the dimensional number. The numerical results of CFLmax with (4,M) 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.
MethodCFLmaxAccuracyD = 2D = 3(4, 2)1.050.86O(Δt4)+O(Δς2)(4, 4)0.900.73O(Δt4)+O(Δς4)(4, 6)0.850.69O(Δt4)+O(Δς6)FDTD(2, 2)0.710.58O(Δt2)+O(Δς2)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∂ζ,ζ∈{x,y,z} 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 M≤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 dipoleP 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 Ezfield emanating along the plane of k=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.
FDTD(2, 2)(4, 4)Physical time200 ns200 nsCFL0.50.6Time step83.33 ps200 psSpatial step0.05 m0.1 mNo. of steps24001000Total run time184 sec131 secAverage CPU time/step0.0766 sec0.1310 secError1.59200.5107Radiation from a dipole with a ten-point PML. (a) The distribution of theEz field. (b) The contour of the Ez amplitude field.
(a)(b)


### 4.2. Scattering of the Dielectric Sphere

Next, consider a dielectric sphere illuminated by a plane wave propagating in thez-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 εr=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(28)Error=|RCS-RCS*||RCS*|,
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.The relative error of computed RCS with different scheme at CFL = 0.5. (a)E-plane  (b) H-plane.
(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Δ=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 computedE-plane RCS with different scheme at CFL = 0.6.Figure 7
The relative error of computedE-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.


## 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 dipoleP 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 Ezfield emanating along the plane of k=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.
FDTD(2, 2)(4, 4)Physical time200 ns200 nsCFL0.50.6Time step83.33 ps200 psSpatial step0.05 m0.1 mNo. of steps24001000Total run time184 sec131 secAverage CPU time/step0.0766 sec0.1310 secError1.59200.5107Radiation from a dipole with a ten-point PML. (a) The distribution of theEz field. (b) The contour of the Ez amplitude field.
(a)(b)


## 4.2. Scattering of the Dielectric Sphere

Next, consider a dielectric sphere illuminated by a plane wave propagating in thez-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 εr=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(28)Error=|RCS-RCS*||RCS*|,
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.The relative error of computed RCS with different scheme at CFL = 0.5. (a)E-plane  (b) H-plane.
(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Δ=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 computedE-plane RCS with different scheme at CFL = 0.6.Figure 7
The relative error of computedE-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,M) 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.



---

*Source: 956431-2011-07-06.xml*