Operator splitting methods for Maxwell’s equations in dispersive media with orientational polarization

Journal of Computational and Applied Mathematics 263 (2014) 160–188

Contents lists available at ScienceDirect

Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam

Operator splitting methods for Maxwell’s equations in dispersive media with orientational polarization V.A. Bokil ∗ , O.A. Keefer, A.C.-Y. Leung Department of Mathematics, Oregon State University, Corvallis, OR 97331-4605, United States

article

info

Article history: Received 17 October 2012 Received in revised form 25 November 2013 MSC: 65N10 65N15 Keywords: Maxwell’s equations Debye media Operator splitting Yee scheme Crank–Nicolson

abstract We present two operator splitting schemes for the numerical simulation of Maxwell’s equations in dispersive media of Debye type that exhibit orientational polarization (the Maxwell–Debye model). The splitting schemes separate the mechanisms of wave propagation and polarization to create simpler sub-steps that are easier to implement. In addition, dimensional splitting is used to propagate waves in different axial directions. We present a sequential operator splitting scheme and its symmetrized version for the Maxwell–Debye system in two dimensions. The splitting schemes are discretized using implicit finite difference methods that lead to unconditionally stable schemes. We prove that the fully discretized sequential scheme is a first order time perturbation, and the symmetrized scheme is a second order time perturbation of the Crank–Nicolson scheme for discretizing the Maxwell–Debye model. Numerical examples are presented that illustrate our theoretical results. © 2013 Elsevier B.V. All rights reserved.

1. Introduction In this paper we consider the numerical simulation of electromagnetic wave propagation in dispersive dielectrics. By a dispersive material we mean a dielectric (insulator) in which planar electromagnetic waves travel with phase velocities that depend on the frequency of the propagating waves. Thus, dispersive materials embody physical dispersion, and incident waves in these media spread and change shape, even when the medium is homogeneous [1]. In particular, we focus on dielectrics whose permittivity is represented by relaxing models. The Debye model for orientational polarization [2] was introduced to model the macroscopic polarization in such relaxing dielectrics like water or biological tissue, for example. This macroscopic polarization models the averaged response of a material to the incident electromagnetic field. The Debye model has been shown to accurately capture the experimental data obtained for the frequency dependent permittivity of several relaxing dielectrics such as biological media, homogeneous dispersive rock and homogeneous sea ice (see [3] and references therein). The combination of the Debye model and Maxwell’s equations, the Maxwell–Debye model, will be considered in this paper for modeling electromagnetic wave propagation in relaxing dielectrics. The Maxwell–Debye model is obtained by appending to Maxwell’s equations a set of ordinary differential equations (ODEs) that model the dynamic evolution of the macroscopic polarization driven by the electric field [4,5]. This approach is known as the auxiliary differential equation (ADE) method [6], for modeling electromagnetic pulse propagation in dispersive media. As shown in [3], several very different wave speeds in well defined spatial and temporal regimes are present in the relaxing dielectric modeled by the (1D) Maxwell–Debye system. The pulse response in such relaxing dielectrics is mainly a diffusion wave travelling with the slowest speed (corresponding to the static or zero-frequency permittivity) of harmonic waves and the fastest speed (corresponding to the infinite frequency permittivity) in the problem is only present in a small

∗

Corresponding author. Tel.: +1 541 737 2609; fax: +1 541 737 0517. E-mail address: [email protected] (V.A. Bokil).

0377-0427/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cam.2013.12.008

V.A. Bokil et al. / Journal of Computational and Applied Mathematics 263 (2014) 160–188

161

boundary layer. Thus, the Maxwell–Debye system is a stiff system of PDEs because of the fast decaying transient in the boundary layer [3]. The Finite Difference Time Domain (FDTD) Yee scheme [7], is an explicit method that numerically solves Maxwell’s equations. The Yee scheme was initially introduced to model wave propagation in non-dispersive media. It was later extended to model wave propagation in dispersive dielectrics (see [6] and references therein). Yee schemes are constructed for the dispersive Maxwell–Debye system by discretizing Maxwell’s equations using a space–time staggering of the electric and magnetic fields as done in the classical Yee scheme for the numerical simulation of non-dispersive wave propagation. Averaging is used to discretize the lower order terms in the Maxwell–Debye system, and the auxiliary ODEs are time discretized using a Crank–Nicolson method which is second order accurate in time. With this approach, the fully discretized Yee method for the dispersive Maxwell–Debye system is fully explicit, second order accurate in space and time and conditionally stable. The Courant–Friedrich–Lewy (CFL) [8] condition for the Yee scheme applied to the Maxwell–Debye system is the same as the CFL condition of the classical Yee scheme in a non-dispersive dielectric which has the fastest speed of wave propagation present in the dispersive Maxwell–Debye system. For the stiff problem modeled by the Maxwell–Debye system the CFL constraint in addition to accuracy requirements leads to very small time steps in numerical computations of the Yee scheme [9]. The Crank–Nicolson (CN) method is an implicit numerical technique that has also been developed for the numerical simulation of the Maxwell–Debye model and demonstrated to be unconditionally stable in [10,11]. Schemes that are unconditionally stable are well suited for problems involving geometries needing different details of discretization such as narrow slots [12]. For such geometries non-uniform meshing techniques can be created by using locally small spatial increments which do not require extremely small time steps in an unconditionally stable scheme. The alternating direction implicit (ADI)–FDTD scheme [13] offers additional reductions in computational time over fully implicit methods like the CN scheme while preserving the property of unconditional stability. The ADI–FDTD scheme has recently been extended to incorporate Debye media [12,14–17]. In this paper we derive new operator splitting methods for the Maxwell–Debye model. Operator splitting schemes decompose operators into sub-operators isolating different characteristics of the physical problem. A thorough understanding of the physical processes being modeled is a necessary component to construct the proper sub-problems in the splitting schemes [18]. Operator splitting schemes for Maxwell’s equations in non-dispersive materials have been considered in [19,20]. We extend the approach used in [19] for non-dispersive materials to construct splitting methods for the Maxwell–Debye model. For the numerical approximation of the Maxwell–Debye system, we enlist the operator splitting schemes to separate the effects of orientational polarization (response of the dispersive material to the incident electromagnetic field) from the spatial derivatives which model wave propagation. One of the main advantages of the new operator splitting schemes introduced in this paper over the FDTD or Yee scheme is their unconditional stability. The lack of a CFL constraint is particularly important for a numerical approximation to the stiff Maxwell–Debye problem. Thus, like the Crank–Nicolson scheme, our new schemes allow one to choose time steps based on accuracy considerations only; however, unlike the Crank–Nicolson scheme, the separation of different physical processes and dimensional splitting leads to numerical methods that are computationally more efficient in comparison. The operator splitting schemes are simple to implement and involve a series of one dimensional problems which either propagate waves or adjust for the material response to the incident field. In Section 2 we present the Maxwell–Debye model and its reduction to two dimensions. A sequential and a symmetrized operator splitting scheme are presented in Section 3. We prove the unconditional stability of the schemes via an energy analysis in Section 4, and in Sections 5 and 6 we consider the accuracy of the schemes. Here, it is shown that the fully discretized sequential and symmetrized splitting schemes are first order, and second order, respectively, in time perturbations of the second order accurate Crank–Nicolson (CN) scheme. Finally, in Section 7 numerical experiments are presented that illustrate the analysis and compare the splitting schemes with the Yee scheme for the Maxwell–Debye model.

2. The Maxwell–Debye model for orientational polarization We consider Maxwell’s equations which govern the electric field E and the magnetic field H in a domain Ω ⊂ R3 from time 0 to T given as

∂B + ∇ × E = 0 in Ω × (0, T ), ∂t ∂D 1 − ∇ × B = 0 in Ω × (0, T ). ∂t µ0 The fields D, B are the electric and magnetic flux densities, respectively. On the boundary, ∂ Ω , the perfect electric conductor (PEC) condition n × E = 0 on ∂ Ω × (0, T ),

(2.1a) (2.1b)

(2.2)

is imposed to terminate the computational domain, where the vector n is the outward unit normal vector to ∂ Ω . Lastly, we add initial conditions to the system.

162

V.A. Bokil et al. / Journal of Computational and Applied Mathematics 263 (2014) 160–188

Within the dielectric medium we have constitutive relations that relate the flux densities D, B to the electric and magnetic fields, respectively, as D = ϵ0 ϵ∞ E + P,

(2.3a)

B = µ0 H,

(2.3b)

where the constants ϵ0 and µ0 are the permittivity and permeability of free space, and are connected to the speed of light √ c0 by c0 = 1/ ϵ0 µ0 . The quantity P is called the (macroscopic) relaxation polarization. The presence of instantaneous polarization in the dielectric is accounted for by the coefficient ϵ∞ , the infinite frequency relative permittivity, in the constitutive relation (2.3a). We neglect any additional magnetic effects and assume that the magnetic constitutive relation (2.3b) for free space is also valid in the Debye medium. To model the macroscopic relaxation polarization we use the Debye model, which is an ordinary differential equation governing the evolution of the relaxation polarization driven by the electric field. In the simplest case of a single-pole Debye medium this model [5] is given as

∂P + P = ϵ0 ϵ∞ (ϵq − 1)E in Ω × (0, T ). (2.4) ∂t In Eq. (2.4) the parameter ϵq := ϵϵs is the ratio of static or zero frequency relative permittivity, denoted as ϵs , to the infinite ∞ frequency relative permittivity ϵ∞ . The parameter τ is the relaxation time associated with the polarization mechanism [1]. In general τ , ϵ∞ , and ϵs can be functions of space, but we assume here that all parameters are constant, and ϵs > ϵ∞ or √ √ ϵq > 1 and τ > 0. Corresponding to the two permittivities we define two speeds, c∞ := c0 / ϵ∞ , and cs := c0 / ϵs . Note that the relative permittivities ϵ∞ and ϵs are unitless parameters; see for example [1,3]. As shown in [3], these speeds with c∞ > cs are valid in well defined and separate spatio-temporal regimes. τ

To construct the 2D TE Maxwell–Debye model we make the assumption that no fields exhibit variation in the z direction, i.e. all partial derivatives with respect to z are zero. The electric field and polarization have two components, E = (Ex , Ey )T , P = (Px , Py )T and the magnetic field has one component Hz = H. Combining (2.4) with the constitutive relations (2.3a) and (2.3b), and substituting in the Maxwell equations (2.1) we get the following three partial differential equations which we call the 2D Maxwell–Debye TE curl equations:

∂H 1 =− curl E, ∂t µ0

(2.5a)

∂E 1 (ϵq − 1) 1 = curl H − E+ P, ∂t ϵ0 ϵ∞ τ ϵ0 ϵ∞ τ ∂P ϵ0 ϵ∞ (ϵq − 1) 1 = E − P, ∂t τ τ

(2.5b) (2.5c) ∂U

where for a vector field, U = (Ux , Uy )T , the scalar curl operator is curl U = ∂ xy − ∂∂Uyx , and for a scalar field V the vector curl operator is curl V =



∂V , − ∂∂Vx ∂y

T

[21]. All the fields in (2.5) are functions of position x = (x, y)T and time t.

System (2.5) along with the PEC boundary conditions (2.2) and initial conditions E(x, 0) = E0 (x), P(x, 0) = P0 (x) and H (x, 0) = H0 (x) for x ∈ Ω ⊂ R2 is well-posed. We define the following two function spaces: H (curl, Ω ) = {u ∈ (L2 (Ω ))2 ; curl u ∈ L2 (Ω )},

(2.6)

H0 (curl, Ω ) = {u ∈ H (curl, Ω ), n × u = 0 on ∂ Ω }.

(2.7)

Let (·, ·) denote the L inner product and ∥ · ∥2 the corresponding norm. Multiplying (2.5a) by v ∈ L (Ω ), (2.5b) by 2

2

2

u ∈ H0 (curl, Ω ), and (2.5c) by w ∈ L2 (Ω ) , integrating over the domain Ω and applying Green’s formula for the curl operator



(curl H , u) = (H , curl u) ,

∀u ∈ H0 (curl, Ω ),

(2.8)

we obtain the weak formulation

 ∂H , v = − (curl E, v) , ∀v ∈ L2 (Ω ), ∂t       ∂E ϵ0 ϵ∞ (ϵq − 1) 1 ϵ0 ϵ∞ , u = (H , curl u) − E, u + P, u , ∀u ∈ H0 (curl, Ω ), ∂t τ τ       1 ∂P 1 1 ,w = E, w − P, w , ∀w ∈ (L2 (Ω ))2 . ϵ0 ϵ∞ (ϵq − 1) ∂ t τ ϵ0 ϵ∞ (ϵq − 1)τ 

µ0

(2.9a)

(2.9b)

(2.9c)

The following theorem shows the stability of the 2D Maxwell–Debye model (2.5) by showing that the model exhibits energy decay.

V.A. Bokil et al. / Journal of Computational and Applied Mathematics 263 (2014) 160–188

163

Theorem 2.1 (Energy Decay). Let Ω ⊂ R2 and suppose that E ∈ C (0, T ; H0 (curl, Ω )) ∩ C 1 (0, T ; (L2 (Ω ))2 ), P ∈ C 1 (0, T ;

2

L2 (Ω ) ), and H ∈ C 1 (0, T ; L2 (Ω )) are solutions of the weak formulation (2.9) for the Maxwell–Debye system of equations (2.5) along with PEC boundary conditions (2.2). Then the system exhibits energy decay,



E (t ) ≤ E (0) ∀t ≥ 0,

(2.10)

where the energy E (t ) is defined as

 2  12 √ 2 √ 2   1       E (t ) =  µ0 H (t ) +  ϵ0 ϵ∞ E(t ) +   P(t )  .  ϵ0 ϵ∞ (ϵq − 1)  2 2 2 

Proof. See [22,23].

(2.11)



In [23], it is also shown that the Gauss laws are satisfied by the Maxwell–Debye system if these laws are also satisfied by the initial fields. In Section 3 we construct operator splitting schemes based on the 2D Maxwell–Debye TE scalar equations derived from (2.5) as

∂H 1 = ∂t µ0 ∂ Ex ∂t ∂ Ey ∂t ∂ Px ∂t ∂ Py ∂t

 ∂ Ex ∂ Ey − , ∂y ∂x   ϵq − 1 1 ∂H 1 = − Ex + Px , ϵ0 ϵ∞ ∂ y τ ϵ0 ϵ∞ τ   ϵq − 1 1 ∂H 1 =− − Ey + Py , ϵ0 ϵ∞ ∂ x τ ϵ0 ϵ∞ τ   ϵ0 ϵ∞ ϵq − 1 1 = Ex − Px , τ τ   ϵ0 ϵ∞ ϵq − 1 1 = Ey − Py . τ τ 

(2.12a)

(2.12b)

(2.12c)

(2.12d) (2.12e)

3. Operator splitting methods for Maxwell’s equations in dispersive media In this section we construct new operator splitting methods for the Maxwell–Debye system (2.12), by extending the ideas from [19] used to construct splitting methods for Maxwell’s equations in free space. The extensions of the operator splitting schemes to dispersive media of Debye type are non-trivial [24] since dispersive media have physical dissipation and dispersion which is absent in the case of wave propagation in free space. 3.1. Discretization Consider the spatial domain Ω = [0, a] × [0, b] ⊂ R2 and time interval [0, T ] with a, b, T > 0 and spatial step sizes 1x > 0 and 1y > 0 and time step 1t > 0. The discretization of the intervals [0, a], [0, b], and [0, T ] is performed as follows. Define L, J , N ∈ N such that L1x = a, J 1y = b and N 1t = T . For ℓ, j, n ∈ N, we consider the discretizations 0 = x 0 ≤ x 1 ≤ · · · ≤ x ℓ ≤ · · · ≤ x L = a, 0 = y0 ≤ y1 ≤ · · · ≤ yj ≤ · · · ≤ yJ = b, 0

