# A Compact Unconditionally Stable Method for  Time-Domain Maxwell's Equations

Authors: Zhuo Su; Yongqin Yang; Yunliang Long
Journal: International Journal of Antennas and Propagation
 (2013)

---

## Abstract

Higher order unconditionally stable methods are effective ways for simulating field behaviors of electromagnetic problems since they are free of Courant-Friedrich-Levy conditions. The development of accurate schemes with less computational expenditure is desirable. A compact fourth-order split-step unconditionally-stable finite-difference time-domain method (C4OSS-FDTD) is proposed in this paper. This method is based on a four-step splitting form in time which is constructed by symmetric operator and uniform splitting. The introduction of spatial compact operator can further improve its performance. Analyses of stability and numerical dispersion are carried out. Compared with noncompact counterpart, the proposed method has reduced computational expenditure while keeping the same level of accuracy. Comparisons with other compact unconditionally-stable methods are provided. Numerical dispersion and anisotropy errors are shown to be lower than those of previous compact unconditionally-stable methods.


---

## Body

## 1. Introduction

Over the past few decades, finite-difference time-domain (FDTD) method has been widely applied to a variety of electromagnetic problems. Nevertheless, the maximum time step size of conventional FDTD is constrained by the Courant-Friedrich-Levy (CFL) condition, thus degrading the efficiency of computations. Recently, unconditionally stable FDTD emerges to be a promising approach since the selection of time step size is no longer restricted by the CFL condition. Computation efficiency can be significantly improved. Various methods, such as alternating-direction-implicit FDTD (ADI-FDTD) and split-step FDTD (SS-FDTD), were presented in the literature [1, 2]. A novel four-step split-step method was proposed in [3]. Unlike conventional split-step FDTD method which is based on the exponential evolution operator, this method simultaneously utilizes symmetric operator and uniform splitting to decompose Maxwell’s matrix into four matrices, then the time step is divided into four substeps accordingly. It has new splitting form, and the execution procedure of this method is symmetric. This method has been proven to outperform the ADI-FDTD method, the SS-FDTD(1, 2), and the SS-FDTD(2, 2) methods in dispersion characteristic [3].Although the time step of unconditionally stable method is no longer bounded by stability, it is still restricted by numerical dispersion error which at worst can render the solution of electromagnetic problems impractical. One of the effective approaches for mitigating dispersion errors is using higher order spatial approximation. Some explicit higher order unconditionally stable methods based on the Taylor central finite-difference scheme were developed [4–6]. Moreover, compact differential operators have been utilized to construct implicit higher order unconditionally stable methods [7, 8]. Besides pursuing higher order spatial accuracy, the higher order techniques in time is also attracting research interest. An ADI-FDTD with fourth-order accuracy in time was developed in [9] to achieve lower phase velocity error for sufficiently fine mesh.The extension of the method proposed in [3] to three-dimensional domains was reported in [10]. Higher order central differential operator based on the Taylor central finite-difference was utilized to approximate the spatial derivatives. The generalization of the arbitrary-order formulation was also discussed. However, following its approach and considering two-dimensional simulations, it is found that the fourth-order four-step split-step method (4OSS-FDTD) involves inversion operation of diagonal band matrix with seven nonzero elements in a row, which may burden the computational complexity.In this paper, a compact fourth-order method is proposed and is denoted as C4OSS-FDTD. A compact scheme with fourth-order accuracy was utilized for spatial discretization. The proposed C4OSS-FDTD method is demonstrated to be unconditionally stable, and better computation efficiency is achieved by compressing the bandwidth of diagonal band matrix to be inversed from seven to five. It is capable of preserving the same level of accuracy as the 4OSS-FDTD method. In addition, compared with previous compact fourth-order unconditionally stable methods, numerical dispersion and anisotropy errors are reduced.This paper is organized as follows. Section2 first outlines the four-step split-step temporal scheme, and then the 4OSS-FDTD method with higher order spatial accuracy is derived by using the fourth-order Taylor central-difference spatial operator. The C4OSS-FDTD method is proposed in Section 3, and the formulation of the proposed method is given in detail. Section 4 proves the unconditional stability of the proposed C4OSS-FDTD method. Numerical dispersion of the C4OSS-FDTD method is analyzed in Section 5. Section 6 illustrates the numerical results, and Section 7 concludes the paper.


## 2. The 4OSS-FDTD Method

Let us first outline the four-step split-step scheme for time advancement [3], which is also adopted in the formulation of the C4OSS-FDTD method. Moreover, applying the fourth-order Taylor central finite-difference scheme to the spatial operator, the explicit 4OSS-FDTD method is derived, which will be involved in the comparison with C4OSS-FDTD later. The generalization of arbitrary-order four-step split-step method in three dimensional cases was discussed in [10]. This section will derive the formulation specifically for two dimensional case, which also can be considered as an enhanced version of the work presented in [3].Consider Maxwell’s equations of two dimensional TM mode in homogenous, isotropic, and sourceless medium with permittivityε and permeability μ. They can be represented in a matrix form as follows:
(1)∂U∂t=[M]U,
where U=[EzHxHy]T(2)[M]=[0-1ε∂∂y1ε∂∂x-1μ∂∂y001μ∂∂x00].Following the approach presented in [3], the matrix [M] can be decomposed into four parts by using symmetric operator and uniform splitting, as indicated in
(3)∂U∂t=[A]2U+[B]2U+[A]2U+[B]2U,
where
(4)[A]=[001ε∂∂x0001μ∂∂x00][B]=[0-1ε∂∂y0-1μ∂∂y00000].Using the Crank-Nicolson scheme, a single time step fromnΔt to (n+1) Δt  is uniformly divided into four substeps:
(5)(I-Δt4[A])Un+1/4=(I+Δt4[A])Un,(I-Δt4[B])Un+2/4=(I+Δt4[B])Un+1/4,(I-Δt4[A])Un+3/4=(I+Δt4[A])Un+2/4,(I-Δt4[B])Un+1=(I+Δt4[B])Un+3/4,
where I is the identity matrix of 3×3 and Δt is the time step size. The equations corresponding to the first substep from nΔt  to (n+1/4) Δt can be generalized as follows:
(6)Ez|I,Jn+1/4=Ez|I,Jn+b2(∂Hy∂x|I,Jn+1/4+∂Hy∂x|I,Jn),(7)Hx|I,J+1/2n+1/4=Hx|I,J+1/2n,(8)Hy|I+1/2,Jn+1/4=Hy|I+1/2,Jn+d2(∂Ez∂x|I+1/2,Jn+1/4+∂Ez∂x|I+1/2,Jn),
where b=Δt/(2ε),  d=Δt/(2μ).To reduce the numerical dispersion error, a fourth-order central finite-difference operator based on the Taylor central finite-difference scheme is utilized to approximate the spatial derivatives, as given in(9)∂f∂ξ|i=p(fi+1/2-fi-1/2Δξ)+q(fi+3/2-fi-3/2Δξ),
where p=9/8, q=-1/24. Equations (6), (8) and (9) yield the following system where Ez is updated from time step n to n+1/4:
(10)R=S+T,
where
(11)R=-bdq24(Δx)2Ez|I+3,Jn+1/4+-bdpq2(Δx)2Ez|I+2,Jn+1/4+bdp(2q-p)4(Δx)2Ez|I+1,Jn+1/4+[1+bd(p2+q2)2(Δx)2]Ez|I,Jn+1/4+bdp(2q-p)4(Δx)2Ez|I-1,Jn+1/4+-bdpq2(Δx)2Ez|I-2,Jn+1/4+-bdq24(Δx)2Ez|I-3,Jn+1/4,S=bdq24(Δx)2Ez|I+3,Jn+bdpq2(Δx)2Ez|I+2,Jn+bdp(p-2q)4(Δx)2Ez|I+1,Jn+[1-bd(p2+q2)2(Δx)2]Ez|I,Jn+bdp(p-2q)4(Δx)2Ez|I-1,Jn+bdpq2(Δx)2Ez|I-2,Jn+bdq24(Δx)2Ez|I-3,Jn,T=bΔx[p(Hy|I+1/2,jn-Hy|I-1/2,Jn)+q(Hy|I+3/2,Jn-Hy|I-3/2,Jn)].It is found that the solution of (10) involves inversion of a diagonal band matrix with seven nonzero elements in a row. The magnetic components are subsequently updated from (7) and (8) that can be solved straightforwardly to advance one substep. The equations corresponding to the remaining substeps are presented in (12). The remaining substeps can be completed in a similar manner as follows:
(12)Ez|I,Jn+2/4=Ez|I,Jn+1/4-b2(∂Hx∂y|I,Jn+2/4+∂Hx∂y|I,Jn+1/4),Hx|I,J+1/2n+2/4=Hx|I,J+1/2n+1/4-d2(∂Ez∂y|I,J+1/2n+2/4+∂Ez∂y|I,J+1/2n+1/4),Hy|I+1/2,Jn+2/4=Hy|I+1/2,Jn+1/4,Ez|I,Jn+3/4=Ez|I,Jn+2/4+b2(∂Hy∂x|I,Jn+3/4+∂Hy∂x|I,Jn+2/4),Hx|I,J+1/2n+3/4=Hx|I,J+1/2n+2/4,Hy|I+1/2,Jn+3/4=Hy|I+1/2,Jn+2/4+d2(∂Ez∂x|I+1/2,Jn+3/4+∂Ez∂x|I+1/2,Jn+2/4),Ez|I,Jn+1=Ez|I,Jn+3/4-b2(∂Hx∂y|I,Jn+1+∂Hx∂y|I,Jn+3/4),Hx|I,J+1/2n+1=Hx|I,J+1/2n+3/4-d2(∂Ez∂y|I,J+1/2n+1+∂Ez∂y|I,J+1/2n+3/4),Hy|I+1/2,Jn+1=Hy|I+1/2,Jn+3/4.