1

n

0 = t ≤ t ≤ ··· ≤ t ≤ ··· ≤ t

N

(3.1) (3.2)

= T,

(3.3) γ

where xℓ = ℓ1x, yj = j1y, and t = n1t for 0 ≤ ℓ ≤ L, 0 ≤ j ≤ J, and 0 ≤ n ≤ N. Define (xα , yβ , t ) = (α 1x, β 1y, γ 1t ) where α is either ℓ or ℓ+ 12 , β is either j or j + 12 , and γ is either n or n + 21 with ℓ, j, n ∈ N. The operator splitting schemes (like the Yee scheme) stagger the electric and magnetic fields in space. Fields Ex , Ey , and H are staggered in the x and y directions. We define the discrete meshes n

τhEx :=









(3.4)

2

τh y :=



τhH :=



E



xℓ+ 1 , yj |0 ≤ ℓ ≤ L − 1, 0 ≤ j ≤ J , xℓ , yj+ 1

| 0 ≤ ℓ ≤ L, 0 ≤ j ≤ J − 1 ,   xℓ+ 1 , yj+ 1 |0 ≤ ℓ ≤ L − 1, 0 ≤ j ≤ J − 1 ,

(3.5)

2

2

2

(3.6)

164

V.A. Bokil et al. / Journal of Computational and Applied Mathematics 263 (2014) 160–188

Fig. 1. Staggered spatial grid showing grid points where fields will be discretized. The horizontal axis corresponds to the x-direction and the vertical axis corresponds to the y-direction.

to be the sets of spatial grid points on which the Ex , Ey , and H fields, respectively, will be discretized. The components Px and Py of the polarization are discretized at the same spatial locations as the fields Ex and Ey , respectively. Fig. 1 shows the staggering in space for the grid. For the time discretization, the components Ex , Ey , Px , Py are all discretized at integer time steps t n for 0 ≤ n ≤ N. In the Yee scheme, the magnetic field H is staggered in time with respect to Ex and Ey and discretized 1

at t n+ 2 for 0 ≤ n ≤ N − 1. However, in the splitting methods there is no time staggering and H is also discretized at integer time steps t n , 0 ≤ n ≤ N. Let U be one of the field variables H , Ex , Ey , Px or Py , and α be either ℓ or ℓ + 12 , β be either j or j + 21 , and γ be either n or n +

with ℓ, j, n ∈ N. We define the grid functions or the numerical approximations

1 2

γ Uα,β

≈ U (xα , yβ , t γ ).

We will also use the notation U (t γ ) to denote the continuous solution on the domain Ω at time t γ , and the notation U γ to denote the corresponding grid function on its discrete spatial mesh at time t γ . The boundary values of the electric field given by the PEC condition (2.2) are Exn

1 ℓ+ 2

= Exn

,0

= Eyn

1 ,J ℓ+ 2

0,j+ 1 2

= Eyn

L,j+ 1 2

= 0,

0 ≤ ℓ ≤ L − 1, 0 ≤ j ≤ J − 1, 0 ≤ n ≤ N .

(3.7)

Next, we define some discrete operators in a standard way (see, for e.g. [19]). The centered temporal difference operator is defined as γ+1

γ δt Uα,β

:=

γ−1

Uα,β 2 − Uα,β 2

1t

,

(3.8)

the centered spatial difference operators in the x and y directions, respectively, as γ δx Uα,β γ δy Uα,β

γ α+ 12 ,β

U

:=

U

:=

γ α,β+ 12

γ

− Uα− 1 ,β 2

1x

,

(3.9)

,

(3.10)

γ

− Uα,β− 1 2

1y

and a discrete time averaging operator as γ U α,β

γ+1

=

γ−1

Uα,β 2 + Uα,β 2 2

.

(3.11)

V.A. Bokil et al. / Journal of Computational and Applied Mathematics 263 (2014) 160–188

165

Lastly, we define the discrete L2 grid norms:

 2  2  J −1  L−1    n    n     ∥ ∥ = Vxℓ+ 1 ,j  + Vyℓ,j+ 1  1x1y, Vn 2E

ℓ=0 j=0

∥H n ∥2H =

2

V = E, P,

(3.12)

2

J −1  L−1  2    n Hℓ+ 1 ,j+ 1  1x1y. 2

ℓ=0 j=0

(3.13)

2

3.2. A sequential operator splitting method We first construct a sequential operator splitting method for the two-dimensional Maxwell–Debye model (2.12). The Cauchy problem for the 2D TE Maxwell–Debye model (2.12) can be written in matrix form as follows. Find U(·, t ) defined on Ω × [0, T ] such that

∂U = ΛU, ∂t U(x, 0) = U0 (x);

(3.14a) x ∈ Ω,

(3.14b)

where



1 ∂

0

  1 ∂   ϵ ϵ ∂y  0 ∞  1 ∂ Λ = −  ϵ0 ϵ∞ ∂ x   0  

 

H Ex    U = Ey  ; P  x Py

−

µ0 ∂ y −(ϵq − 1) τ

0

0

µ0 ∂ x 0

−(ϵq − 1) τ

0

ϵ0 ϵ∞ (ϵq − 1) τ

1 ∂

0

ϵ0 ϵ∞ (ϵq − 1) τ

0 1

ϵ0 ϵ∞ τ 0

−

1

τ

0

0



  0    1  . ϵ0 ϵ∞ τ    0   1 − τ

(3.15)

We will assume that the PEC boundary conditions (2.2) hold. To construct the sequential operator splitting scheme, we first decompose operator Λ so that (3.14a) becomes

∂U 1 = AU + BU, ∂t τ

(3.16)

where the sub-operators A, B in (3.16) are defined as



1 ∂

0

   1 ∂   A =  ϵ0 ϵ∞ ∂ y  − 1 ∂  ϵ0 ϵ∞ ∂ x 

−

µ0 ∂ y

0 0

0

1 ∂

µ0 ∂ x 0

0 0

0

0

0

0 0

0 0

0 0



0

   0  ,  0  

(3.17)

0 0

and

0 0   B= 0  

0 0

0

0

−(ϵq − 1)

0

ϵ0 ϵ∞

0

−(ϵq − 1)

0

0

−1

ϵ0 ϵ∞ (ϵq − 1) 0

ϵ0 ϵ∞ (ϵq − 1)

0 1

0

0



0

   1 .  ϵ0 ϵ∞  

(3.18)

0

−1

This decomposition separates terms containing the relaxation parameter τ from the other terms. For the sequential operator splitting method, we further decompose A into A1 and A2 so that

∂U 1 = A1 U + A2 U + BU, ∂t τ

(3.19)

166

V.A. Bokil et al. / Journal of Computational and Applied Mathematics 263 (2014) 160–188

Fig. 2. A visualization of the SQ-DB scheme over the time interval [t n , t n+1 ]. Each step communicates with its subsequent steps through initial conditions.

where



0

  1 ∂  A1 =   ϵ0 ϵ∞ ∂ y  0  0 0

1 ∂





µ0 ∂ y

0

0

0

0

0

0

0 

0 0 0

0 0 0

0 0 0

0  0 0

  ;

0

  0  A2 =  − 1 ∂  ϵ ϵ ∂x  0 ∞ 0 0

0 0

−

1 ∂

0

0

0

0 

0

0

0 

µ0 ∂ x

0 0 0



0

0 0

0 0

 .

(3.20)



0 0

Each Ai (i = 1, 2) operator represents wave propagation along a one-dimensional direction (x or y) in a physically non-dispersive dielectric and thus constitutes dimensional splitting of the operator A [25]. The terms involving τ remain separated from the terms involving spatial derivatives. 3.2.1. The semi-discrete sequential operator splitting scheme SQ-DB The sequential operator splitting scheme SQ-DB on discrete time intervals [t n , t n+1 ], 0 ≤ n ≤ N is defined as follows. (SQ-DB): Given U(t 0 ) = H (t 0 ), Ex (t 0 ), Ey (t 0 ), Px (t 0 ), Py (t 0 ) for n = 0 . . . N − 1 



T

on Ω

Step A1 : Solve for U˜ (t n+1 ) = H˜ (t n+1 ), E˜ x (t n+1 ), E˜ y (t n+1 ), P˜x (t n+1 ), P˜y (t n+1 ) ∂ U˜ (t ) = A1 U˜ (t ); ∂t

T

on Ω × (t n , t n+1 ] in:

˜ (t n ) = U(t n ). U

(3.21)

 T Step A2 : Solve for Uˆ (t n+1 ) = Hˆ (t n+1 ), Eˆ x (t n+1 ), Eˆ y (t n+1 ), Pˆx (t n+1 ), Pˆy (t n+1 ) on Ω × (t n , t n+1 ] in: ∂ Uˆ (t ) = A2 Uˆ (t ); ∂t

ˆ (t n ) = U˜ (t n+1 ). U

(3.22)

 T Step B: Solve for U(t n+1 ) = H (t n+1 ), Ex (t n+1 ), Ey (t n+1 ), Px (t n+1 ), Py (t n+1 ) on Ω × (t n , t n+1 ] in: ∂ U(t ) 1 = BU(t ); ∂t τ

ˆ (t n+1 ). U(t n ) = U

(3.23)

end. The perfect conductor boundary conditions (2.2) are assumed to hold in the above scheme. The sub-steps in the sequential splitting scheme SQ-DB communicate through their initial conditions. Fig. 2 shows the sequence of steps for the interval [t n , t n+1 ]. The main relaxation of water at room temperature in the microwave range of frequencies is τ = 8.1 × 10−12 s, and the two relative permittivities are ϵ∞ = 1, and ϵs = 78.2 [9]. In [3] it is shown that the Maxwell–Debye model for water has a boundary layer of size O (c∞ τ ) ≈ O (10−3 ) m in which the model solution has a fast decaying transient term travelling with the fastest speed of harmonic waves c∞ . Thus, the small value of τ results in a small boundary layer contributing to the stiffness of the system. It is discussed in [26] that for operator splitting schemes with one stiff and one non-stiff sub-operator, the non-stiff operator must precede the stiff operator. This is because the exact solution and the splitting schemes for which the stiff

V.A. Bokil et al. / Journal of Computational and Applied Mathematics 263 (2014) 160–188

167

operator ends the iteration always preserve the reduced (up to first order in τ ) structure of the coupled system. The stiff operator alone determines this reduced structure regardless of what the non-stiff operator is and acts as a projection onto the reduced model. The integration of the non-stiff operator does not ensure that the solution remains near the reduced model unless the two operators commute [26]. We consider the operator B, containing all the τ related terms to be the stiff-operator for this splitting and the operator A to be the non-stiff operator. Thus, the steps with the sub-operators A1 and A2 precede the step involving the sub-operator B. 3.2.2. The fully-discrete sequential operator splitting scheme SQ-DB-D We define the discrete mesh E

E

Ωh = τhH × τhEx × τh y × τhEx × τh y ,

(3.24)

and use the implicit Crank–Nicolson scheme to fully discretize steps A1 , A2 and B. The resulting discretized sequential operator splitting scheme SQ-DB-D is as follows. (SQ-DB-D):

T

Given U0 = H 0 , Ex0 , Ey0 , Px0 , Py0 on Ωh for n = 0 . . . N − 1 



Step A1 : Solve for U˜ n+1 = H˜ n+1 , E˜ xn+1 , E˜ yn+1 , P˜xn+1 , P˜yn+1 ˜ n+11 H

n − Hℓ+ 1 ,j + 1

ℓ+ 2 ,j+ 21

2

2

1t E˜ xn+11 − Exn

1 ,j ℓ+ 2

ℓ+ 2 ,j

=

1t

=

1 2µ0

1 2ϵ0 ϵ∞

 δy E˜ xn+11

ℓ+ 2 ,j+ 1 2

T



+ Exn

1 ,j+ 1 ℓ+ 2 2

on Ωh in:

,

(3.25a)

  n +1 n + H , δy H˜ ℓ+ 1 ℓ+ 1 ,j ,j

(3.25b)

2

2

E˜ yn+1 1 = Eyn

,

(3.25c)

P˜ xn+11 = Pxn

,

(3.25d)

P˜ yn+1 1 = Pyn

.

(3.25e)

ℓ,j+ 1 2

ℓ,j+ 2

,j ℓ+ 1 2

ℓ+ 2 ,j

ℓ,j+ 1 2

ℓ,j+ 2

 T Step A2 : Solve for Uˆ n+1 = Hˆ n+1 , Eˆ xn+1 , Eˆ yn+1 , Pˆxn+1 , Pˆyn+1 on Ωh in: n +1 − H˜ ℓ+ 1 ,j + 1

ˆ n+11 H

ℓ+ 2 ,j+ 21

2

2

1t Eˆ yn+1 1 − E˜ yn+1 1 ℓ,j+ 2

ℓ,j+ 2

1t

=−

=−

1 2µ0

1 2ϵ0 ϵ∞

δx



Eˆ yn+1

1 ,j+ 1 ℓ+ 2 2

+ E˜ yn+1



1 ,j+ 1 ℓ+ 2 2

,

(3.26a)

  δx Hˆ ℓ,n+j+1 1 + H˜ ℓ,n+j+1 1 , 2

(3.26b)

2

Eˆ xn+11 = E˜ xn+11 ,

(3.26c)

Pˆ xn+11 = P˜ xn+11 ,

(3.26d)

Pˆ yn+1 1 = P˜ yn+1 1 .

(3.26e)

ℓ+ 2 ,j

ℓ+ 2 ,j

ℓ+ 2 ,j

ℓ+ 2 ,j

ℓ,j+ 2

ℓ,j+ 2

 T Step B: Solve for Un+1 = H n+1 , Exn+1 , Eyn+1 , Pxn+1 , Pyn+1 on Ωh in: Exn+11 − Eˆ xn+11 ℓ+ 2 ,j

ℓ+ 2 ,j

1t Pxn+11 − Pˆ xn+11 ℓ+ 2 ,j

ℓ+ 2 ,j

1t

=

ϵ0 ϵ∞ (ϵq − 1) n+1 1 = (Ex 1 + Eˆ xn+11 ) − ℓ+ 2 ,j ℓ+ 2 ,j 2τ 2τ

Eyn+1 1 − Eˆ yn+1 1 ℓ,j+ 2

ℓ,j+ 2

1t

−(ϵq − 1) n+1 1 (Ex 1 + Eˆ xn+11 ) + ℓ+ 2 ,j ℓ+ 2 ,j 2τ 2τ ϵ0 ϵ∞

=



−(ϵq − 1) n+1 1 (Ey 1 + Eˆ yn+1 1 ) + ℓ,j+ 2 ℓ,j+ 2 2τ 2τ ϵ0 ϵ∞



Pxn+11 + Pˆ xn+11 ℓ+ 2 ,j

Pxn+11 ℓ+ 2 ,j



ℓ+ 2 ,j

+ Pˆxn+1



ℓ+ 1 ,j 2

Pyn+1 1 + Pˆ yn+1 1 ℓ,j+ 2



ℓ,j+ 2

,

, 

(3.27a)

(3.27b)

,

(3.27c)

168

V.A. Bokil et al. / Journal of Computational and Applied Mathematics 263 (2014) 160–188

Pyn+1 1 − Pˆ yn+1 1 ℓ,j+ 2

ϵ0 ϵ∞ (ϵq − 1) n+1 1 = (Ey 1 + Eˆ yn+1 1 ) − ℓ,j+ 2 ℓ,j+ 2 2τ 2τ

ℓ,j+ 2

1t H n+11

ℓ+ 2 ,j+ 21



Pyn+1 1 ℓ,j+ 2



+ Pˆyn+1

ℓ,j+ 1 2

,

(3.27d)

n +1 = Hˆ ℓ+ . 1 ,j + 1 2

(3.27e)

2

end The perfect conductor boundary conditions (3.7) are assumed to hold for the above scheme. Below we provide details of the implementations of all three steps of the SQ-DB-D scheme. In particular, we show that the update equations for Step A1 and Step A2 consist of a tridiagonal matrix solve for the electric field followed by an explicit update for the magnetic field using the most recently computed values of the electric field. The implementation of Step B involves a matrix–vector product computation. We define two Courant numbers

ηw :=

c∞ 1t

1w

,

w = x, y.

(3.28)

3.2.3. Implementation of the sequential scheme SQ-DB-D

Step A1 : On the interval (t n , t n+1 ] solve for E˜ xn+1 in the tridiagonal system   1  0 ··· 0  0    E˜ xn+11 2 2 2 0 ℓ+ , 0 ηy −ηy  −ηy  2    n + 1 ··· 0   E˜ 1+   b1   4 x 1    2 4    ℓ+ 2 ,1   ..    .     . .. .. ..   ..  =  . , .  . .  0 0         ˜ n +1   ..   Ex 1 −ηy2 ηy2 −ηy2        0   ℓ+ 2 ,J −1  b J −1 ··· 1+  4 2 4  ˜ n +1 ···

0

0

0

Ex

1

(3.29)

0

,J ℓ+ 1 2

where 0 ≤ ℓ ≤ L − 1, and for 0 ≤ j ≤ J − 1, bj =

Exn 1 ℓ+ 2 ,j





+ c∞ ηy µ0 Hℓ+ 1 ,j+ 1 − Hℓ+ 1 ,j− 1 + n

n

2

2

2

ηy2

2



4

Exn 1 ℓ+ 2 ,j+1

−

2Exn 1 ℓ+ 2 ,j

−

Exn 1 ℓ+ 2 ,j

+

Exn 1 ℓ+ 2 ,j−1



.

(3.30)

˜ n+1 in Next, solve explicitly for H ˜ n+11 H

n

ℓ+ 2 ,j+ 21

= Hℓ+ 1 ,j+ 1 2

2

1t + 2µ0 1y



E˜ xn+1

1 ,j+1 ℓ+ 2



,

(3.31)

  0 0,j+ 2    ˆ n+1   b1  0  Ey 1   .  1,j+  2   .   ..  =  .  ,     .   0   ..  2   ˆ n +1    . −ηx  Ey     L−1,j+ 12  bL−1  4 ˆEyn+1 0

(3.32)

− E˜ xn+11 + Exn

,j+1 ℓ+ 1 2

ℓ+ 2 ,j

where 0 ≤ ℓ ≤ L − 1 and 0 ≤ j ≤ J − 1.

Step A2 : On the interval (t n , t n+1 ] solve for Eˆ yn+1 in the tridiagonal system  1 2  −ηx   4    0    0

0

 1+

..

η

2 x



−η

2

..

. −ηx2

···

.. 

4 0

···

0

···

4

.

 n+1  0  Eˆ y 1

···

0 2 x

1+

.  ηx2 2

0

1

L,j+ 1 2

with 0 ≤ j ≤ J − 1, and for 0 ≤ ℓ ≤ L − 1, bℓ = E˜ yn+1 1 − c∞ ηx µ0 ℓ,j+ 2



˜ n+11 H

ℓ+ 2 ,j+ 12

n +1 − H˜ ℓ− 1 ,j + 1 2

 +

ηx2

2

4



E˜ yn+1

1 ℓ+1,j+ 2

− 2E˜ yn+1 1 + E˜ yn+1

ℓ−1,j+ 1 2

ℓ,j+ 2



.

(3.33)

ˆ n+1 in Next, solve explicitly for H ˆ n+11 H

ℓ+ 2 ,j+ 21

= H˜ n+11

ℓ+ 2 ,j+ 12

1t − 2µ0 1x

for 0 ≤ ℓ ≤ L − 1 and 0 ≤ j ≤ J − 1.



Eˆ yn+1

1 ℓ+1,j+ 2

− Eˆ yn+1 1 + E˜ yn+1 ℓ,j+ 2

ℓ+1,j+ 1 2

− E˜ yn+1

1 ℓ,j+ 2



,

(3.34)

V.A. Bokil et al. / Journal of Computational and Applied Mathematics 263 (2014) 160–188

169

We note that the tridiagonal matrices in Step A1 and Step A2 are constant in time, and thus an LU factorization can be utilized outside the time loop. Further, on a uniform mesh in which 1x = 1y the two Courant numbers are equal ηx = ηy , and the two tridiagonal matrices in these two steps are the same.

Step B: Using Eqs. (3.27a) and (3.27b), we can write Exn+11 and Pxn+11 explicitly as ℓ+ 2 ,j

1t (ϵq − 1) 1t 2 (ϵq − 1) − 2τ 2τ (2τ + 1t )

ℓ+ 2 ,j

 1t (ϵq − 1) 1t 2 (ϵq − 1) + Eˆ xn+11 ℓ+ 2 ,j ℓ+ 2 ,j 2τ 2τ (2τ + 1t )   2 1t 1t 1t + + − Pˆ n+1 , 2τ ϵ0 ϵ∞ ϵ0 ϵ∞ (2τ + 1t ) 2τ ϵ0 ϵ∞ (2τ + 1t ) xℓ+ 12 ,j     1t 1t 2 (ϵq − 1) 21t ϵ0 ϵ∞ (ϵq − 1) 1+ Pxn+11 − = Eˆ xn+11 ℓ+ 2 ,j ℓ+ 2 ,j 2τ 2τ (2τ + 1t (ϵq − 1)) 2τ + 1t (ϵq − 1)   1t 1t 2 (ϵq − 1) + 1− − Pˆ xn+11 . ℓ+ 2 ,j 2τ 2τ (2τ − 1t (ϵq − 1)) 

Exn+11

1+





= 1−

(3.35)

(3.36)

With some algebraic manipulation, we arrive at the system of equations

   1 n +1 ˆ = 1 − α(−(ϵq − 1)) Ex 1 + −α Pˆ xn+11 , ℓ+ 2 ,j ℓ+ 2 ,j ϵ0 ϵ∞   n+1 n + 1 = −αϵ0 ϵ∞ (ϵq − 1) Eˆ x + (1 − α(−1)) Pˆx ,

Exn+11 ℓ+ 2 ,j



Pxn+11



1 ,j ℓ+ 2

ℓ+ 2 ,j

1 ,j ℓ+ 2

−1

with 0 ≤ ℓ ≤ L − 1, and 0 ≤ j ≤ J, and α = −21t 2τ + ϵq 1t . An identical approach is used to show that similar results hold for the Ey and Py components in Step B. Thus, Step B of the SQ-DB-D scheme on the interval (t n , t n+1 ] can be written as the computation of the matrix–vector product





Exn+11



ℓ+ 2 ,j



Pxn+11



=C 

ℓ+ 2 ,j

Eˆ xn+11

ℓ+ 2 ,j

Pˆ xn+11





 ,



Eyn+1 1



Pyn+1 1

=C

ℓ,j+ 2



ℓ,j+ 2

ℓ+ 2 ,j

Eˆ yn+1 1

ℓ,j+ 2

Pˆ yn+1 1

 ,

(3.37)

ℓ,j+ 2

for 0 ≤ ℓ ≤ L − 1, and 0 ≤ j ≤ J − 1, where C = I − α A, I is the 2 × 2 identity matrix, and

 A=

−(ϵq − 1)

(ϵq − 1)ϵ0 ϵ∞

1



ϵ0 ϵ∞ . −1

(3.38)

The update equations for the SQ-DB-D scheme consist of Eqs. (3.29)–(3.38). In Section 4 we prove the unconditional stability of the SQ-DB-D scheme, and in Section 6 we show that the scheme is first order in time and second order in space. 3.3. A symmetrized operator splitting scheme To obtain a second order accurate scheme in time, we use the ideas of symmetrization of splitting schemes as first proposed in [27]. To construct a symmetrized operator splitting scheme, we decompose (3.16) as

∂U 1 1 1 1 1 = A1 U + A2 U + BU + A2 U + A1 U. ∂t 2 2 τ 2 2 The ordering 21 A1 →

1 A 2 2

(3.39)

→ B → 12 A2 → 21 A1 results in a scheme that is second order accurate in time [25].

3.3.1. The semi-discrete symmetrized operator splitting scheme SY-DB The symmetrized operator splitting scheme SY-DB for the Maxwell–Debye system is defined as follows. (SY-DB):

T

Given U(t 0 ) = H (t 0 ), Ex (t 0 ), Ey (t 0 ), Px (t 0 ), Py (t 0 ) for n = 0 . . . N − 1



on Ω

 T ˜ (t n+ 12 ) = H˜ (t n+ 12 ), E˜ x (t n+ 21 ), E˜ y (t n+ 21 ), P˜x (t n+ 12 ), P˜y (t n+ 12 ) on Ω × (t n , t n+ 12 ] in: Step 21 A1 : Solve for W ˜ (t ) ∂W ˜ (t ); = A1 W ∂t

˜ (t n ) = U(t n ). W

(3.40)

170

V.A. Bokil et al. / Journal of Computational and Applied Mathematics 263 (2014) 160–188

Fig. 3. A visualization of the SY-DB scheme over the time interval [t n , t n+1 ]. Each step communicates with its subsequent steps through initial conditions.

 T 1 1 1 1 1 1 1 Step 12 A2 : Solve for W(t n+ 2 ) = H (t n+ 2 ), Ex (t n+ 2 ), Ey (t n+ 2 ), Px (t n+ 2 ), Py (t n+ 2 ) on Ω × (t n , t n+ 2 ] in: ∂ W(t ) = A2 W(t ); ∂t



1

˜ t n+ 2 W(t n ) = W



.

(3.41)

 T Step B: Solve for U˜ (t n+1 ) = H˜ (t n+1 ), E˜ x (t n+1 ), E˜ y (t n+1 ), P˜x (t n+1 ), P˜y (t n+1 ) on Ω × (t n , t n+1 ] in: ∂ U˜ (t ) 1 = BU˜ (t ); ∂t τ



1

˜ ( t n ) = W t n+ 2 U



.

(3.42)

 T 1 Step 12 A2 : Solve for Uˆ (t n+1 ) = Hˆ (t n+1 ), Eˆ x (t n+1 ), Eˆ y (t n+1 ), Pˆx (t n+1 ), Pˆy (t n+1 ) on Ω × (t n+ 2 , t n+1 ] in: ∂ Uˆ (t ) = A2 Uˆ (t ); ∂t



1

ˆ t n+ 2 U



= U˜ (t n+1 ).

(3.43)

 T 1 Step 21 A1 : Solve for U(t n+1 ) = H (t n+1 ), Ex (t n+1 ), Ey (t n+1 ), Px (t n+1 ), Py (t n+1 ) on Ω × (t n+ 2 , t n+1 ] in: ∂ U(t ) = A1 U(t ); ∂t



1

U t n+ 2



= Uˆ (t n+1 ).

(3.44)

end. Fig. 3 displays this sequence of steps. 3.3.2. The fully-discrete symmetrized operator splitting scheme SY-DB-D The fully discretized symmetrized operator splitting scheme SY-DB-D uses an implicit time discretization similar to the SQ-DB-D scheme. (SY-DB-D): Given U0 = H 0 , Ex0 , Ey0 , Px0 , Py0 for n = 0 . . . N − 1



T

on Ωh

  1 1 1 1 T ˜ n+ 21 = H˜ n+ 12 , E˜ xn+ 2 , E˜ yn+ 2 , P˜xn+ 2 , P˜yn+ 2 on Ωh in: Step 12 A1 : Solve for W ˜ H

n+ 21

ℓ+ 12 ,j+ 21

n − Hℓ+ 1 ,j + 1 2

1t n+ 12

E˜ x

1 ,j ℓ+ 2

− Exn

1 ,j ℓ+ 2

1t

=

2

=

1 4µ0

1 4ϵ0 ϵ∞

 n+ 1 δy E˜ x 21

ℓ+ 2 ,j+ 1 2

+ Exn

1 ,j+ 1 ℓ+ 2 2

  n+ 1 n δy H˜ ℓ+ 12 ,j + Hℓ+ , 1 ,j 2

2



,

(3.45a)

(3.45b)

V.A. Bokil et al. / Journal of Computational and Applied Mathematics 263 (2014) 160–188 n+ 12

E˜ y

ℓ,j+ 1 2

n+ 12

P˜ x

ℓ+ 1 ,j 2

n+ 21

P˜ y

ℓ,j+ 1 2

171

= Eyn

,

(3.45c)

= Pxn

,

(3.45d)

= Pyn

.

(3.45e)

ℓ,j+ 1 2

ℓ+ 1 ,j 2

ℓ,j+ 1 2

 T 1 1 n+ 1 n+ 1 n+ 1 n+ 1 Step 12 A2 : Solve for Wn+ 2 = H n+ 2 , Ex 2 , Ey 2 , Px 2 , Py 2 on Ωh in: n+ 21

H

1

n+ − H˜ ℓ+ 12 ,j+ 1

ℓ+ 12 ,j+ 21

2

2

1t n+ 21

Ey

ℓ,j+ 1 2

n+ 21 ℓ+ 1 ,j 2

n+ 12

Px

,j ℓ+ 1 2

n+ 21

Py

ℓ,j+ 1 2

(3.46a)

1

n+ − E˜ y 2 1

  −1 n+ 21 n+ 12 ˜ = δx Hℓ,j+ 1 + Hℓ,j+ 1 , 4ϵ0 ϵ∞ 2 2

ℓ,j+ 2

1t Ex

  −1 n+ 21 n+ 21 ˜ = δ x Ey 1 1 + Ey 1 1 , ℓ+ 2 ,j+ 2 ℓ+ 2 ,j+ 2 4µ0

(3.46b)

1

n+ = E˜ x 21 ,

(3.46c)

ℓ+ 2 ,j

n+ 12

= P˜x

1 ,j ℓ+ 2

,

(3.46d)

1

n+ = P˜y 2 1 .

(3.46e)

ℓ,j+ 2



Step B: Solve for U˜ n+1 = H˜ n+1 , E˜ xn+1 , E˜ yn+1 , P˜xn+1 , P˜yn+1 n+ 12

E˜ xn+11 − Ex ℓ+ 2 ,j

1 ,j ℓ+ 2

1t ℓ+ 2 ,j

,j ℓ+ 1 2

1t n+ 12

E˜ yn+1 1 − Ey ℓ,j+ 2

ℓ,j+ 1 2

1t n+ 12

P˜ yn+1 1 − Py ℓ,j+ 2

ℓ,j+ 1 2

1t

on Ωh in:

=

    1 −(ϵq − 1) ˜ n+1 n+ 1 n+ 1 + , Ex 1 + Ex 21 P˜ xn+11 + Px 21 ℓ+ 2 ,j ℓ+ 2 ,j ℓ+ 2 ,j ℓ+ 2 ,j 2τ 2τ ϵ0 ϵ∞

(3.47a)

=

    ϵ0 ϵ∞ (ϵq − 1) ˜ n+1 1 n+ 1 n+ 1 − , Ex 1 + Ex 21 P˜ xn+11 + Px 21 ℓ+ 2 ,j ℓ+ 2 ,j ℓ+ 2 ,j ℓ+ 2 ,j 2τ 2τ

(3.47b)

    −(ϵq − 1) ˜ n+1 1 n+ 21 n+ 12 n +1 ˜ = Ey 1 + Ey 1 + Py 1 + Py 1 , ℓ,j+ 2 ℓ,j+ 2 ℓ,j+ 2 ℓ,j+ 2 2τ 2τ ϵ0 ϵ∞