## 3.  Formulations of the C4OSS-FDTD Method

The C4OSS-FDTD method is proposed based on the previous four-step split-step approach for time advancement. However, the spatial derivatives are approximated by different scheme. The benefits of this approach will be explained at the end of this section by comparing with the noncompact 4OSS-FDTD method.For the C4OSS-FDTD method, a fourth-order accurate compact finite-difference scheme will be used as spatial differential operator, as given in(13)α∂f∂ξ|I-1+∂f∂ξ|I+α∂f∂ξ|I+1=βfI+1/2-fI-1/2Δξ,
where α=1/22, β=12/11, and f is one field component of either Ez, Hx, or Hy.  Δξ(ξ=x,y) denotes the discretization interval in either x or y direction.For simplicity and clarity, the following notations are used to represent various operations on the field components:(14)Qx(fI,J)=αfI-1,J+fI,J+αfI+1,J,Qy(fI,J)=αfI,J-1+fI,J+αfI,J+1,Rx(fI,J)=βfI+1/2,J-fI-1/2,JΔx,Ry(fI,J)=βfI+1/2,J-fI-1/2,JΔy.According to (13), the following relation holds
(15)Qξ(∂f∂ξ|I,J)=Rξ(fI,J).ApplyingQx operation twice to (6) on both sides, combining (15), and using (8) to eliminate Hy at time step n+1/4, the following matrix equation is derived after some algebraic manipulations:
(16)(Qx2-bdRx24)(Ez|I,Jn+1/4)=(Qx2+bdRx24)(Ez|I,Jn)+bQxRx(Hy|I,Jn),
where Qx2 represents applying the operator Qx twice and Rx2 represents applying the operator Rx twice. QxRx indicates applying Rx operator first, and then the Qx operation is following. Expanding (16) yields a system of(17)R′=S′+T′,
where
(18)R′=α2Ez|I+2,Jn+1/4+(2α-bdβ24Δx2)Ez|I+1,Jn+1/4+(1+2α2+bdβ22Δx2)Ez|I,Jn+1/4+(2α-bdβ24Δx2)Ez|I-1,Jn+1/4+α2Ez|I-2,Jn+1/4,S′=α2Ez|I+2,Jn+(2α+bdβ24Δx2)Ez|I+1,Jn+(1+2α2-bdβ22Δx2)Ez|I,Jn+(2α+bdβ24Δx2)Ez|I-1,Jn+α2Ez|I-2,Jn,T′=bβΔx(αHy|I+3/2,Jn+(1-α)Hy|I+1/2,Jn+(α-1)Hy|I-1/2,Jn-αHy|I-3/2,Jn).E
z can be updated from (17), the solution of (17) involves inversion of diagonal band matrix with five nonzero elements in a row. However, as analyzed in the previous section, the 4OSS-FDTD method involves a diagonal band matrix with seven nonzero elements in a row during the same procedure. The discrepancy of bandwidth of the matrix will affect the computational performance, especially for solving large-scale problems.By similar approach, the other matrix equations for the remaining three substeps are obtained as follows:(1)
substep 2, fromn+1/4 to n+2/4(19)(Qy2-bdRy24)(Ez|I,Jn+2/4)=(Qy2+bdRy24)(Ez|I,Jn+1/4)-bQyRy(Hx|I,Jn+1/4),(2)
substep 3, fromn+2/4 to n+3/4(20)(Qx2-bdRx24)(Ez|I,Jn+3/4)=(Qx2+bdRx24)(Ez|I,Jn+2/4)+bQxRx(Hy|I,Jn+2/4),(3)
substep 4, fromn+3/4 to n+1(21)(Qy2-bdRy24)(Ez|I,Jn+1)=(Qy2+bdRy24)(Ez|I,Jn+3/4)-bQyRy(Hx|I,Jn+3/4).At each substep, both 4OSS-FDTD and C4OSS-FDTD methods are required to solve a diagonal band matrix. Because the inversion of aL band N×N diagonal band matrix yields complexity of O(LN2), the bandwidth of the diagonal band matrix is a crucial factor for computation efficiency. Consequently the proposed C4OSS-FDTD method is advantageous with a computation efficiency gain.