(3.47c)

    ϵ0 ϵ∞ (ϵq − 1) ˜ n+1 1 n+ 21 n+ 12 n +1 ˜ = Ey 1 + Ey 1 − Py 1 + Py 1 , ℓ,j+ 2 ℓ,j+ 2 ℓ,j+ 2 ℓ,j+ 2 2τ 2τ

(3.47d)

n+ 12

P˜ xn+11 − Px

T

n+ 1

= Hℓ+ 12 ,j+ 1 .

˜ n+11 H

ℓ+ 2 ,j+ 12

2

(3.47e)

2



Step 12 A2 : Solve for Uˆ n+1 = Hˆ n+1 , Eˆ xn+1 , Eˆ yn+1 , Pˆxn+1 , Pˆyn+1 ˆ n+11 H

n +1 − H˜ ℓ+ 1 ,j + 1

ℓ+ 2 ,j+ 21

2

1t Eˆ yn+1 1 − E˜ yn+1 1 ℓ,j+ 2

ℓ,j+ 2

1t

=

2

T

  −1 n+1 n +1 ˆ ˜ = δx Ey 1 1 + Ex 1 1 , ℓ+ 2 ,j+ 2 ℓ+ 2 ,j+ 2 4µ0

  −1 δx Hˆ ℓ,n+j+1 1 + H˜ ℓ,n+j+1 1 , 4ϵ0 ϵ∞ 2 2

on Ωh in:

(3.48a)

(3.48b)

Eˆ xn+11 = E˜ xn+11 ,

(3.48c)

Pˆ xn+11 = P˜ xn+11 ,

(3.48d)

Pˆ yn+1 1 = P˜ yn+1 1 .

(3.48e)

ℓ+ 2 ,j

ℓ+ 2 ,j

ℓ,j+ 2

ℓ+ 2 ,j

ℓ+ 2 ,j

ℓ,j+ 2

172

V.A. Bokil et al. / Journal of Computational and Applied Mathematics 263 (2014) 160–188

T  Step 12 A1 : Solve for Un+1 = H n+1 , Exn+1 , Eyn+1 , Pxn+1 , Pyn+1 on Ωh in: H n+11

ℓ+ 2 ,j+ 21

n+1 − Hˆ ℓ+ 1 ,j + 1 2

2

1t Exn+11 − Eˆ xn+11 ℓ+ 2 ,j

ℓ+ 2 ,j

=

1t

1

=

4µ0

 δy Exn+11

ℓ+ 2 ,j+ 1 2



1

δy Hℓ+ 1 ,j

4ϵ0 ϵ∞

n+1

ℓ+ 2 ,j+ 1 2

+ Hˆ n+11



ℓ+ 2 ,j

2



+ Eˆ xn+11

,

(3.49a)

,

(3.49b)

Eyn+1 1 = Eˆ yn+1 1 ,

(3.49c)

Pxn+11 = Pˆ xn+11 ,

(3.49d)

Pyn+1 1 = Pˆ yn+1 1 .

(3.49e)

ℓ,j+ 2

ℓ,j+ 2

ℓ+ 2 ,j

ℓ+ 2 ,j

ℓ,j+ 2

ℓ,j+ 2

end. We assume the PEC boundary conditions (3.7) hold in the above scheme. 3.4. Implementation of the symmetrized splitting scheme SY-DB-D In this section we explain how the fully discrete scheme SY-DB-D is implemented in detail. As in the case of the SQ-DB-D scheme, we note that in the case of a uniform mesh the tridiagonal matrices involved in all the steps below are identical. Thus, an LU factorization of these matrices can be utilized outside the time loop. 1

1 n+ Step 12 A1 : On the interval (t n , t n+ 2 ], solve for E˜ x 2 in the tridiagonal system

 1  −ηy2   16     0    0 

0 1+



−ηy2

8

..

..

.

. −ηy2

.. 1+

0





n+ 12

E˜ x



ℓ+ 1 ,0 2

   0    n+ 1    b˜ 1  0    E˜ xℓ+ 21 ,1      ..  2      ..  = .  ,  0  . . ..       2   n+ 1  −ηy  E˜ x 21  b˜  ℓ+ 2 ,J −1  J −1 16     0 n+ 12 1 E˜ x 1

.  ηy2



16

···

0

···

16

···

0

···

0

ηy2



8

0

(3.50)

ℓ+ 2 ,J

where 0 ≤ ℓ ≤ L − 1, and for 0 ≤ j ≤ J − 1, b˜ j = Exn

ℓ+ 1 ,j 2

+

c∞ ηy µ0 

 η2  y n + Exn − Hℓ+ 1 ,j − 1

Hn

ℓ+ 21 ,j+ 12

2

2

n+ 12

Next, using the newly computed values of E˜ x

˜ H

n+ 12

ℓ+ 12 ,j+ 21

n = Hℓ+ + 1 ,j + 1 2

2

1t 4µ0 1y



n+ 21

E˜ x

1 ,j+1 ℓ+ 2

ℓ+ 1 ,j+1 2

16

2

− 2Exn

ℓ+ 1 ,j 2



+ Exn

ℓ+ 1 ,j−1 2

.

(3.51)

1

˜ n+ 2 in solve explicitly for the values of H n+ 12



,

(3.52)

 0        b˜ 1  0  E   .     .     .    ..  =  . ,  .  0     n+ 1    ..  −ηx2   2     E   yL−1,j+ 1  b˜ L−1 2  16  1

(3.53)

− E˜ x

ℓ+ 1 ,j 2

+ Exn

ℓ+ 1 ,j+1 2

− Exn

ℓ+ 1 ,j 2

with 0 ≤ ℓ ≤ L, 0 ≤ j ≤ J. 1

1 n+ Step 12 A2 : On the interval (t n , t n+ 2 ], solve for E˜ x 2 in the tridiagonal system

  1 2  −ηx   16    0    0 0

0

 1+

..

ηx2 8

.

··· ···



0

···

−ηx2

···

16

..

. −ηx2 16 0

..  1+

0 

.  ηx2

0

8

1

n+ 12

Ey



0,j+ 1 2 n+ 12 y 1,j+ 1 2

n+ 2

Ey

L,j+ 1 2

0

V.A. Bokil et al. / Journal of Computational and Applied Mathematics 263 (2014) 160–188

173

where 0 ≤ j ≤ J − 1, and for 0 ≤ ℓ ≤ L − 1 n+ 21

b˜ ℓ = E˜ y

ℓ,j+ 1 2

−

c∞ ηx µ0



˜ H

n+ 12

n+ 1

− H˜ ℓ− 12 ,j+ 1

ℓ+ 12 ,j+ 12

2

2

n+ 21

Next, using the computed values of Ey H

n+ 12

n+ 1

ℓ+ 21 ,j+ 12

= H˜ ℓ+ 12 ,j+ 1 − 2

2

1t 4µ0 1x



ηx2

 +



16

2

n+ 12

E˜ y

n+ 12

ℓ+1,j+ 1 2



n+ 12

− 2E˜ y

ℓ,j+ 1 2

+ E˜ y

ℓ−1,j+ 1 2

.

(3.54)

1

, solve for H n+ 2 explicitly in

n+ 21

n+ 21

Ey

ℓ+1,j+ 1 2

− Ey

n+ 12

+ E˜ y

1 ℓ,j+ 2

ℓ+1,j+ 1 2



n+ 12

− E˜ y

ℓ,j+ 1 2

,

(3.55)

with 0 ≤ ℓ ≤ L, 0 ≤ j ≤ J. Step B: As in the case of the SQ-DB-D scheme, on the time interval (t n , t n+1 ] we compute the matrix–vector product  n+ 1   n+ 1   n +1   n +1  E˜ y Ey 2 E˜ x Ex 2 (3.56) ˜ n+1 = C n+ 1 , ˜ n+1 = C n+ 1 , Px

Py

2

Px

2

Py

with matrix C = I − α A, I the 2 × 2 identity matrix, A as defined in (3.38), and α = −21t 2τ + ϵq 1t



Step 12 A2 : On the interval (t  1 2  −ηx   16    0    0

n+ 21

···

ηx

−ηx2

8

16

···

 2

1+

..

. −ηx2

···



16 0

···

0

..

..

.

.

, t n+1 ] solve for E˜ yn+1 in the tridiagonal system

0

0



−1

1+

 n +1  0  Eˆ y 1

  0 0,j+ 2   n + 1  b1  0  Eˆ y 1    .  1,j+  2   .    .    ..   =  . , 0  .     ˆ n +1    ..  −ηx2   Ey     1   L−1,j+ 2  bL−1 16 ˆEyn+1 0

.  ηx2 8

0

1

(3.57)

L,j+ 1 2

where 0 ≤ j ≤ J − 1, and for 0 ≤ ℓ ≤ L − 1 bℓ = E˜ yn+1 1 −

c∞ ηx µ0



n +1 − H˜ ℓ− 1 ,j + 1

ℓ+ 2 ,j+ 12

2

ℓ,j+ 2

˜ n+11 H

2

ηx2

 +



ℓ+1,j+ 1 2

16

2

E˜ yn+1

− 2E˜ yn+1 1 + E˜ yn+1 ℓ,j+ 2



ℓ−1,j+ 1 2

.

(3.58)

ˆ n+1 explicitly in Using the computed values of Eˆ yn+1 , solve for H ˆ n+11 H

ℓ+ 2 ,j+ 12

= H˜ n+11

ℓ+ 2 ,j+ 21

1t − 4µ0 1x



Eˆ yn+1

1 ℓ+1,j+ 2

− Eˆ yn+1 1 + E˜ yn+1

ℓ+1,j+ 1 2

ℓ,j+ 2

− E˜ yn+1



ℓ,j+ 1 2

,

(3.59)

with 0 ≤ ℓ ≤ L, 0 ≤ j ≤ J. 1

Step 12 A1 : On the interval (t n+ 2 , t n+1 ], solve for Exn+1 in the tridiagonal system   1 0 ··· 0  E n +1    0   x 1 2 2 2 0 ℓ+ 2 ,0 ηy −ηy   −ηy     n + 1 1+ ··· 0  E   b1   16  8 16   xℓ+ 12 ,1     ..      .   . .. .. ..   = .. , . . .   0   .  0   ..       ηy2 −ηy2 −ηy2  Exn+11     ℓ+ 2 ,J −1  bJ −1   0 1+ ···  n+1 16 8 16  Ex

···

0

0

0

0

ℓ+ 1 ,J 2

1

(3.60)

where 0 ≤ ℓ ≤ L − 1, and for 0 ≤ j ≤ J − 1 bj = Eˆ xn+11 +

c∞ ηy µ0



ℓ+ 2 ,j+ 12

2

ℓ+ 2 ,j

ˆ n+11 H

n +1 − Hˆ ℓ+ 1 ,j − 1 2

2

 +

ηy2 16



Eˆ xn+11

ℓ+ 2 ,j+1

− 2Eˆ xn+11 + Eˆ xn+11 ℓ+ 2 ,j

ℓ+ 2 ,j−1



.

(3.61)

Using the computed values of Exn+1 , for 0 ≤ ℓ ≤ L, 0 ≤ j ≤ J solve for H n+1 explicitly in H n+11

ℓ+ 2 ,j+ 12

n +1 = Hˆ ℓ+ + 1 ,j+ 1 2

2

1t 4µ0 1y



Exn+11

ℓ+ 2 ,j+1

− Exn+11 + Eˆ xn+11 ℓ+ 2 ,j

ℓ+ 2 ,j+1

− Eˆ xn+11

The update equations for the SY-DB-D scheme consist of (3.50)–(3.62).

ℓ+ 2 ,j



.

(3.62)

174

V.A. Bokil et al. / Journal of Computational and Applied Mathematics 263 (2014) 160–188

4. Stability analysis In this section we prove one of the main advantages of the operator splitting schemes over the conditionally stable Yee scheme, namely their unconditional stability. To show stability of our new methods, we will use the energy technique to show that the energy decay in the continuous problem is maintained in the discrete schemes. Before we prove stability, we need the following lemma [19]. Ey

Lemma 4.1. Let the field variables U , V , and W be defined on the discrete meshes, τh x , τh , and τhH , respectively. Also, assume that the variables U and V satisfy the perfect conducting boundary condition (3.7), i.e. E

Uℓ+ 1 ,0 = Uℓ+ 1 ,J = 0, 2

0 ≤ ℓ ≤ L − 1,

2

V0,j+ 1 = VL,j+ 1 = 0, 2 2

(4.1)

0 ≤ j ≤ J − 1.

(4.2)

Then for any integers 0 < ℓ < L − 1 and 0 < j < J − 1, we have J −1 

2

j =0 I −1 

ℓ=0

J −1 

Wℓ+ 1 ,j+ 1 δy Uℓ+ 1 ,j+ 1 = − 2

2

2

I −1 

Wℓ+ 1 ,j+ 1 δx Vℓ+ 1 ,j+ 1 = − 2

2

2

Uℓ+ 1 ,j δy Wℓ+ 1 ,j , 2

j =1

2

ℓ=1

2

Vℓ,j+ 1 δx Wℓ,j+ 1 . 2

Proof. This is summation by parts. See [19].

2



To prove unconditional stability of the operator splitting schemes constructed in this paper we prove a discrete equivalent of the (continuous) energy decay (2.11) with respect to the discrete energy Ehn defined as

 12   2       1 √ √ 2 2   Ehn =  ϵ0 ϵ∞ En E +   Pn  +  µ0 H n H  ,  ϵ0 ϵ∞ (ϵq − 1) 

(4.3)

E

where 0 ≤ n ≤ N − 1. The discrete grid norms ∥ · ∥E , and ∥ · ∥H are defined in (3.12), and (3.13) respectively. We now prove that the SQ-DB-D scheme is unconditionally stable. Theorem 4.1 (Discrete Energy Decay for Scheme SQ-DB-D). For all 0 ≤ n ≤ N − 1, and for all 1t , 1x, 1y > 0, the sequential operator splitting scheme SQ-DB-D (3.25)–(3.27) satisfies the energy decay

Ehn ≤ Eh0 ,

(4.4)

for the discrete energy Ehn defined in (4.3). Since the discrete energy decay (4.4) is satisfied unconditionally, the sequential operator splitting scheme SQ-DB-D (3.25)–(3.27) is unconditionally stable.

˜ n+11 Proof. In Step A1 , we multiply (3.25a) by µ0 1t (H

ℓ+ 2 ,j+ 12

µ0



˜ n+11 H

2

ℓ+ 2 ,j+ 12

 2 n − Hℓ+ 1 ,j+ 1 2



1t

=

Multiplying (3.25b) by ϵ0 ϵ∞ 1t (E˜ xn+11 + Exn

1 ,j ℓ+ 2

ϵ0 ϵ∞



E˜ xn+1

2

 −

ℓ+ 1 ,j 2

Exn 1 ℓ+ 2 ,j

2  =

1t

2

˜ n+11 H

2

n + Hℓ+ 1 ,j + 1

ℓ+ 2 ,j+ 12

2

2

ℓ+ 2 ,j



n + Hℓ+ ) to obtain 1 ,j + 1

2

2



 δy E˜ xn+11

1 ℓ+ 2 ,j+ 2

+ Exn

ℓ+ 1 ,j+ 1 2 2



.

(4.5)

) gives us



2



E˜ xn+11 + Exn

ℓ+ 1 ,j 2

ℓ+ 2 ,j

  n +1 δy H˜ ℓ+ + Hn 1 . 1 ℓ+ ,j ,j

(4.6)

2

2

Adding together (4.5) and (4.6) and summing over all spatial nodes, we obtain J −1 L−1  

 ϵ0 ϵ∞



ℓ=0 j=0

=

Exn+11 ℓ+ 2 ,j

˜

2

 −

Exn 1 ℓ+ 2 ,j

J −1  L −1  1t 

2 ℓ=0 j=0

 +

˜ n+11 H

ℓ+ 2 ,j+ 12

E˜ xn+11 + Exn

1 ,j ℓ+ 2

ℓ+ 2 ,j

n + Hℓ+ 1 ,j+ 1 2



2



2 

+ µ0



˜ H

2

n +1

ℓ+ 12 ,j+ 12



n

− Hℓ+ 1 ,j+ 1 2

2



2

  n +1 δy H˜ ℓ+ + Hn 1 1 ℓ+ ,j ,j 2

2



δy E˜ xn+11

1 ℓ+ 2 ,j+ 2

+ Exn

1 ℓ+ 1 ,j+ 2 2



.

(4.7)

V.A. Bokil et al. / Journal of Computational and Applied Mathematics 263 (2014) 160–188

175

By Lemma 4.1, the right-hand side vanishes and we see that

ϵ0 ϵ∞

 J −1 L−1  

2

E˜ xn+1

−

1 ,j ℓ+ 2

ℓ=0 j=0



Exn 1 ℓ+ 2 ,j

 J −1 L−1  

2  = −µ0

˜ n+11 H

2

2



2



.

(4.8)

.

(4.9)

, and H n+11

n +1 = Hˆ ℓ+ , respectively, 1 ,j + 1

n

− Hℓ+ 1 ,j

ℓ+ 2 ,j

ℓ=0 j=0



2

From (3.26c), using Eˆ xn+11 ,j = E˜ xn+11 ,j we arrive at the identity ℓ+ 2

ℓ+ 2

ϵ0 ϵ∞

 J −1 L−1  

Exn+11 ℓ+ 2 ,j

ℓ=0 j=0

ˆ

2

 −

Exn 1 ℓ+ 2 ,j

 J −1 L−1  

2  = −µ0

˜ H

ℓ=0 j=0

n +1

2



n

− Hℓ+ 1 ,j

ℓ+ 12 ,j

2

We repeat this process for Step A2 . From (3.25c) and (3.27e), using E˜ yn+1 1 = Eyn we obtain the identity

ϵ0 ϵ∞

 J −1 L−1  

Eyn+1 1 ℓ,j+ 2

ℓ=0 j=0

ℓ,j+ 1 2

ℓ,j+ 2

ˆ

2



Eyn 1 ℓ,j+ 2

−

2 

 J −1 L−1  

= −µ0

H

2

n +1

ℓ+ 12 ,j+ 12

ℓ=0 j=0

ℓ+ 2 ,j+ 12



2 

n+1

− H˜ ℓ+ 1 ,j+ 1 2

2

.

2

(4.10)

2

Adding together (4.9) and (4.10), multiplying both sides by 1x1y, and using the definition of the discrete grid norms (3.12), and (3.13) we have

   2 √ √  2 √ 2  √ n +1  n +1  2  ˆ ϵ ϵ E + µ H −  ϵ 0 ϵ ∞ E n  E +  µ0 H n  H = 0 .  0 ∞  0 H

(4.11)

E

Next, we multiply (3.27a) by ϵ0 ϵ∞ 1t (Exn+11 + Eˆ xn+11 ) and (3.27b) by ϵ ϵ 1(ϵt −1) (Pxn+11 + Pˆ xn+11 ) to obtain 0 ∞ q ℓ+ ,j ℓ+ ,j ℓ+ ,j ℓ+ ,j 2

ϵ0 ϵ∞



Exn+11 ℓ+ 2 ,j

2



Exn+11 ℓ+ 2 ,j

− ˆ

2

2



2  = −

21t ϵ0 ϵ∞ (ϵq − 1) 

τ

 +



1

ϵ0 ϵ∞ (ϵq − 1)

Pxn+11

2

ℓ+ 2 ,j



− Pˆxn+11

21t 

τ

ℓ+ 2 ,j

21t 

τ



ℓ+ 2 ,j

Pxn+11 + Pˆ xn+11 ℓ+ 2 ,j

ℓ+ 2 ,j

21t   τ ϵ0 ϵ∞ (ϵq − 1)



ℓ+ 2 ,j

2



ℓ+ 2 ,j

2



 

Exn+11 + Eˆ xn+11

 

 −

ℓ+ 2 ,j

2

2

 =

ℓ+ 2 ,j

ℓ+ 2 ,j

Pxn+11 + Pˆ xn+11



2 

Exn+11 + Eˆ xn+11



2

2

 ,

Exn+11 + Eˆ xn+11

 

ℓ+ 2 ,j

2

Pxn+11 + Pˆ xn+11 ℓ+ 2 ,j

ℓ+ 2 ,j

ℓ+ 2 ,j

2

(4.12)

  

2   .

(4.13)

Adding (4.12) and (4.13) together and summing over all spatial nodes we have

 J −1  L−1    ϵ0 ϵ∞  ℓ=0 j=0

 

Exn+11

ℓ+ 2 ,j

2

 − Eˆ xn+11

2  +

ℓ+ 2 ,j

Pxn+11 ℓ+ 2 ,j

2



− Pˆxn+11

ℓ+ 2 ,j

ϵ0 ϵ∞ (ϵq − 1)

2   21t  X1 , =−  τ ϵ0 ϵ∞ (ϵq − 1)

(4.14)

where X1 is defined as

    2 Exn+11 + Eˆ xn+11 Pxn+11 + Pˆ xn+11 J −1 L −1   ℓ+ 2 ,j  ℓ+ 2 ,j    ℓ+ 2 ,j  ℓ+ 2 ,j X1 := ϵ0 ϵ∞ (ϵq − 1)  −  . ℓ=0 j=0

2

2

(4.15)

176

V.A. Bokil et al. / Journal of Computational and Applied Mathematics 263 (2014) 160–188

Repeating this process, from Eqs. (3.27c) and (3.27d) we obtain





J −1  L−1    ϵ0 ϵ∞  ℓ=0 j=0



Eyn+1 1

2

ℓ,j+ 2



− Eˆ yn+1 1

2  +

ℓ,j+ 2

Pyn+1 1 ℓ,j+ 2

2



− Pˆyn+1 1

2 

ℓ,j+ 2

ϵ0 ϵ∞ (ϵq − 1)

 21t  X2 , =−  τ ϵ0 ϵ∞ (ϵq − 1)

(4.16)

where X2 is defined as

    2 Eyn+1 1 + Eˆ yn+1 1 Pyn+1 1 + Pˆ yn+1 1 J −1 L−1   ℓ,j+ 2  ℓ,j+ 2   ℓ,j+ 2   ℓ,j+ 2 X2 := ϵ0 ϵ∞ (ϵq − 1)  −  . 2

ℓ=0 j=0

Adding (4.14) and (4.16), using Pˆ xn+11 = Pxn

ℓ+ 1 ,j 2

ℓ+ 2 ,j

(4.17)

2

from (3.25d) and (3.26d) and Pˆ yn+1 1 = Pyn ℓ,j+ 2

multiplying both sides by 1x1y we obtain the identity

1 ℓ,j+ 2

from (3.25e) and (3.26e), and



  2   2    √ 2   √  1 1  n +1   n  ˆ n +1   ϵ0 ϵ∞ En+1 2 +     P − ϵ ϵ E + P       0 ∞ E  ϵ0 ϵ∞ (ϵq − 1)  ϵ0 ϵ∞ (ϵq − 1)   E E

=−

E

21t 1x1y (X1 + X2 ). τ ϵ0 ϵ∞ (ϵq − 1)

(4.18)

Adding (4.11) and (4.18), we have the identity

Ehn+1



2

 2 = Ehn −

21t 1x1y

τ ϵ0 ϵ∞ (ϵq − 1)

(X1 + X2 ),

(4.19)

where the discrete energy Ehn is defined in (4.3). Since X1 ≥ 0 and X2 ≥ 0, we have that

Ehn+1 ≤ Ehn .

(4.20)

Since we have showed this for arbitrary n, the result holds for all n, i.e. we arrive at the energy decay given in (4.4). Since this energy decay does not depend on any stability condition, the energy decay holds for all 0 ≤ n ≤ N − 1 and for all 1t , 1x, 1y > 0. Thus, the SQ-DB-D scheme is unconditionally stable.  Remark 4.1. We can rewrite Eq. (4.19) as



Ehn+1 − Ehn



Ehn+1 + Ehn

1t

2

 =−

1x1y (X1 + X2 ). τ ϵ0 ϵ∞ (ϵq − 1)

(4.21)

From the definitions of X1 and X2 in (4.15) and (4.17), respectively, and the definitions of the grid norms (3.12), and (3.13), the term 1x1y(X1 + X2 ) can be written as the norm

    2  ˆ n+1  En+1 + Eˆ n+1 Pn+1 + P   1x1y(X1 + X2 ) = ϵ0 ϵ∞ (ϵq − 1) −  .   2 2

(4.22)

E

Thus, from Eqs. (4.21) and (4.22) and the definitions of the temporal difference operator in (3.8) and the averaging operator in (3.11), we can rewrite Eq. (4.19) as

 n+ 21

δt Eh

=

−1 n+ 12

Eh

(τ ϵ0 ϵ∞ (ϵq − 1))

    2  ˆ n +1  En+1 + Eˆ n+1 Pn+1 + P    ϵ0 ϵ∞ (ϵq − 1) −  .   2 2

(4.23)

E

For the continuous case the time derivative of the continuous energy (2.11) is dE (t ) dt

 =

 −1 ∥ϵ0 ϵ∞ (ϵq − 1)E − P∥2 ≤ 0. E (t )(τ ϵ0 ϵ∞ (ϵq − 1))

(4.24)

We observe that Eq. (4.23) is a discrete approximation to the time derivative of the (continuous) energy given in (4.24). In our numerical simulations we will illustrate that this approximation is first order accurate. We show using a similar approach that the SY-DB-D scheme is also unconditionally stable.

V.A. Bokil et al. / Journal of Computational and Applied Mathematics 263 (2014) 160–188

177

Theorem 4.2 (Discrete Energy Decay for Scheme SY-DB-D). For all 0 ≤ n ≤ N − 1, and for all 1t , 1x, 1y > 0, the symmetrized operator splitting scheme SY-DB-D (3.45)–(3.49) satisfies the energy decay

Ehn ≤ Eh0 ,

(4.25)

for the discrete energy Ehn defined in (4.3). Since the discrete energy decay (4.25) is satisfied unconditionally, the symmetrized operator splitting scheme SY-DB-D (3.45)–(3.49) is unconditionally stable. Proof. We define X1 and X2 as





L −1 J −1

X1 :=

  ϵ0 ϵ∞ (ϵq − 1)   

n+ 21

E˜ xn+11 + Ex ℓ+ 2 ,j

2

ℓ=0 j=0

 X2 :=



J −1 L −1     ϵ0 ϵ∞ (ϵq − 1)   

,j ℓ+ 1 2

n+ 21

E˜ yn+1 1 + Ey ℓ,j+ 2

ℓ,j+ 1 2

2

ℓ=0 j=0





  −   



  −  

n+ 12

P˜ xn+11 + Px ℓ+ 2 ,j

,j ℓ+ 1 2

2 n+ 12

P˜ yn+1 1 + Py

1 ℓ,j+ 2

ℓ,j+ 2

2

2   , 

(4.26a)

2   . 

Using the same process from the proof of Theorem 4.1, we obtain from the first two steps, identity

(4.26b)

Step 12 A1 and Step 12 A2 the

 √    √ 2 √ 2  1 2 1 2  √  ϵ0 ϵ∞ En+ 2  +  µ0 Hn+ 2  −  ϵ0 ϵ∞ En E +  µ0 Hn H = 0. E

(4.27)

H

From Step B we obtain the identity



  2   2  √ √   2    1 1 1 2        n + 1 n + 1 n + n  ϵ0 ϵ∞ E˜  +   P   −  ϵ0 ϵ∞ E 2  +   P    ϵ0 ϵ∞ (ϵq − 1)  ϵ0 ϵ∞ (ϵq − 1)   E E E

=−

E

21t 1x1y (X1 + X2 ). τ ϵ0 ϵ∞ (ϵq − 1)

(4.28)

Finally, from the last two steps, Step 12 A2 and Step 21 A1 we obtain the identity

 2 √ 2  √     √  n+1  n+ 12   ϵ0 ϵ∞ En+1 2 + √µ0 Hn+1 2 −  ˜  ϵ0 ϵ∞ E  +  µ0 H  = 0. E

H

E

(4.29)

H

Adding (4.27)–(4.29), we have the following energy decay for the SY-DB-D scheme



Ehn+1

2

 2 = Ehn −

21t 1x1y (X1 + X2 ), τ ϵ0 ϵ∞ (ϵq − 1)

(4.30)

where the discrete energy Ehn is defined in (4.3). Using the fact that X1 ≥ 0 and X2 ≥ 0, we have that

Ehn+1 ≤ Ehn .

(4.31)

Since we have showed this for arbitrary n, the result holds for all n, i.e. we arrive at the energy decay given in (4.25), which shows that the SY-DB-D scheme is unconditionally stable.  Remark 4.2. Using the definitions of the terms X1 , X2 from (4.26), we can rewrite Eq. (4.30) as

 n+ 21

δt Eh

=

    1  ˜ n+1 + Pn+ 12 E˜ n+1 + En+ 2 P   ϵ0 ϵ∞ (ϵq − 1) −    2 2 τ ϵ0 ϵ∞ (ϵq − 1) −1

n+ 21

Eh

2    . 

(4.32)

E

In our numerical simulations we will illustrate that this approximation to the time derivative of the (continuous) energy (4.24) is second order accurate for the SY-DB-D scheme.

178

V.A. Bokil et al. / Journal of Computational and Applied Mathematics 263 (2014) 160–188

5. The Crank–Nicolson scheme for the Maxwell–Debye model To compute the order of accuracy in space/time of the splitting schemes, we will show in Section 6 that the splitting schemes are time perturbations of the Crank–Nicolson (CN) scheme. Like the Yee scheme, the CN scheme also staggers the electric and magnetic fields in space. However it does not stagger the fields in time. Instead, like the operator splitting schemes constructed in this paper, the CN scheme uses an implicit time discretization based on the time averaging operator defined in (3.11). For the 2D Maxwell–Debye system (2.12), the CN scheme is as follows: (CN): 1

n+ 21

δt Hℓ+ 1 ,j+ 1 = 2

n+ 21

δt Ex

,j ℓ+ 1 2

n+ 21

δt Ey

ℓ,j+ 1 2

n+ 12 t Px 1 ℓ+ 2 ,j

δ

n+ 21

δt Py

ℓ,j+ 1 2

µ0

2

=

  n+ 1 δy E x 21



1 ℓ+ 2 ,j+ 2

− δx



n+ 1 E y 21 1 ℓ+ 2 ,j+ 2



,

(5.1a)

      (ϵq − 1) n+ 21 1 n+ 1 n+ 1 + , δy H ℓ+ 12 ,j − Ex 1 P x 21 ℓ+ 2 ,j ℓ+ 2 ,j ϵ0 ϵ∞ τ τ ϵ0 ϵ∞ 2 1

=−

(5.1b)

      (ϵq − 1) n+ 21 1 n+ 1 n+ 1 Ey 1 + Py 2 1 , δx H ℓ,j+2 1 − ℓ,j+ 2 ℓ,j+ 2 ϵ0 ϵ∞ τ τ ϵ0 ϵ∞ 2 1

(5.1c)

    ϵ0 ϵ∞ (ϵq − 1) n+ 21 1 n+ 12 = − , Ex 1 Px 1 ℓ+ 2 ,j ℓ+ 2 ,j τ τ

(5.1d)

    1 ϵ0 ϵ∞ (ϵq − 1) n+ 12 n+ 1 Ey 1 − Py 2 1 . ℓ,j+ 2 ℓ,j+ 2 τ τ

(5.1e)

=