## 4. Stability Analysis ofthe C4OSS-FDTD Method

To analyze the stability of C4OSS-FDTD method, the conventional Fourier method is used.kξ is defined as spatial frequency along the ξ direction. The field components in spectral domain at the nth time step can be written as
(22)Ez|I,Jn=Ezne-j(kxIΔx+kyJΔy),Hx|I,Jn=Hxne-j(kxIΔx+kyJΔy),Hy|I,Jn=Hyne-j(kxIΔx+kyJΔy).Defining a vectorFn=[EznHxnHyn]T to represent the amplitude of field components at the nth time step, the field components can be written as
(23)F|I,Jn=Fne-j(kxIΔx+kyJΔx).Substituting (23) into (13), the spatial derivative operator is represented as
(24)∂∂ξ=jγξ=j-2(β/Δξ)sin(kξΔx/2)2αcos(kξΔx)+1.The time marching is represented in terms of a growth matrix as follows:(25)Fn+1=[Λ]Fn=[ΛB][ΛA][ΛB][ΛA]Fn,
where
(26)[ΛA]=[CxDx04jbγxDx0104jdγxDx0CxDx],[ΛB]=[CyDy-4jbγyDy0-4jdγyDyCyDy0001],Cξ=4-bdγξ2Dξ=4+bdγξ2.The eigenvalues of growth matrixΛ can be obtained as follows:
(27)λ1=1,λ2=X+jY2-X2Y,λ3=X-jY2-X2Y,
where
(28)X=256-384bd(γx2+γy2)+16b2d2×(γx4+γy4-4γx2γy2)+b3d3γx2γy2(γx2+γy2)+b4d4γx4γy4,Y=256+128bd(γx2+γy2)+16b2d2×(γx4+γy4+4γx2γy2)+8b3d3γx2γy2(γx2+γy2)+b4d4γx4γy4.Note thatγξ is derived from a compact finite-difference operator and is a real number for α=1/22, β=12/11. Thus the magnitude of each eigenvalue of the growth matrix Λ equals to unity, |λ1|=|λ2|=|λ3|=1. The proposed C4OSS-FDTD method is proven to be unconditionally stable.Whenα=0, β=1, (13) reduces to the second-order Taylor central finite-difference operator. In such case, (24) reduces to the form that is coincident with the analysis of [3].


## 5. Numerical Dispersion ofthe C4OSS-FDTD Method

An expression of monochromatic wave with frequencyω is considered as follows:
(29)Fn+1=ejωΔtFn.Substitution of (29) into (25) gives rise to:
(30)(ejωΔtI-Λ)Fn=0.To obtain nontrivial solution of (30), the determinant of the coefficient matrix should be 0 which leads to
(31)det{ejωΔtI-Λ}=0.Specifically, in two dimensional simulation with uniform cells ofΔx=Δy=λ/Nλ, the dispersion relation can be written in a simplified form as given in (32). Where Nλ is the grid sampling density per wavelength, c0 is the speed of light, and vp is the numerical velocity. S=(c0Δt)/Δx is the Courant number which is defined in [11]. θ is the propagation direction of the numerical wave with respect to the x axis. One has
(32)tan2(πSNλ)=PQ,
where
(33)P=256S2(ζx2+ζy2)+64S4ζx2ζy2,Q=256-128S2(ζx2+ζy2)+16S4(ζx4+ζy4)+8S6ζx2ζy2(ζx2+ζy2)+S8ζx4ζy4,ζx=βsin((πcosθ/Nλ)(c0/vp))2αcos((2πcosθ/Nλ)(c0/vp))+1,ζy=βsin((πsinθ/Nλ)(c0/vp))2αcos((2πsinθ/Nλ)(c0/vp))+1.