We next compute the local truncation errors (LTE) of the CN scheme and show that they are second order accurate. The n n LTE is defined by replacing the finite difference solution Uℓ, j with the exact solution U (xℓ , yj , t ) for U = H , Ex , Ey , Px , Py in the finite difference scheme (5.1). Since, the exact solution does not satisfy the finite difference scheme exactly, the discrepancy is the LTE which is denoted by ξU , U = H , Ex , Ey , Px , Py for the CN scheme (5.1). Lemma 5.1 (Crank–Nicolson Truncation Errors). Suppose that the solutions to the two dimensional Maxwell–Debye   ¯ )]2 , P ∈ system (2.12) (or equivalently (2.5)) are smooth enough and satisfy the regularity conditions E ∈ C 3 [0, T ]; [C 3 (Ω n+ 12

¯ )]2 and H ∈ C 3 [0, T ]; C 3 (Ω ¯ ) . Let ξH C 3 [0, T ]; [C (Ω Crank–Nicolson scheme (5.1). Then 







n+ 21

, ξEx

n+ 12

, ξEy

n+ 21

, ξPx

           n+ 21   n+ 12   n+ 21   n+ 21   n+ 21             max ξH  , ξEx  , ξEy  , ξPx  , ξPy  ≤ C 1x2 + 1y2 + 1t 2 ,   where C = C ϵ0 , µ0 , ϵ∞ , ϵq , τ does not depend on the mesh sizes 1x, 1y, and 1t.

n+ 12

, ξPy

be the truncation errors for the

(5.2)

Proof. The local truncation errors for the CN scheme (5.1) have the form n+ 1

(ξH )ℓ+ 12 ,j+ 1 = 2

2

   1   H xℓ+ 1 , yj+ 1 , t n+1 − H xℓ+ 1 , yj+ 1 , t n 2 2 2 2 1t −

+

1  2µ0

     δy Ex xℓ+ 1 , yj+ 1 , t n+1 + Ex xℓ+ 1 , yj+ 1 , t n 2

2

2

2

1  2µ0

     δx Ey xℓ+ 1 , yj+ 1 , t n+1 + Ey xℓ+ 1 , yj+ 1 , t n , 2

2

2

(5.3a)

2

    n+ 12 1   ξEx ℓ+ 1 ,j = Ex xℓ+ 1 , yj , t n+1 − Ex xℓ+ 1 , yj , t n 2 2 1t 2       1 − δy H xℓ+ 1 , yj , t n+1 + H xℓ+ 1 , yj , t n 2 2 2ϵ0 ϵ∞      (ϵq − 1) + Ex xℓ+ 1 , yj , t n+1 + Ex xℓ+ 1 , yj , t n 2 2 2τ      1 − Px xℓ+ 1 , yj , t n+1 + Px xℓ+ 1 , yj , t n , 2 2 2τ ϵ0 ϵ∞

(5.3b)

V.A. Bokil et al. / Journal of Computational and Applied Mathematics 263 (2014) 160–188







ξEy

ξPx

ξPy

n+ 21 ℓ,j+ 21

n+ 12 ℓ+ 12 ,j

n+ 21 ℓ,j+ 21

179

   1   Ey xℓ , yj+ 1 , t n+1 − Ey xℓ , yj+ 1 , t n = 2 2 1t       1 δx H xℓ , yj+ 1 , t n+1 + H xℓ , yj+ 1 , t n+1 + 2 2 2ϵ0 ϵ∞    (ϵq − 1)   + Ey xℓ , yj+ 1 , t n+1 + Ey xℓ , yj+ 1 , t n 2 2 2τ      1 Py xℓ , yj+ 1 , t n+1 + Py xℓ , yj+ 1 , t n , − 2 2 2τ ϵ0 ϵ∞    1   = Px xℓ+ 1 , yj , t n+1 − Px xℓ+ 1 , yj , t n 2 2 1t      ϵ0 ϵ∞ (ϵq − 1) − Ex xℓ+ 1 , yj , t n+1 + Ex xℓ+ 1 , yj , t n 2 2 2τ      1 Px xℓ+ 1 , yj , t n+1 + Px xℓ+ 1 , yj , t n , + 2 2 2τ      1 = Py xℓ , yj+ 1 , t n+1 − Py xℓ , yj+ 1 , t n 2 2 1t    ϵ0 ϵ∞ (ϵq − 1)   − Ey xℓ , yj+ 1 , t n+1 + Ey xℓ , yj+ 1 , t n 2 2 2τ      1 Py xℓ , yj+ 1 , t n+1 + Py xℓ , yj+ 1 , t n . + 2 2 2τ

(5.3c)

(5.3d)

(5.3e)

By performing Taylor expansions around the appropriate points and using (2.12), we arrive at the local truncation errors for the Crank–Nicolson scheme: n+ 1

  ∂ 3 Ex  1 ∂ 3 Ey  , y , t , t + x x , y 1 11 12 13 11 j+ 12 2 2 24 ∂ t 3 8µ0 ∂ t 2 ∂ y ℓ+ 2 8µ0 ∂ t 2 ∂ x 2       3 3 1 ∂ Ey ∂ Ex 1 1 + 1x2 3 x12 , yj+ 1 , t n+ 2 − 1y2 3 xℓ+ 1 , y21 , t n+ 2 , 2 2 24µ0 ∂x ∂y    (ϵ − 1) ∂ 2 E   1 1 ∂ 3 Ex  ∂ 3H  q x − + = 1t 2 x x x 1 , yj , t21 1 , y21 , t22 1 , yj , t23 ℓ+ 2 ℓ+ 2 ℓ+ 2 2 24 ∂ t 3 8τ ∂t2 8ϵ0 ϵ∞ ∂ t2 ∂ y 3     2 1 ∂ Px 1y ∂ H 1 − xℓ+ 1 , yj , t24 − x21 , yj , t n+ 2 , 2 8τ ϵ0 ϵ∞ ∂ t 2 24ϵ0 ϵ∞ ∂ y3    (ϵ − 1) ∂ 2 E   ∂ 3H  1 ∂ 3 Ey  1 q y 2 = 1t x , y + x , y + x , y 1 , t31 1 , t32 1 , t33 ℓ 31 ℓ j + j + j + 2 2 2 24 ∂ t 3 8ϵ0 ϵ∞ ∂ t 2 ∂ x 8τ ∂t2   ∂ 2 Py  1x2 ∂ 3 H  1 1 − xℓ , yj+ 1 , t34 + x31 , yj+ 1 , t n+ 2 , 2 2 8τ ϵ0 ϵ∞ ∂ t 2 24ϵ0 ϵ∞ ∂ x3   ϵ ϵ (ϵ − 1) ∂ 2 E    1 ∂ 3 Px  1 ∂ 2 Px  0 ∞ q x = 1t 2 x − x + x , 1 , yj , t41 1 , yj , t42 1 , yj , t43 ℓ+ 2 ℓ+ 2 ℓ+ 2 24 ∂ t 3 8τ ∂t2 8τ ∂ t 2   ϵ ϵ (ϵ − 1) ∂ 2 E    1 ∂ 3 Py  1 ∂ 2 Py  0 ∞ q y 2 = 1t xℓ , yj+ 1 , t51 − xℓ , yj+ 1 , t52 + xℓ , yj+ 1 , t53 , 2 2 2 24 ∂ t 3 8τ ∂t2 8τ ∂ t 2

(ξH )ℓ+ 12 ,j+ 1 = 1t 2 2

n+ 1

(ξEx )ℓ+ 12 ,j 2

n+ 1

(ξEy )ℓ,j+2 1 2

n+ 1

(ξPx )ℓ+ 12 ,j 2

n+ 1

(ξPy )ℓ,j+2 1 2



1 ∂ 3H 



xℓ+ 1 , yj+ 1 , t11 −

1

where xℓ− 1 ≤ xi1 , xi2 ≤ xℓ+ 1 , yj− 1 ≤ yi1 , yi2 ≤ yj+ 1 , and t n ≤ tij ≤ t n+1 for i ∈ {1, 2, 3, 4, 5}, and j ∈ {1, 2, 3}. 2

2

2

2

(5.4a)

(5.4b)

(5.4c)

(5.4d)

(5.4e)



6. Equivalent operator splitting schemes In this section we analyze the accuracy of the operator splitting schemes SQ-DB-D and SY-DB-D. To this end, we construct equivalent schemes for the sequential and symmetrized splitting methods that eliminate the intermediate variables and present the splitting schemes in un-split form. We show that the equivalent sequential operator splitting scheme, EQ-SQ-DB-D, is a first order time perturbation of the CN scheme, whereas the equivalent symmetrized operator splitting scheme, EQ-SY-DB-D, is a second order time perturbation of the CN scheme.

180

V.A. Bokil et al. / Journal of Computational and Applied Mathematics 263 (2014) 160–188

˜ and Uˆ produces the following equivalent sequential Rearranging Eqs. (3.25)–(3.27) to eliminate intermediate variables U operator splitting scheme. (EQ-SQ-DB-D):



n+ 12

δt Hℓ+ 1 ,j+ 1 = 2

n+ 21

δt Ex

,j ℓ+ 1 2

n+ 21

δt Ey

ℓ,j+ 1 2

n+ 21

δt Px

ℓ+ 1 ,j 2

n+ 21

δt Py

ℓ,j+ 1 2

2

1

µ0

  n+ 1 δy E x 21



ℓ+ 2 ,j+ 1 2

− δx



n+ 1 E y 21 1 ℓ+ 2 ,j+ 2

 +

2 c∞ 1t

2



 n+ 1 δy δt Px 21

ℓ+ 2 ,j+ 1 2



 n+ 1 − δx δt Py 21



ℓ+ 2 ,j+ 1 2

      (ϵq − 1) n+ 21 1 n+ 1 n+ 1 δy H ℓ+ 12 ,j − + Ex 1 P x 21 ℓ+ 2 ,j ϵ0 ϵ∞ τ ϵ0 ϵ∞  ℓ+ 2 ,j  2  τ1  2 2 1tc∞ 1t 2 (ϵq − 1)c∞ n+ 2 − δy δx P y 1 − δy2 Exn 1 ℓ+ 2 ,j ℓ+ 2 ,j 4ϵ0 ϵ∞ (ϵq − 1) 8τ       2 τ 1tc∞ 2τ + 1t (ϵq − 1) 1t 2 n+ 21 n+ 12 − δy δt Px 1 − δy δx δt Py 1 ℓ+ 2 ,j ℓ+ ,j 2ϵ0 ϵ∞ 4τ 2 (ϵq − 1)    2    2 2 2 1t c∞ 2 n+ 21 1t (ϵq − 1) 1 tc ∞ n − + δy P x 1 − δy Hℓ+ δy δx Eyn 1 , 1 ,j ℓ+ 2 ,j ℓ+ 2 ,j 8τ ϵ0 ϵ∞ 2τ ϵ0 ϵ∞ 4 2        1 −1 (ϵq − 1) n+ 21 n+ 1 n+ 1 = δx H ℓ,j+2 1 − Ey 1 + Py 2 1 ℓ,j+ 2 ϵ0 ϵ∞ τ ϵ0 ϵ∞  ℓ,j+ 2  2  τ1  2 2 1tc∞ 1t 2 (ϵq − 1)c∞ n+ 2 − δx δy P x 1 + δx2 Eyn 1 ℓ,j+ 2 ℓ,j+ 2 4ϵ0 ϵ∞ (ϵq − 1) 8τ       2 2 1t 2 c∞ 1t (ϵq − 1) 1tc∞ n+ 12 n 2 δx δy Ex 1 + δx P y 1 + δx Hℓ,n+j+1 1 − ℓ,j+ 2 ℓ,j+ 2 4 8τ ϵ0 ϵ∞    2τ ϵ0 ϵ∞  2 1   2 1tc∞ 2τ + 1t (ϵq − 1) 1t 2 τ n+ n+ 12 + δx δt Py 1 − δx δy δt Px 2 1 , ℓ,j+ 2 ℓ,j+ 2 2ϵ0 ϵ∞ 4τ 2 (ϵq − 1)         1t 2 c 2 ϵ0 ϵ∞ (ϵq − 1) n+ 12 1 1t (ϵq − 1)  n n+ 1 n+ 1 ∞ 2 = Ex 1 − P x 21 + δy Hℓ+ 1 ,j + δy P x 21 ℓ+ 2 ,j ℓ+ ,j ℓ+ 2 ,j τ 2τ  τ 2 22   2  8τ 1t c∞ 2τ + 1t (ϵq − 1) 2 1t 2 (ϵq − 1) 2 n n+ 12 + δy Ex 1 + δy δt Px 1 , ℓ+ 2 ,j ℓ+ 2 ,j 8τ µ0 4 4τ          2 ϵ0 ϵ∞ (ϵq − 1) n+ 12 1 1t (ϵq − 1) 1t 2 c∞ n+ 21 n+ 21 2 n = Ey 1 − Py 1 − δx Hℓ,j+ 1 − δx P y 1 ℓ,j+ ℓ,j+ ℓ,j+ 2 τ 2τ   2  τ 2 2 2  2 1  8τ 2 1t (ϵq − 1) 2 n 1t c∞ 2τ + 1t (ϵq − 1) 2 n+ 2 − δ x Ey 1 − δx δt Py 1 . ℓ,j+ 2 ℓ,j+ 2 8τ µ0 4 4τ 

=

,

(6.1a)

1

(6.1b)

(6.1c)

(6.1d)

(6.1e)

We have the following lemma. Lemma 6.1 (EQ-SQ-DB-D Truncation Errors). Suppose that the solutions to the two dimensional  Maxwell–Debye  sys¯ )]2 , P ∈ tem (2.12) (or equivalently (2.5)) are smooth enough and satisfy the regularity conditions E ∈ C 3 [0, T ]; [C 3 (Ω n+ 1

n+ 1

n+ 1

n+ 1

n+ 12

¯ )]2 , and H ∈ C 3 [0, T ]; C 3 (Ω ¯ ) . Let ξˆH 2 , ξˆE 2 , ξˆE 2 , ξˆP 2 , ξˆP C 3 [0, T ]; [C (Ω x y x y EQ-SQ-DB-D scheme. Then, for sufficiently small mesh step sizes, 1x, 1y and 1t we have 



 

1

   



1

   

1

   



1

   

1

 

n+ n+ n+ n+ n+ max ξˆH 2  , ξˆEx 2  , ξˆEy 2  , ξˆPx 2  , ξˆPy 2 

  ≤ C1 1x2 + 1y2 + 1t ,

be the truncation errors for the

(6.2)

where C1 = C1 ϵ0 , µ0 , ϵ∞ , ϵq , τ does not depend on the mesh sizes 1x, 1y, and 1t.





Proof. The terms in square brackets on the right hand side in each sub-equation of (6.1) form the right hand side of the corresponding equation of the CN scheme (5.1). The remaining terms have at least a factor of 1t in them. We observe that we can write the LTE for the EQ-SQ-DB-D as first order in time perturbations of the local truncation errors for the CN scheme (5.1). Thus,

 n+ 12 n+ 1 ξˆV = (ξV )α,β2 + O (1t ), α,β

(6.3)

where the variable ξˆ is used to denote truncation errors for the EQ-SQ-DB-D scheme, whereas the variable ξ is used to 1

denote truncation errors for the CN scheme. The letter V denotes one of the variables H , Ex , Ey , Px , Py , and (xα , yβ , t n+ 2 ) is the space–time point around which the Taylor expansions are taken to obtain the required truncation errors. 

V.A. Bokil et al. / Journal of Computational and Applied Mathematics 263 (2014) 160–188

181

We can similarly construct the equivalent symmetrized operator splitting scheme (EQ-SY-DB-D) by eliminating intermediate variables from Eqs. (3.45)–(3.49). The equivalent scheme is as follows (EQ-SY-DB-D): n+ 12

δt Hℓ+ 1 ,j+ 1 2

2

          2 1t 2 c∞ 1 2 1 2 n+ 21 n+ 21 n+ 12 n+ 12 δy E − δx E y 1 1 + δx δy E x 1 1 + δ δx E y 1 1 = ℓ+ 2 ,j+ 2 ℓ+ 2 ,j+ 2 ℓ+ 2 ,j+ 2 µ0  xℓ+ 12 ,j+12 16 µ µ0 y 0    2 2 1 1 1 1 t c n + n + n + ∞ 2 2 − δx2 δt Hℓ+ 12 ,j+ 1 − δy2 δt Hℓ+ 12 ,j+ 1 − δx δy δt Hℓ+ 12 ,j+ 1 , (6.4a) 

1

2

δt E

n+ 12 x 1 ℓ+ 2 ,j

n+ 12

δ t Ey

ℓ,j+ 1 2

n+ 21

δ t Px

ℓ+ 1 ,j 2

n+ 12

δ t Py

ℓ,j+ 1 2

2

2

16

2

2

2

1 (ϵq − 1) 1 − E + P δy H ϵ0 ϵ∞  τ τ ϵ 0 ϵ∞       2 2 1t 2 c∞ 1t 2 c∞ n+ 1 n+ 1 + δx2 δt Ex 21 − δy2 δt E + δx2 δy2 δt Ex 21 ℓ+ 2 ,j ℓ+ ,j 16 16   2     2 2 1 1 1 1 t c 1 n + n + n + ∞ 2 2 + δ2 P 2 − − δ2 P 2 δx δy P x 21 ℓ+ ,j τ ϵ0 ϵ∞ x xℓ+ 12 ,j  y xℓ+ 12 ,j  16  2 1  2 2 1 (ϵq − 1) 2 n+ 12 1 t c n + n+ ∞ 2 2 + δx E x 1 − δy2 E x 21 + δx δy E x 21 ℓ+ 2 ,j ℓ+ 2 ,j ℓ+ ,j τ 16   2      2 1 4(ϵq − 1)µ0 1t 2 c∞ 1 n+ 12 n + n+ 1 + − , δy δt Hℓ+ 1 ,j − δx2 δy δt Hℓ+ 12 ,j δx2 δy H ℓ+ 12 ,j τ 16 ϵ0 ϵ∞ 2 2 2        −1 (ϵq − 1) n+ 21 1 n+ 1 n+ 1 Ey 1 + Py 2 1 = δx H ℓ,j+2 1 − ℓ,j+ ℓ,j+ 2 ϵ0 ϵ∞  2 τ    21  τ ϵ0 ϵ∞  2 2 1t 2 c∞ 1t 2 c∞ n+ 21 n+ 2 n+ 1 2 2 2 2 δy δt Ey 1 − δx δt Ey 1 + δx δy δt Ey 2 1 + ℓ,j+ 2 ℓ,j+ 2 ℓ,j+ 16 16      2   2 2 1 1 1 1 1 t c n + n + n + ∞ 2 2 δ2 P 2 δx δy P y 2 1 − + δ2 P 2 − ℓ,j+ τ ϵ0 ϵ∞ x yℓ,j+ 21  y yℓ,j+ 12  16  21  2 2 1 (ϵq − 1) 2 n+ 12 1 t c n + n+ ∞ 2 2 + δ E − δx2 E y 2 1 + δx δy E y 2 1 1 +2 ℓ,j+ 2 ℓ,j+ τ  y yℓ,j 16    2 1  2 1 1t 2 c∞ 1 n + n+ 12 n+ δy2 δx δt Hℓ,j+21 − δx δy E x 2 1 − 4µ0 δx δt Hℓ,j+ 1 + ℓ,j+ 2 16  µ0 2  2   1 1 1 n + n + δ 2 δx H ℓ,j+2 1 + 4δx δy δt Ex 2 1 , − ℓ,j+ 2 ϵ0 ϵ∞ y 2         2 1 1t 2 c∞ ϵ0 ϵ∞ (ϵq − 1) n+ 21 n+ 1 n+ 1 − + = Ex 1 P x 21 δx2 δt Px 21 ℓ+ ,j ℓ+ 2 ,j τ ℓ+ 2 ,j  16    2   τ 1  2 2 1 1 1 t c 1 n + n + n + n+ 12 ∞ 2 2 2 2 2 2 2 2 + δy δt Px 1 − δx δy δt Px 1 + δ P x 1 + δy P x 1 ℓ+ 2 ,j ℓ+ 2 ,j ℓ+ ,j ℓ+ ,j 16 τ x     2   2 2 2 1 1t c∞ 2 2 n+ 12 ϵ0 ϵ∞ (ϵq − 1) 2 n+ 12 n + − δx δy P x 1 + δy E x 1 − δx2 E x 21 ℓ+ 2 ,j ℓ+ 2 ,j ℓ+ 2 ,j 16 τ       2 2 1 1 1 1t 2 c∞ 1 t n + n + n + − δx2 δy2 E x 21 − 4µ0 δy δt Hℓ+ 12 ,j + δx2 δy δt Hℓ+ 12 ,j , ℓ+ 2 ,j 16 4ϵ0 ϵ∞ 2 2         2 1 1t 2 c∞ ϵ0 ϵ∞ (ϵq − 1) n+ 12 n+ 1 n+ 1 = Ey 1 − Py 2 1 + δx2 δt Py 2 1 ℓ,j+ ℓ,j+ 2 τ ℓ,j+ 2  16 τ 1     2   2 1t 2 c∞ 1 n+ 2 n+ 12 n+ 21 n+ 1 2 2 2 2 + δy δt Py 1 − δx δy δt Py 1 + δx P y 1 + δy2 P y 2 1 ℓ,j+ 2 ℓ,j+ 2 ℓ,j+ ℓ,j+ 16   τ   2   2 2 1 1 1 1t 2 c∞ ϵ ϵ (ϵ − 1 ) n + n + n + 0 ∞ q − δx2 δy2 P y 2 1 + δx2 E y 2 1 − δy2 E y 2 1 ℓ,j+ 2 ℓ,j+ 2 ℓ,j+ 2 16 τ      2 2 1 1 1 1t 2 c∞ 1 t n + n + n + 2 2 2 2 2 2 − δx δy E y 1 + 4µ0 δx δt Hℓ,j+ 1 + δy δx δt Hℓ,j+ 1 . ℓ,j+ 2 16 4ϵ0 ϵ∞ 2 2





n+ 12 ℓ+ 21 ,j

=







n+ 21 x 1 ℓ+ 2 ,j n+ 12 x 1 ℓ+ 2 ,j



n+ 12 x 1 ℓ+ 2 ,j



(6.4b)

(6.4c)

(6.4d)

(6.4e)

Lemma 6.2 (EQ-SY-DB-D Truncation Errors). Suppose that the solutions  to the two dimensional  Maxwell–Debye2  sys¯ )]2 , P ∈ C 3 [0, T ]; [C (Ω ¯ )] and tem (2.12) are smooth enough and satisfy the regularity conditions E ∈ C 3 [0, T ]; [C 3 (Ω n+ 12

¯ ) . Let ξ˜H H ∈ C 3 [0, T ]; C 3 (Ω 



n+ 21

, ξ˜Ex

n+ 21

, ξ˜Ey

n+ 12

, ξ˜Px

n+ 12

, ξ˜Py

be the truncation errors for the EQ-SY-DB-D scheme. Then, for

182

V.A. Bokil et al. / Journal of Computational and Applied Mathematics 263 (2014) 160–188

sufficiently small mesh step sizes, 1x, 1y and 1t we have

 

1

   

1

   

1

   

1

   

1

 

n+ n+ n+ n+ n+ max ξ˜H 2  , ξ˜Ex 2  , ξ˜Ey 2  , ξ˜Px 2  , ξ˜Py 2 

  ≤ C2 1x2 + 1y2 + 1t 2 ,

(6.5)

where C2 = C2 ϵ0 , µ0 , ϵ∞ , ϵq , τ does not depend on the mesh sizes 1x, 1y, and 1t.





Proof. The terms in the square brackets on the right hand sides in each sub-equation in (6.4) form the right hand side of the corresponding equation of the CN scheme (5.1). We again observe that we can write the local truncation errors for the EQ-SY-DB-D as second order in time perturbations of the local truncation errors for the Crank–Nicolson scheme (5.1). Thus, we have

 n+ 12 n+ 1 = (ξV )α,β2 + O (1t 2 ), ξ˜V α,β

(6.6)

where the variable ξ˜ is used to denote truncation errors for the EQ-SY-DB-D scheme, whereas the variable ξ is used to 1

denote truncation errors for the CN scheme. The letter V denotes one of the variables H , Ex , Ey , Px , Py , and (xα , yβ , t n+ 2 ) is the point around which the Taylor expansions are taken to obtain the required truncation errors.  We note that the error constants in the truncation errors for both the equivalent operator splitting schemes EQ-SQ-DB-D, and EQ-SY-DB-D, given in Lemmas 6.1, and 6.2, respectively, depend on the stiffness parameter τ . However, the parameter τ is fixed in our model. In conclusion, we have established energy estimates and unconditional stability for the SQ-DB-D and SY-DB-D schemes and consistency for the EQ-SQ-DB-D and EQ-SY-DB-D schemes. In Section 7 we illustrate our theoretical results using numerical simulations. 7. Numerical experiments We perform numerical simulations on the domain Ω = [0, a] × [0, b]. For our simulation we assume a uniform mesh with 1x = 1y = h. We introduce an exact solution to the Maxwell–Debye system (2.12) which we use to initialize our simulations. This exact solution is

  E=

Ex Ey

 −θ ky e−θ t cos(kx π x) sin(ky π y) = , θ kx e−θ t sin(kx π x) cos(ky π y) 

|k|2 −θ t e cos(kx π x) cos(ky π y), µ0 π     Px ky α (θ , |k|) e−θ t cos(kx π x) sin(ky π y) P= = , Py −kx α (θ , |k|) e−θ t sin(kx π x) cos(ky π y)

H =

(7.1a)

(7.1b) (7.1c)

where the parameter θ is a real number. In the exact solution, the wave vector is k = π (kx , ky )T , and the corresponding 2 τ |k|2 . The wave number |k| and wave number is |k|. The function α is defined by α(θ , |k|) := ϵ0 ϵ∞ τ θ 2 − (ϵq − 1)θ + c∞ parameter θ are related by the equation

θ3 −

c 2 |k|2 ϵq 2 2 θ− ∞ θ + |k|2 c∞ = 0. τ τ

(7.2)

We note that for kx a and ky b integers, the exact solution (7.1) satisfies the perfect conductor conditions on the boundary of the domain Ω , and the electric and polarization fields are divergence free on Ω . The energy defined in (2.11) for the exact solution (7.1) can be computed to be

|k|e−θ t E (t ) = 2π



 |k|2 α(θ , |k|)2 2 ab + ϵ0 ϵ∞ θ + . µ0 ϵ0 ϵ∞ (ϵq − 1) 

(7.3)

In our simulations we use various values of the Courant number (3.28) ηx = ηy = η, and various values for the wave vector components kx = ky = k. The exact dispersion relation for Debye media [9,28] relating the wave number |k| to the angular frequency ω is

|k| =

ω c0



ϵs − iωτ ϵ∞ . 1 − iωτ

(7.4)

V.A. Bokil et al. / Journal of Computational and Applied Mathematics 263 (2014) 160–188

183

Table 1 Comparison of relative errors for the SQ-DB-D and SY-DB-D schemes, for k = 1, and η = 0.5.

1t

N

1x = 1y = h

SQ-DB-D

SY-DB-D

Relative error 25 50 100 200 400

0.04 0.02 0.01 0.005 0.0025

5.26 × 10−2 2.62 × 10−2 1.31 × 10−2 6.52 × 10−3 3.26 × 10−3

0.08 0.04 0.02 0.01 0.005

Rate

Relative error

Rate

1.00 1.00 1.00 1.00

3.07 × 10−3 7.94 × 10−4 1.98 × 10−4 4.96 × 10−5 1.24 × 10−5

1.94 2.00 2.00 2.00

Squaring both sides of the dispersion relation (7.4) we can rewrite as

(iω)3 −

c 2 |k|2 ϵq 2 |k|2 iω − ∞ (iω)2 + c∞ = 0. τ τ

(7.5)

Comparing the dispersion relation (7.5) to the relation (7.2) for real θ , we note that the exact solution (7.1) corresponds to a solution for the Maxwell–Debye system (2.12) for a purely imaginary angular frequency ω = −iθ . For the discrete solution produced we compute relative errors defined as

 relative error = max 0≤n≤N

E r (t n ) E (t n )



,

(7.6)

where the absolute error E r (t n ) is defined as

 E r (t ) = ϵ0 ϵ∞ ∥E(t n ) − En ∥2E + µ0 ∥H (t n ) − H n ∥2H +

1

n

∥P(t ) − n

ϵ0 ϵ∞ (ϵq − 1)

Pn 2E

∥

 12

,

(7.7)

and E (t ) is defined in (7.3). We also define the energy error for the discrete solutions as

   dE n+ 1  n+ 1  dt t 2 − δt Eh 2   1 energy error = max   0≤n≤N −1 dE t n+ 2  dt where the discrete energy Ehn is defined in (4.3) and the time point t

n+ 21

dE dt

     ,   

(7.8)

1

t n+ 2



is the time derivative of the exact energy (7.3) computed at

.

7.1. Numerical example 1 We use a = 1, b = 1, T = 1, parameter values µ0 = 1, ϵ0 = 1 (i.e. c0 = 1), ϵ∞ = 1, ϵq = 2 (ϵs = 2), and τ = 1. The

√

real root of Eq. (7.2) depends on the value of the wave number |k| = 2kπ . In particular, for k = 1, θ ≈ 1.0532, so that θτ ≈ 1.0532. The Yee scheme for the Maxwell–Debye model is known to be second order accurate in space and time and conditionally stable with the stability condition η < √1 (in the case of a uniform mesh 1x = 1y = h) [9,29]. In order to provide 2

comparison of the operator splitting schemes with the Yee scheme we first perform simulations with η restricted to values for which the Yee scheme is stable. Table 1 presents the relative errors (7.6) and confirms the first order accuracy of the SQ-DB-D scheme and the second order accuracy of the SY-DB-D scheme for k = 1, Courant number η = 0.5, and various values of 1t and h. The variable N ∈ N is the number of time steps performed, thus N 1t = T . We note that the largest value of 1t chosen (0.04) is such that 1t /τ (τ = 1 in this example) is O (10−2 ) or lower. This is in agreement with results obtained in [9] which indicate that to resolve all time scales in the problem we must choose 1t = O (10−2 τ ) for Debye media. Table 2 presents the energy errors (7.8) and also confirms the first order accuracy of the SQ-DB-D scheme and the second order accuracy of the SY-DB-D schemes for k = 1, Courant number η = 0.5, and various values of 1t and h. Fig. 4 is a loglog plot of the relative errors for the Yee scheme compared to √the relative errors for the SQ-DB-D and SY-DB-D schemes using η = 0.5 and different values of the wave number |k| = 2kπ . We note that, since the Yee scheme staggers the electric and magnetic components in time, the relative error for the Yee scheme is computed using a time 1

1

average of the magnetic field at t n+ 2 and t n− 2 . We refer the reader to [30] for an energy analysis of the Yee scheme for the Maxwell–Debye and other dispersive models. The plots in Fig. 4 clearly indicate the first order accuracy of the SQ-DB-D scheme and the second order accuracy of the Yee and SY-DB-D schemes. We also note that the relative errors of the SY-DB-D schemes are higher than those of the Yee scheme and this difference in errors and the errors themselves grow larger as the wave number is increased. However, the accuracy of the schemes is maintained across all wave numbers that are considered

184

V.A. Bokil et al. / Journal of Computational and Applied Mathematics 263 (2014) 160–188 Table 2 Comparison of energy errors for the SQ-DB-D and SY-DB-D schemes, for k = 1, and η = 0.5. N

1t

1x = 1y = h

SQ-DB-D Energy error

25 50 100 200 400

0.04 0.02 0.01 0.005 0.0025

0.08 0.04 0.02 0.01 0.005

5.19 × 10−2 2.64 × 10−2 1.33 × 10−2 6.67 × 10−3 3.34 × 10−3

SY-DB-D Rate

Energy error

Rate

0.98 0.99 0.99 1.00

4.00 × 10−3 1.06 × 10−3 2.65 × 10−4 6.64 × 10−5 1.66 × 10−5

1.92 1.99 2.00 2.00

√

Fig. 4. A comparison of the relative errors for different wave numbers |k| = 2kπ of the sequential splitting scheme SQ-DB-D, and the symmetrized splitting scheme SY-DB-D with the Yee scheme for example 1. Here the Courant number is fixed at η = 0.5.

here. We also observe that the growth in the errors of the SQ-DB-D scheme as the wave number is increased is negligible on the scale chosen in this figure. Fig. 5 provides similar comparisons as in Fig. 4, but with η = 0.7. Fig. 6 provides similar comparisons with η = 1.0. For this choice of η the Yee scheme is unstable. However, the splitting schemes maintain their accuracy across all wave numbers. Since, the splitting schemes are unconditionally stable, this is a distinct advantage that these schemes have over the conditionally stable Yee scheme. In particular, these simulations imply that in heterogeneous media the splitting schemes can use a large η, or in other words a refined mesh with a much larger value of 1t than is allowed in the Yee scheme, to resolve small heterogeneities. In Fig. 7 we present a loglog plot of the energy errors for the splitting schemes with k = 1, 5, 10 and different values of η including values for which the Yee scheme is unstable (η = 1). The figures confirm the first order accuracy of the SQ-DB-D scheme and the second order accuracy of the SY-DB-D scheme as discussed in Remarks 4.1 and 4.2 for the time derivative of the continuous energy and the discrete approximations to this derivative computed for the SQ-DB-D and SY-DB-D schemes, respectively. For the sequential scheme SQ-DB-D the energy errors (7.8) corresponding to different values of η lie on top of each other. As η is increased these energy errors decrease, but the decrease is not observable on the chosen scale of the plot. We observe that the energy errors for the symmetrized scheme SY-DB-D decrease as the Courant number η is increased. 7.2. Numerical example 2 In this example we use more realistic values of the parameters to illustrate convergence. We use a = b = c∞ τ , T = τ , parameter values ϵ0 = 8.85 × 10−12 , c0 = 3 × 108 , µ0 = 21 , ϵ∞ = 1, ϵq = 80.35 (ϵs = 80.35), and τ = 8 × 10−12 . The c0 ϵ 0

parameter values describe the main relaxation of water in the microwave range of frequencies [9]. We use ka = kb = 10, and we choose θ to be the real root of Eq. (7.2). For this root θ τ ≈ 1.0438. The largest value of 1t is chosen so that 1t /τ = 0.04

V.A. Bokil et al. / Journal of Computational and Applied Mathematics 263 (2014) 160–188

185

√

Fig. 5. A comparison of the relative errors for different wave numbers |k| = 2kπ of the sequential splitting scheme SQ-DB-D, and the symmetrized splitting scheme SY-DB-D with the Yee scheme for example 1. Here the Courant number is fixed at η = 0.7.

√

Fig. 6. A comparison of the relative errors for different wave numbers |k| = 2kπ of the sequential splitting scheme SQ-DB-D, and the symmetrized splitting scheme SY-DB-D for example 1. Here the Courant number is fixed at η = 1.0, for which the Yee scheme is unstable.

as was done in example 1. The spatial step size 1x = 1y = h is chosen so that the Courant number η = c∞ 1t /h = 1.0. For this CFL number the Yee scheme is unstable. In Fig. 8 we plot the relative errors (left plot) and the energy errors (right plot) of the SQ-DB-D and SY-DB-D schemes. The plots confirm the first order accuracy of the SQ-DB-D scheme and the second order accuracy of the SY-DB-D scheme. We note that the domain modeled in this example, which is of length c∞ τ , is O (10−3 ) m. As discussed in Section 3.2, in [3] it was shown for the 1D Maxwell–Debye model that this length is the size of the small boundary layer that develops

186

V.A. Bokil et al. / Journal of Computational and Applied Mathematics 263 (2014) 160–188

Fig. 7. A comparison of the energy errors of the sequential splitting scheme SQ-DB-D, and the symmetrized splitting scheme SY-DB-D for different Courant numbers η and different wave numbers for example 1.

Fig. 8. Relative (left plot) and energy (right plot) errors of the sequential splitting scheme SQ-DB-D, and the symmetrized splitting scheme SY-DB-D for example 2 with η = 1.0 and ka = kb = 10.

in the model in which the solution has a fast decaying transient term travelling with the fastest speed of harmonic waves c∞ . The small value of T = τ is the amount of time required to travel through this boundary layer. Beyond the skin depth the pulse is mainly a diffusion wave travelling with the zero-frequency phase velocity cs of harmonic waves and exhibits algebraic decay. Thus, we are simulating the evolution of the fields in very small time and space intervals and we do this mainly to test convergence of the solution in the region of fastest decay. There are important applications in which one may be interested in the response of electromagnetic pulses in the boundary layer where most of the energy of the pulse attenuates. However, beyond the skin depth important transient features of the pulse develop, such as Brillouin precursors [31], that are also important in other applications. For example, for the case of water considered here, as well as many biological tissues, the boundary layer is ≈2.4 mm. In modeling the health risk of the increased use of millimeter wave detection devices, one may want to quantify the energy absorbed by human skin, or estimate the exposure to internal organs.

V.A. Bokil et al. / Journal of Computational and Applied Mathematics 263 (2014) 160–188

187

The accuracy limitation 1t ≤ τ can be severe for long time integration. We note that this restriction is also present in the Yee scheme [9] as well as time-domain finite element methods [32] for the Maxwell–Debye model. Violating this accuracy condition will not cause instability. However, if 1t ≤ τ is not satisfied then the relevant physical scales are not fully resolved and the solution will be very inaccurate. One technique to overcome this restriction beyond the boundary layer may involve domain decomposition methods, which is a topic for our future work. 8. Conclusions In this paper, extending the ideas used in [19], we have constructed and analyzed two operator splitting schemes for Maxwell’s equations in dispersive media of Debye type (the Maxwell–Debye model). The motivation for constructing splitting schemes comes from work done in [3] which indicates that the (1D) Maxwell–Debye system for modeling wave propagation in relaxing dielectrics is stiff with at least two very different time scales. The relaxation phenomenon associated with the Debye model gives rise to a fast decaying transient in a small boundary layer of length O (c∞ τ ), and a very small time step is required in order to properly resolve all time scales [9,28]. The additional constraint of conditional stability due to the explicit nature of the Yee scheme requires an unnecessarily small time step. By removing stability constraints the time step is chosen based on accuracy requirements alone. The splitting schemes constructed in this paper, a sequential method based on the Marchuk–Yanenko scheme [18,33], and a symmetrized splitting method based on the ideas by Strang and Marchuk [18,27], are shown to be first order and second order, respectively, in time perturbations of the Crank–Nicolson method for the Maxwell–Debye model. Thus, like the Crank–Nicolson method, the splitting schemes are unconditionally stable. However, the splitting of the spatial derivative operators (wave propagation) from the lower order terms containing the relaxation time τ (polarization) gives rise to simpler numerical implementations of the Maxwell–Debye model than that of the Crank–Nicolson scheme, since we also employ dimensional splitting in addition to splitting of different physical mechanisms. By removing stability constraints the time step in the splitting methods have to be chosen to satisfy accuracy constraints. The error analysis, in Section 6, shows that the ratio 1t /τ will be an essential parameter influencing the accuracy of the splitting schemes (see local error representation), and it may be expected that theoretical error bounds will depend on this parameter, guaranteeing reasonable accuracy if this parameter is O (1) (or small). We have introduced an exact solution for the Maxwell–Debye model which we have used to provide numerical comparisons of the new operator splitting schemes for the Maxwell–Debye model with the Yee scheme. Our numerical simulations illustrate that the convergence rate of the symmetrized operator splitting scheme compares well with the convergence rate of the Yee scheme for time steps at which the Yee scheme is stable. A dispersion analysis of the operator splitting schemes presented in this paper (see [34] for a 1D analysis) and a convergence analysis of these schemes is in preparation and will be forthcoming. To prove convergence of the operator splitting schemes presented in this paper we will need to utilize the property of unconditional stability of these schemes that we proved in Section 4. In addition we will also need to compute the local truncation errors of the operator splitting schemes. Here we have only computed the truncation errors of the corresponding equivalent schemes. In [19], in which sequential and symmetrized operator splitting schemes were constructed for Maxwell’s equations in non-dispersive media (like vacuum), the authors defined new intermediate variables that transformed the equivalent splitting schemes back to the original splitting schemes. This transformation allowed for the truncation errors of the equivalent scheme to be redistributed to obtain the truncation errors of the original operator splitting schemes. Based on such a truncation error analysis and unconditional stability the energy method was used to prove convergence of the schemes. We will attempt a similar strategy for the operator splitting schemes for the dispersive Maxwell system constructed in this paper. Finally, we point out that similar splitting schemes can be constructed for other dispersive models such as the Maxwell–Lorentz model for wave propagation in media that exhibit electronic polarization [28]. Acknowledgments The authors would like to acknowledge the National Science Foundation’s Computational Mathematics program grant number DMS-0811223. The second author was supported as a Graduate Research Fellow from this grant and the third author was supported by an REU supplement to this grant. References [1] H.T. Banks, M.W. Buksas, T. Lin, Electromagnetic Material Interrogation Using Conductive Interfaces and Acoustic Wavefronts, in: Frontiers in Applied Mathematics, vol. FR21, SIAM, Philadelphia, PA, 2000. [2] P.J.W. Debye, Polar Molecules, Dover, New York, 1929. [3] P.G. Petropoulos, The wave hierarchy for propagation in relaxing dielectrics, Wave Motion 21 (3) (1995) 253–262. [4] T. Kashiwa, I. Fukai, A treatment by the FD–TD method of the dispersive characteristics associated with electronic polarization, Microw. Opt. Technol. Lett. 3 (6) (1990) 203–205. [5] T. Kashiwa, N. Yoshida, I. Fukai, A treatment by the finite-difference time domain method of the dispersive characteristics associated with orientational polarization, IEICE Trans. 73 (8) (1990) 1326–1328. [6] A. Taflove, S.C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method, Wiley, 2005.

188

V.A. Bokil et al. / Journal of Computational and Applied Mathematics 263 (2014) 160–188

[7] K. Yee, Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media, IEEE Trans. Antennas Propag. 14 (3) (1966) 302–307. [8] J.C. Strikwerda, Finite Difference Schemes and Partial Differential Equations, SIAM, 2004. [9] P.G. Petropoulos, Stability and phase error analysis of FDTD in dispersive dielectrics, IEEE Trans. Antennas Propag. 42 (1) (1994) 62–69. [10] H.K. Rouf, F. Costen, S.G. Garcia, 3D Crank–Nicolson finite difference time domain method for dispersive media, Electron. Lett. 45 (19) (2009) 961–962. [11] H.K. Rouf, F. Costen, S.G. Garcia, S. Fujino, On the solution of 3-D frequency dependent Crank–Nicolson FDTD scheme, J. Electromagn. Waves 23 (16) (2009) 2163–2175. [12] S. González García, R. Godoy Rubio, A. Rubio Bretones, R. Gómez Martín, Extension of the ADI–FDTD method to Debye media, IEEE Trans. Antennas Propag. 51 (11) (2003) 3183. [13] T. Namiki, A new FDTD algorithm based on alternating-direction implicit method, IEEE Trans. Microw. Theory Tech. 47 (10) (1999) 2003–2007. [14] J.A. Pereda, A. Grande, O. Gonzalez, A. Vegas, An efficient 2D ADI–FDTD modeling of high-order, frequency-dependent media, IEEE Antennas Wirel. Propag. Lett. 8 (2009) 724–727. [15] O. Ramadan, DSP-based ADI–PML formulations for truncating linear Debye and Lorentz dispersive FDTD domains, in: Computational Science and its Applications—ICCSA 2005, Springer, 2005, pp. 926–935. [16] O. Ramadan, General ADI–FDTD formulations for multi-term dispersive electromagnetic applications, IEEE Microw. Wirel. Compon. Lett. 21 (10) (2011) 513–515. [17] H.K. Rouf, F. Costen, S.G. Garcia, Reduction of numerical errors in frequency dependent ADI–FDTD, Electron. Lett. 46 (7) (2010) 489–490. [18] G.I. Marchuk, Some application of splitting-up methods to the solution of mathematical physics problems, Appl. Math. 13 (2) (1968) 103–132. [19] W. Chen, X. Li, D. Liang, Energy conserved splitting FDTD methods for Maxwell’s equations, Numer. Math. 108 (3) (2008) 445–485. [20] J. Lee, B. Fornberg, A split step approach for the 3-D Maxwell’s equations, J. Comput. Appl. Math. 158 (2) (2003) 485–505. [21] P. Monk, A comparison of three mixed methods for the time-dependent Maxwell’s equations, SIAM Sci. Stat. Comput. 13 (1992) 1097–1122. [22] S. Lanteri, C. Scheid, Convergence of a discontinuous Galerkin scheme for the mixed time-domain Maxwell’s equations in dispersive media, IMA J. Numer. Anal. (2012). [23] J. Li, Unified analysis of leap-frog methods for solving time-domain Maxwell’s equations in dispersive media, J. Sci. Comput. 47 (1) (2011) 1–26. [24] O.A. Keefer, Operator splitting methods for Maxwell’s equations in dispersive media, Master’s Thesis, Oregon State University, 2012. [25] I. Faragó, A. Havasiy, Operator Splittings and their Applications, Nova Science Publishers, 2009. [26] B. Sportisse, An analysis of operator splitting techniques in the stiff case, J. Comput. Phys. 161 (1) (2000) 140–168. [27] G. Strang, On the construction and comparison of difference schemes, SIAM J. Numer. Anal. 5 (3) (1968) 506–517. [28] V.A. Bokil, N.L. Gibson, Analysis of spatial high-order finite difference methods for Maxwell’s equations in dispersive media, IMA J. Numer. Anal. 32 (3) (2012) 926–956. [29] B. Bidégaray-Fesquet, Stability of FD–TD schemes for Maxwell–Debye and Maxwell–Lorentz equations, SIAM J. Numer. Anal. 46 (5) (2008) 2551–2566. [30] V.A. Bokil, N.L. Gibson, Convergence analysis of Yee schemes for Maxwell’s equations in Debye and Lorentz dispersive media, Int. J. Numer. Anal. Model. (2014) in press. [31] R.A. Albanese, R.L. Medina, J.W. Penn, Mathematics, medicine and microwaves, Inverse Probl. 10 (1994) 995–1007. [32] H.T. Banks, V.A. Bokil, N.L. Gibson, Analysis of stability and dispersion in a finite element method for Debye and Lorentz dispersive media, Numer. Methods Partial Differential Equations 25 (4) (2009). [33] N.N. Yanenko, The Method of Fractional Steps: The Solution of Problems of Mathematical Physics in Several Variables, Springer, 1971. [34] V.A. Bokil, A.C.-Y. Leung, A Sequential Operator Splitting Method for Electromagnetic Wave Propagation in Dispersive Media, Technical Report ORSTMATH-11-03, Oregon State University, 2010, http://hdl.handle.net/1957/23059.