## 6.  Numerical Results and Discussion

Numerical results are given and discussed in this section. The proposed C4OSS-FDTD, 4OSS-FDTD, compact ADI-FDTD [7], and CSS-FDTD [8] methods are included in the comparison. Note that all the methods involved have fourth-order accuracy in space. Normalized phase velocity vp/c0 is shown, as propagation direction varies. Uniform cells are used in the computation. Nλ is selected to be 20, and the CFL limit for Δt  is denoted by Δtc=Δx/2c0. Time step size is selected as Δt=CFLN×Δtc. From Figure 1, it can be observed that the numerical dispersion of the C4OSS-FDTD method is almost identical to that of the 4OSS-FDTD method at different CFLN values. Therefore, the C4OSS-FDTD method preserves the same level of accuracy in comparison with the 4OSS-FDTD method while benefiting from the improvement of computation efficiency as analyzed before. Figure 2 presents the normalized phase velocity of the C4OSS-FDTD, compact ADI-FDTD, and CSS-FDTD methods for CFLN = 1, 3. As can be seen from Figure 2, the dispersion error of the C4OSS-FDTD method is lower than that of the compact ADI-FDTD and CSS-FDTD methods.Figure 1
Normalized phase velocity of C4OSS-FDTD and 4OSS-FDTD (Nλ=20, CFLN = 1, 2, 4).Figure 2
Normalized phase velocity of fourth-order compact methods (Nλ=20, CFLN = 1, 3).The variation of anisotropy with respect toNλ is illustrated in Figure 3. Here, normalized velocity-anisotropy error is defined as
(34)ERRaniso=20·lg(max[vp(θ)]-min[vp(θ)]min[vp(θ)]).Figure 3
Normalized velocity-anisotropy error against sample density (CFLN = 2).For the C4OSS-FDTD and 4OSS-FDTD methods, their discrepancy of velocity-anisotropy error is negligibly small. Only at coarser meshes, tiny difference can be noticed where the C4OSS-FDTD method performs slightly better. It becomes apparent that the C4OSS-FDTD method outperforms the other two compact methods. Since the C4OSS-FDTD method has been proved to be accurate and less computationally expensive, it is potentially advantageous to electrically large simulations such as ground wave propagation. In that application, the whole system is sensitive to computational expenditure because the FDTD computation space is very large.From conclusions made in [7, 8], the compact ADI-FDTD and CSS-FDTD methods have almost the same dispersion characteristic as their noncompact counterparts for different time-step size. The numerical results in this section also indicate that the dispersion errors of the C4OSS-FDTD and 4OSS-FDTD methods are essentially at the same order. It is observed that the introduction of compact operators with fourth-order spatial accuracy in these unconditionally stable methods does not significantly improve their dispersion characteristics. The benefits of applying compact spatial operator are mainly on the reduction of computational expenditure. Further investigation of other compact unconditionally stable methods and their non-compact higher order counterparts needs to be performed to verify this observation in more general. However, noting that when compact spatial operator is applied to other time-marching schemes such as symplectic integrator, improvement in dispersion may be significant; [12] provides an example that the dispersion characteristic of compact symplectic FDTD method is apparently better than that of non-compact symplectic method.


## 7. Conclusion

This paper presented the C4OSS-FDTD method based on the combination of (1) four-step split-step approach for time advancement and (2) higher order compact difference scheme to the spatial differential operators. Compared with the 4OSS-FDTD method based on the Taylor central finite-difference scheme, the proposed method preserves the same level of accuracy while achieving higher computation efficiency. The stability and dispersion analysis is carried out. Moreover, better accuracy is demonstrated in comparison with previous compact unconditionally stable methods. It is also found that the introduction of higher order compact spatial operator in some unconditionally stable methods contributes only to reduced computational expenditure rather than improvement of the dispersion characteristic of their non-compact counterparts.



---

*Source: 689327-2013-03-27.xml*