PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Comput Phys Commun. Author manuscript; available in PMC 2010 June 16.
Published in final edited form as:
Comput Phys Commun. 2008 December 1; 179(11): 791–800.
doi:  10.1016/j.cpc.2008.07.001
PMCID: PMC2886521
NIHMSID: NIHMS121792

On the immersed interface method for solving time-domain Maxwell’s equations in materials with curved dielectric interfaces

Abstract

This paper deals with accurate numerical simulation of two-dimensional time-domain Maxwell’s equations in materials with curved dielectric interfaces. The proposed fully second-order scheme is a hybridization between the immersed interface method (IIM), introduced to take into account curved geometries in structured schemes, and the Lax–Wendroff scheme, usually used to improve order of approximations in time for partial differential equations. In particular, the IIM proposed for two-dimensional acoustic wave equations with piecewise constant coefficients [C. Zhang, R.J. LeVeque, The immersed interface method for acoustic wave equations with discontinuous coefficients, Wave Motion 25 (1997) 237–263] is extended through a simple least squares procedure to such Maxwell’s equations. Numerical results from the simulation of electromagnetic scattering of a plane incident wave by a dielectric circular cylinder appear to indicate that, compared to the original IIM for the acoustic wave equations, the augmented IIM with the proposed least squares fitting greatly improves the long-time stability of the time-domain solution. Semi-discrete finite difference schemes using the IIM for spatial discretization are also discussed and numerically tested in the paper.

Keywords: Immersed interface method, Lax-Wendroff scheme, Least squares fitting, Maxwell’s equations, Acoustic wave equations

1. Introduction

The second-order finite difference, Cartesian grid-based immersed interface method (IIM), originally developed by LeVeque and Li for elliptic equations with discontinuous coefficients and/or singular source terms [1,2], has been applied to a broad range of equations with discontinuous coefficients and practical problems with irregular boundaries; see [3,4] and references therein. In particular, the method has been extended by Zhang and LeVeque to solve hyperbolic systems of partial differential equations with piecewise constant coefficients, including two-dimensional acoustic wave equations and elasticity equations [5,6]. In this note, we shall further extend the IIM for the acoustic wave equations to two-dimensional time-domain Maxwell’s equations in materials with curved dielectric interfaces.

A number of second-order finite difference methods have been proposed in the past for modeling time-domain Maxwell’s equations in inhomogeneous media with irregular material interfaces. The usual and straightforward approach is to introduce appropriate local modifications into the famous staggered-grid Yee scheme [7] but still keep the staggered grid [810]. It should be mentioned that there are also some studies of fourth-order staggered-grid finite difference time-domain (FDTD) schemes for the Maxwell’s equations with material interfaces [1115]. Subtle interface techniques with one-sided difference approximations and extrapolations are employed in these methods. The significant accomplishment of these fourth-order schemes lies in the fact that the physical jump conditions at the material interfaces are correctly enforced up to fourth-order, so that high-order convergence is uniformly assured over the entire computational domain. The implementation of these high-order FDTD methods for Maxwell’s equations in inhomogeneous media with complex material interfaces, however, has not been addressed with satisfactory results in the literature and therefore remains a great challenge.

Previous second-order finite difference methods making use of Cartesian grids for the approximation of time-domain Maxwell’s equations in media with material interfaces can be found in [16,17]. In addition to one-sided difference approximations and extrapolations, solutions at both sides of an interface are kept track of and updated in a manner that the upwind characteristic information is unaffected by the material interface. When employed to solve the two-dimensional Maxwell’s equations in inhomogeneous media, albeit they are shown to be stable for time integration long enough to resolve the wave propagation in the media, for very long time integration instability has been observed for some applications. In addition, a higher-order Cartesian grid-based FDTD method via the so-called hierarchical implicit derivative matching is presented in [18]. Again no implementation of this higher-order FDTD method for Maxwell’s equations in inhomogeneous media with complex material interfaces, however, has yet been published so how efficient the method will be when solving such equations remains an interesting open question.

In this note, we shall investigate the IIM proposed by Zhang and LeVeque [5] to time-domain Maxwell’s equations in media with material interfaces. The direct adoption of the method to such Maxwell’s equations, however, shows long-time instability of the time-domain solution for some applications. To ameliorate this unsatisfactory stability degradation, we propose an augmented IIM based on a simple least squares procedure. Numerical results from the simulation of electromagnetic scattering of a plane incident wave by a dielectric circular cylinder then appear to indicate that, compared to the original IIM for the acoustic wave equations, the augmented IIM with the proposed least squares fitting greatly improves the stability of the time-domain solution. In addition, semi-discrete finite difference schemes using the least squares-based IIM for spatial discretization are discussed and numerically tested as well.

2. IIM for two-dimensional Maxwell’s equations

In [5], the IIM has been developed for the two-dimensional acoustic wave equations with piecewise constant coefficients, which can be written as a linear first-order hyperbolic system

Ut+AUx+BUy=0,
(1)

where

U=(u(x,y,t)v(x,y,t)p(x,y,t)),A=(001/ρ000K00),B=(000001/ρ0K0).

Here ρ(x, y), the density, and K(x, y), the bulk modulus of elasticity, are piecewise constant with jump discontinuities at some interface(s). The pressure perturbation p(x, y, t) and the velocity u = (u(x, y, t), v(x, y, t))T are therefore generally non-smooth or even discontinuous across such an interface. For example, typically the tangential velocity will have a jump discontinuity, while the normal velocity and the pressure will have jumps in their derivatives, across the interface, leading to the following jump conditions on the interface:

[n·u]=0,[t·ρu]=0,[p]=0,
(2)

where n and t are the unit vectors in the directions normal and tangential to the interface, respectively, and the notation [ ] represents the difference of the values of a quantity on both sides of the interface, usually named as the jump of the quantity.

Under the general assumption that uyvx [equivalent] 0 (the deformation is irrotational) which remains true provided it is true of the initial data, by following the same procedure as proposed by LeVeque and Li [1] for elliptic equations with discontinuous coefficients and/or singular source terms, a second-order IIM on a uniform Cartesian grid has been developed in [5] for the acoustic wave equations (1) with the jump conditions (2).

In this note, we consider the two-dimensional z-transverse magnetic (TMz) set of Maxwell’s equations

Hxt+1μEzy=0,Hyt1μEzx=0,Ezt+1ε(HxyHyx)=0,
(3)

where Ez and H = (Hx, Hy)T represent the scalar electric field and the vector magnetic field, respectively, and ε and μ represent the material electric permittivity and the material magnetic permeability, respectively. In addition, the solution domain Ω is divided by a material/dielectric interface into two disjoint regions, Ω and Ω+, representing two distinct dielectric materials. The boundary conditions at the material interface are

[t·H]=0,[n·μH]=0,[Ez]=0.
(4)

Now note that if we let

u=Hy,v=Hx,p=Ez,ρ=μ,K=1/ε,

then the two-dimensional Maxwell’s equations (3) can be written in the form of the two-dimensional acoustic wave equations (1), and the boundary conditions (4) become the jump conditions (2). Also, from [nabla] · H = 0, we have uyvx [equivalent] 0. Therefore, in principle the IIM developed for the acoustic wave equations with piecewise constant coefficients can be directly adopted for us to solve the time-domain Maxwell’s equations in media with material interfaces, which is summarized below.

2.1. Lax–Wendroff scheme-based IIM

Let the solution domain Ω be discretized into a uniform Cartesian grid as illustrated in Fig. 1. Let further, for simplicity, h be the spatial grid size in both the x- and y-directions. Away from the interface standard finite difference methods can be used, for example the second-order Lax–Wendroff scheme which, for a general hyperbolic system (1), takes the form

Fig. 1
A diagram of a two-dimensional interface embedded in a uniform Cartesian grid and a coordinate system (ξ, η) defined locally at a point P on the interface.
Uijn+1=Uijnk2h(A(Ui+1,jnUi1,jn)+B(Ui,j+1nUi,j1n))+12(kh)2(A2(Ui+1,jn2Uijn+Ui1,jn)+B2(Ui,j+1n2Uijn+Ui,j1n)+14(AB+BA)(Ui+1,j+1nUi+1,j1nUi1,j+1n+Ui1,j1n)).
(5)

Here and in the sequel, k denotes the step size in time. The scheme (5) is second-order accurate at regular grid points. Here a grid point (xi, yj) is said to be regular if all nine grid points in the standard 9-point stencil (5) centered at (xi, yj) lie on the same side of the interface. On the contrary, a grid point (xi, yj) is said to be irregular if the nine grid points in the standard 9-point stencil (5) centered at (xi, yj) are in the different sides of the interface. For example, under these definitions, Point 2 of Fig. 1 is a regular grid point and Points 1, 3, 4, 5, 6 and A all are irregular grid points, respectively.

Moreover, for the Maxwell’s equations (3), the cross derivative terms in (5) can be eliminated since

Utt=A2Uxx+(AB+BA)Uxy+B2Uyy=c2(Uxx+Uyy),

Where c = (k/ρ)1/2 = 1/(εμ)1/2. This allows us to use instead the 5-point second-order accurate Lax–Wendroff scheme

Uijn+1=Uijnk2h(A(Ui+1,jnUi1,jn)+B(Ui,j+1nUi,j1n))+12(kh)2c2(Ui+1,jn+Ui1,jn+Ui,j+1n+Ui,j1n4Uijn)
(6)

at regular grid points. Note that now a grid point (xi, yj) is said to be regular if all five grid points in the standard 5-point stencil (6) centered at (xi, yj) lie on the same side of the interface, and irregular if otherwise. Therefore, this not only simplifies the formulation hereafter but also reduces the number of irregular grid points where special treatment is needed. For instance, Point A shown in Fig. 1 is irregular under (5) but becomes regular under (6).

On the other hand, at irregular grid points the difference equations are modified to take into account the discontinuity in the coefficients and the resulting non-smoothness in the solution. More precisely, at each of these irregular grid points, say (xi, yj), the following 6-point difference scheme

Uijn+1=Uijn+khl=1MΓij,lUij,ln
(7)

with M = 6 is used, where the coefficients Γij,l are 3 × 3 matrices, and Uij,ln for l = 1, 2,…, M are solutions at some set of grid points neighboring the point (xi, yj). In order for (7) to reduce to the standard Lax–Wendroff scheme (6) when no interface is present, the six grid points in the stencil (7) are chosen to always include those in the standard 5-point stencil (namely, Points 1, 2, 3, 4 and 5) together with one additional corner grid point neighboring (xi, yj), as indicated in Fig. 1.

The Γ matrices are so determined that the second-order accuracy is maintained, which can be achieved by ensuring the local truncation error to satisfy

T=1hl=1MΓij,lUij,l(Ut(xi,yj)+12kUtt(xi,yj))+O(k2).
(8)

In more detail, to calculate the coefficient matrices Γij,l for an irregular grid point (xi, yj), we first let P = (x*, y*) be a point on the interface near (xi, yj) as shown in Fig. 1. Basically we can take P as any point on the interface that is in the neighborhood of (xi, yj), but usually we take this point either as the orthogonal projection of (xi, yj) on the interface or the intersection of the interface and one of the two lines x = xi and y = yj. Then a local coordinate system (ξ, η) centering at P is defined by the coordinate transformation

(ξη)=(cosθsinθsinθcosθ)(xxyy),
(9)

where ξ and η are in the directions normal and tangential to the interface at P, respectively, and θ is the rotation angle. We assume the interface is smooth in a neighborhood of the point P and can be described locally by the relation ξ = χ(η), with χ(0) = χ′(0) = 0, and χ″(0), the curvature, to be finite.

The system (1) can then be written in the local coordinates and the rotated equations have the same form due to isotropy

U¯t+AU¯ξ+BU¯η=0,
(10)

Where

U¯(ξ,η,t)=(u¯(ξ,η,t)v¯(ξ,η,t)p(ξ,η,t))=Q01U(x,y,t),Q0=(cosθsinθ0sinθcosθ0001).

Here ū = n · u and v = t · u are the components of the velocity in the ξ- and η-directions, respectively. Consequently, the jump conditions in terms of the transformed quantities are given by

[u¯]=0,[ρv¯]=0,[p]=0.
(11)

Like in the derivation of any IIM, a knowledge of jump conditions for the solution and its derivatives at the interface is required. Using the jump conditions (11) together with the differential equations (10), the interface conditions between + and − values of the rotated solution and its derivatives at the interface are found to be as follows [5]:

U¯+=Q1U¯,U¯η+=Q1U¯η,U¯ξ+=Q2U¯ξ+Q3U¯η,U¯ηη+=Q1U¯ηηχ(0)((Q2Q1)U¯ξ+Q3U¯η),U¯ξη+=Q2U¯ξη+Q3U¯ηη+χ(0)Q4((Q2Q1)U¯ξ+Q3U¯η),U¯ξξ+=(cc+)2Q1U¯ξξ+((cc+)21)Q1U¯ηη+χ(0)((Q2Q1)U¯ξ+Q3U¯η),
(12)

where Qi, i = 1, 2, 3, 4, are constant matrices given by

Q1=diag(1,ρ/ρ+,1),Q2=diag(K/K+,1,ρ+/ρ),Q3=ρρ+((cc+)21)(010000000),Q4=(010100000).
(13)

Then by following the standard procedure of an IIM one can determine the unknown coefficient matrices Γij,l. Passing the details which can be found in [5], we have that the coefficient matrices Γij,l in the 6-point stencil (7) are calculated by

Γij,l=Q0Γ^ij,lQ01,

where [Gamma]ij,l are obtained by solving the following system of six matrix equations

l=1MΓ^ij,lQ^1l=0,l=1MΓ^ij,lQ^4l=khc2I2ξ1hA,l=1MΓ^ij,lQ^2l=A,l=1MΓ^ij,lQ^5l=2h(η1A+ξ1B),l=1MΓ^ij,lQ^3l=B,l=1MΓ^ij,lQ^6l=khc2I2η1hB.
(14)

Here c = c, or c+ depending on which side the center point of the stencil sits, (ξ1, η1) are the local coordinates of the center point (xi, yj), the Q’s are 3 × 3 matrices, and I is the 3 × 3 identity (xi, yj), the matrix. If the lth grid point in the stencil (7) is on the same side as the center point (xi, yj), then

Q^1l=I,Q^2l=ξlhI,Q^3l=ηlhI,Q^4l=(ξlh)2I,Q^5l=2ξlηlh2I,Q^6l=(ηlh)2I.
(15)

Otherwise, we have

Q^1l=Q1,Q^2l=ξlhQ2+12hχ(0)(ξl2ηl2+2ξlηlQ4)(Q2Q1),Q^3l=ηlhQ1+ξlhQ3+12hχ(0)(ξl2ηl2+2ξlηlQ4)Q3,Q^4l=(ξlchc+)2Q1,Q^5l=2ξlηlh2Q2,Q^6l=(ηlh)2Q1+2ξlηlh2Q3+(ξlh)2((cc+)21)Q1.
(16)

Here, (ξ1, η1) are the local coordinates of the lth grid point in the stencil (7).

Finally, for completeness, let us also mention the two-dimensional z-transverse electric (TEz) set of Maxwell’s equations,

Ext1εHzy=0,Eyt+1εHzx=0,Hzt1μ(ExyEyx)=0,
(17)

where Hz and E = (Ex, Ey)T represent the magnetic scalar field and the vector electric field, respectively. The corresponding boundary conditions at the material interface become

[t·E]=0,[n·εE]=0,[Hz]=0.
(18)

Likely, if we let

u=Ey,v=Ex,p=Hz,ρ=ε,K=1/μ,

then the Maxwell’s equations (17) can be written again in the form of the acoustic wave equations (1), and the boundary conditions (18) become the jump conditions (2). Also, uyvx [equivalent] 0 since [nabla] · E [equivalent] 0. Therefore, in principle the foregoing IIM developed for the acoustic wave equations can also be utilized to solve the TEz set of Maxwell’s equations (17).

2.2. Lax–Wendroff scheme-based IIM with least squares fitting

When M = 6, Eq. (14) forms a determined linear system of equations for 54 unknowns. First of all, the existence of the solution to such a small linear system cannot always be guaranteed, although in general it should not be a problem. But most importantly, even though the linear system (14) is solvable for each irregular grid point, the resulting finite difference scheme (7) may be unstable, in particular when the curvature of the interface and/or the difference between the material constants are relatively large. As a matter of fact, like to be demonstrated in Section 3, the Lax–Wendroff scheme-based IIM using the 6-point stencil (7) indeed becomes unstable when used to simulate electromagnetic scattering of a plane incident wave by a dielectric circular cylinder with ε = 2.25 and μ = 2 embedded in the vacuum.

In hoping to resolve these potential problems at least to some extent, instead of using only six grid points, we can choose to use M > 6 grid points in the stencil (7). It should be pointed out that this basic idea is not new as it has been applied in the IIM for elliptic equations [19]. In this case, (14) is an underdetermined linear system, and therefore will have an infinite number of solutions. We could choose the one with some minimum 2-norm. To this end, we form the following constrained quadratic optimization problem to determine the coefficient matrices Γij,l of the finite difference scheme (7)

min12l=1M(Γij,lDl)2
(19)

subject to (14). Additional constraints may be imposed for further improvements. We also want to calculate the matrices Γij,l in such a way that, if no material interface is present, the finite difference scheme (7) reduces to the standard 5-point Lax–Wendroff scheme (6) which can be rewritten as

Uijn+1=Uijn+kh[(k2hc2I+12A)Ui1,jn+(k2hc2I12A)Ui+1,jn+(k2hc2I+12B)Ui,j1n+(k2hc2I12B)Ui,j+1n2khc2Ui,jn].

This can be done by always including the five grid points of the standard 5-point stencil in the general M-point stencil, numbering them as indicated in Fig. 1, and then selecting the matrices Dl in the optimization problem (19) as

D1=2khc2I,D2=k2hc2I+12A,D3=k2hc2I12A,D4=k2hc2I+12B,D5=k2hc2I12B,Dl=0,l=6,,M.
(20)

One natural concern with the IIM with the proposed least squares fitting is the existence of the solution to the constrained quadratic optimization problem (19). For problems with general curved interfaces, there seems to be little hope of proving existence rigorously. The situation here is more complicated than in [19], partially because here we are looking for unknown 3 × 3 matrices instead of just unknown scalars in [19]. Nevertheless, numerous computational experiments have shown existence of the solution to this minimization problem, and in all numerical experiments we have conducted, we have experienced no case whatsoever in which the corresponding minimization problem fails to return a solution when it is solved by QL0001, a package developed by K. Schittkowski to solve constrained quadratic optimization problems.

Similarly, a full-blown theoretical analysis of the stability of the IIM seems to be unrealistic, particularly for problems with general curved interfaces. However, numerical results from the simulation of the same electromagnetic scattering problem appear to indicate that, compared to the original IIM using only six points in the stencil, the augmented IIM with the proposed least squares fitting greatly improves the stability of the time-domain solution.

2.3. Semi-discrete scheme-based IIM

The basic idea of the IIM shall also allow us to develop semi-discrete finite difference schemes, either the Cartesian grid-based or the staggered grid-based, for time-domain Maxwell’s equations. This approach may be particularly needed in developing higher-order FDTD methods for Maxwell’s equations in inhomogeneous media with complex material interfaces, which so far have not yet been addressed with substantial success in the literature. However, it is not our intent in this note to develop such higher-order methods. Rather, we simply want to give an example of how to apply the basic idea of the IIM with least squares fitting in discretizing the solution in space and then to show that using more grid points in the discretization of space again leads to better long-time stability of the overall finite difference scheme.

At a regular grid point in a Cartesian grid, the standard central finite difference is utilized to discretize the solution in space

(AUx+BUy)ij=A(Ui+1,jUi1,j2h)+B(Ui,j+1Ui,j12h).
(21)

On the other hand, at an irregular grid point (xi, yj) the difference equation

(AUx+BUy)ij=1hl=1MΓij,lUij,l
(22)

is used, where the coefficients Γij,l are 3 × 3 matrices, and Uij,l for l = 1, 2, …, M are solutions at some set of grid points neighboring the point (xi, yj). The Γ matrices are so determined that the second-order accuracy in the space increment is maintained, which can be accomplished by demanding the local truncation error to satisfy

T=1hl=1MΓij,lUij,l(AUx(xi,yj)+BUy(xi,yj))+O(h).
(23)

Then following the same procedure as outlined in [5] leads to

Γij,l=Q0Γ^ij,lQ01,

where in general the [Gamma]ij,l matrices are obtained by solving the constrained quadratic optimization problem (19) subject to

l=1MΓ^ij,lQ^1l=0,l=1MΓ^ij,lQ^4l=2ξ1hA,l=1MΓ^ij,lQ^2l=A,l=1MΓ^ij,lQ^5l=2h(η1A+ξ1B),l=1MΓ^ij,lQ^3l=B,l=1MΓ^ij,lQ^6l=2η1hB.
(24)

Here the Q matrices are identical to those in (15) or (16). And for (22) to reduce to the standard central finite difference (21) when there is no material interface, again the five grid points in the standard 5-point stencil (21) shall be included in (22), and the corresponding Dl matrices in the optimization problem (19) shall be chosen as

D2=12A,D3=12A,D4=12B,D5=12B,Dl=0,otherwise.

3. Numerical results

To evaluate the proposed IIMs for time-domain Maxwell’s equations in media with curved material interfaces, we shall consider a typical electromagnetic scattering problem as illustrated in Fig. 2, namely, scattering of a plane incident wave by a dielectric circular cylinder embedded in the vacuum. The Maxwell’s equations (3), together with the boundary conditions (4), can be used to solve the problem.

Fig. 2
Scattering of a plane incident wave by a dielectric circular cylinder.

We assume that the cylinder is illuminated by a time-harmonic incident plane wave of the form

Eincz=ei(k1xωt),Hincy=ei(k1xωt),

where k1 = ω(μ1ε1)1/2 is the propagation constant for the homogeneous, isotropic free space medium. An analytical solution of this problem due to Mie can be found in [20] and is also available in [17].

We would like to verify the second-order convergence of the proposed IIMs for solving the above practical scattering problem. To this end, the analytical solution of the problem is imposed as the initial condition as well as the Dirichlet boundary condition at the artificial boundary of the computational domain. Also, in our simulations we let the cylinder radius a be 0.5 and the angular frequency ω incident wave be 2π. The computational domain is the square Ω = [−1, 1] × [−1, 1], and the time step size is imposed as k=h/(22λmax) and k = h/(10λmax) for the Lax–Wendroff scheme-based IIM and the semi-discrete scheme-based one, respectively, where λmax = max{c1, c2} and ci = 1/(εiμi)1/2, for i = 1, 2 in our tests. Please note that λmax = 1 in our simulations since the material exterior to the cylinder is assumed to be the vacuum (ε1 = 1 and μ1 = 1).

3.1. Lax–Wendroff scheme-based IIM

Let us begin by considering the situation in which ε2 = 2.25 and μ2 = 1, i.e., the cylinder is non-magnetic. Fig. 3 shows the contour and the x = 0.2 cross-section of the computed field component Hx at the time t = 100 obtained by the Lax–Wendroff scheme-based IIM with M = 25 and a 400 × 400 mesh. Please note that Hx is continuous but its derivative is discontinuous across the material interface. Actually, in this case, it is clear from the interface conditions (4) that all three field components are continuous across the material interface. The derivative of the Ez component is also continuous, but the derivatives of the Hx and Hy components are discontinuous, across the interface.

Fig. 3
Scattering of a plane incident wave by a dielectric circular cylinder with the material parameters ε2 = 2.25 and μ2 = 1. (a) The contour of the computed field component Hx at the time t = 100 when using a 400 × 400 mesh, and (b) ...

Table 1 shows the error analysis results for the Lax–Wendroff scheme-based IIM with M = 6, M = 9, and M = 25, respectively. In our simulations, at an irregular grid point (xi, yj), the sixth grid point in the 6-point stencil is chosen as the closest corner point to (x*, y*) which is always taken as the orthogonal projection of (xi, yj) on the interface. The nine or twenty-five grid points in the 9- or 25-point stencil are those defined in a typical compact 9- or 25-point stencil on a regular square grid. ||E|| denotes the maximal relative error in the numerical solution measured in L norm over all grid points at the time t = 100. Studying this error as a function of the grid/time step sizes, the global second-order accuracy is clearly observed for all three stencils, and moreover it can be noted that the three stencils have comparable accuracy.

Table 1
Convergence analysis for the Lax–Wendroff scheme-based IIM (ε2 = 2.25 and μ2 = 1)

We then consider the situation in which ε2 = 2.25 and μ2 = 2. Fig. 4 shows the contour and the y = 0.2 cross-section of the computed field component Hx at the time t = 100 obtained by the Lax–Wendroff scheme-based IIM with M = 25 and a 400 × 400 mesh. Please note that in this case, both Hx and its derivative are discontinuous across the interface. Actually, it is clear again from the interface conditions (4) that in this case only the field component Ez is continuous, and the field components Hx and Hy and the derivatives of all three field components are discontinuous across the material interface.

Fig. 4
Scattering of a plane incident wave by a dielectric cylinder with the material parametersε2 = 2.25 and μ2 = 2. (a) The contour of the computed field component Hx at the time t = 100 when using a 400 × 400 mesh, and (b) the y = ...

Table 2 shows the error analysis results again at the time t = 100 for the Lax–Wendroff scheme-based IIM with M = 6, M = 9, and M = 25, respectively, where ∞ indicates a value greater than10300 and “div” stands for “divergent”. First as clearly shown in Table 2, unlike for the μ2 = 1 case, the Lax–Wendroff scheme-based IIM with M = 6 becomes numerically unstable for the case of μ2 = 2. The errors grow very rapidly without limit with time. Several other ways to choose the sixth grid point in the 6-point stencil and to impose the boundary condition at the artificial boundary are also tested, but none could result in a stable 6-point stencil.

Table 2
Convergence analysis for the Lax–Wendroff scheme-based IIM (ε2 = 2.25 and μ2 = 2)

Nevertheless, the Lax–Wendroff scheme-based IIM with M = 9 or M = 25 still maintains the long-time stability. As can be seen from Table 2, first the global second-order accuracy is clearly maintained for the IIM with M = 9 or M = 25. Second, compared to for the situation of μ2 = 2, the numerical solution for the situation of μ2 = 1 retains better accuracy, which is because the wave length inside the cylinder with μ2 = 2 is a factor of 2 smaller than that with μ1 = 1, and thus a smaller grid size shall be needed to resolve the wave propagation to the same degree of accuracy.

In order to gain a clearer picture of how the errors of the numerical solutions evolve, the time histories of the errors of the Lax–Wendroff scheme-based IIM using different numbers of stencil points are recorded and displayed in Fig. 5. Here, N denotes the mesh size in both the x- and y-directions. Note that for the Lax–Wendroff scheme-based IIM using the least squares procedure with M = 9 and M = 25, for both μ2 = 1 and μ2 = 2 cases there is no noticeable growth of the error after running the simulation for 100 time units (see Fig. 5(c), (d)), indicating the scheme will remain stable and second-order accurate for long time computations. For the Lax–Wendroff scheme-based IIM with M = 6 (namely, the original IIM proposed for acoustic wave equations [5]), for the case of μ2 = 1, the method shows the same long-time stability behavior without noticeable growth of the errors (see Fig. 5(a)). For the case of μ2 = 2, however, the method becomes numerically unstable and the errors grow very rapidly without limit with time. As indicated by Fig. 5(b), even for a very short time integration, the numerical solution shall quickly diverge from the actual one. In summary, compared to the original IIM for the acoustic wave equations, the augmented IIM with the proposed least squares fitting greatly improves the long-time stability of the time-domain solution.

Fig. 5
The maximal relative error in L norm in the numerical solution as a function of time for the Lax–Wendroff scheme-based IIM. (a) M = 6 (t ε [0, 100]), (b) M = 6 (t ε [0, 1]), (c) M = 9 (t ε [0, 100]), (d) M = 25 ...

3.2. Semi-discrete scheme-based IIM

Fig. 6 shows the time histories of the errors of the numerical solution for the semi-discrete scheme-based IIM using different numbers of stencil points in (22). In all numerical experiments, the two-stage second-order mid-point Runge–Kutta method is used for time integration, which can be described by

Fig. 6
The maximal relative error in L norm in the numerical solution as a function of time for the mid-point Runge–Kutta scheme-based IIM. (a) M = 6 (t ε [0, 1]), (b) M = 9 (t ε [0, 10]), (c) M = 9 (t ε [0, 1]), (d) ...

un+1=un+Δtf(tn+Δt2,un+Δt2f(tn,un))

when solving the following initial value problem

dudt=f(t,u),u(t0)=u0.

It should be mentioned that two other second-order time integration schemes, the improved Euler method and the truncation error-optimal Runge–Kutta method [21], are also employed in our numerical tests, and the same conclusions to be discussed below follow.

First of all, for the semi-discrete scheme-based IIM, it is again observed that using more grid points in the spatial stencil (22) improves the stability of the overall semi-discrete finite difference scheme. First as illustrated in Fig. 6(a), the semi-discrete scheme using M = 6 clearly shows non-convergent behavior for both the μ2 = 1 and μ2 = 2 cases. Then as shown in Fig. 6(b), for the case of μ2 = 1, the numerical solution obtained by the semi-discrete scheme using M = 9 converges at t = 1, and after that the error grows gradually but not as fast as that produced by the semi-discrete scheme using M = 6. Next as evidenced by Fig. 6(c), for the case of μ2 = 2, the semi-discrete scheme using M = 9 also exhibits non-convergent behavior. And finally as demonstrated by Fig. 6(d), for the case of μ2 = 1, the semi-discrete scheme-based IIM using M = 25 clearly holds the long-time stability since there is no noticeable growth of the error over a long time integration to t = 100. On the other hand, for the case of μ2 = 2, the semi-discrete scheme-based IIM using M = 25 still maintains relatively long-time stability because there is still no noticeable growth of the error until t = 85 or so but after that the error ||E|| starts to grow and keeps growing for t > 100. It should be pointed out that the growing errors are not spatially localized only close to irregular grid points due to the wave propagation nature of the problem, and the phenomenon remains even if smaller time step sizes are used.

Next, the error analysis results for the mid-point Runge–Kutta scheme-based IIM using M = 25 are listed in Table 3. For the case of μ2 = 1, the second-order convergence rate with comparable accuracy to the Lax–Wendroff scheme-based IIM is clearly observed. For the case of μ2 = 2, the last order corresponding to the 400 × 400 mesh is only 1.1. The reason is that, as shown in Fig. 6(d) and discussed above, the semi-discrete scheme-based IIM with M = 25 has only relatively long-time stability. In particular, the error ||E|| starts to grow after t = 85 and those errors listed in Table 3 are measured at t = 100.

Table 3
Convergence analysis for the mid-point Runge–Kutta scheme-based IIM with M = 25 (ε2 = 2.25 and μ2 = 1 or 2)

Finally, by comparing the numerical results for the Lax–Wendroff scheme-based IIM to those for the semi-discrete scheme-based one, it is readily seen that the former is more stable if the same number of grid points are included in the corresponding finite difference stencils. Whether this is true or not when they are extended to higher-order methods, if possible at all, will remain an open question.

4. Conclusions

In this note, the Lax–Wendroff scheme-based IIM proposed for two-dimensional acoustic wave equations with piecewise constant coefficients is applied to time-domain Maxwell’s equations in materials with curved dielectric interfaces. To improve the long-time stability, we have proposed an augmented IIM based on a least squares procedure. Numerical results from the simulation of electromagnetic scattering of a plane incident wave by a dielectric circular cylinder have demonstrated that, compared to the original IIM for the acoustic wave equations, the augmented Lax–Wendroff scheme-based IIM with the proposed least squares fitting greatly improves the stability of the time-domain solution. The author thus hopes that this least square-based IIM can be extended for us to further derive both higher-order accurate and long-time stable FDTD methods for time-domain Maxwell’s equations in inhomogeneous media with complex material interfaces. The extension will be a subject of the author’s future work, and any significant progresses in this endeavor will be reported in future publications.

Acknowledgments

The author thanks the support of the National Institutes of Health (grant number: 1R01GM083600-01) for the work reported in this note. This work was also supported, in part, by funds provided by the University of North Carolina at Charlotte. The author also thanks the anonymous referee for his many productive suggestions on the earlier draft of this paper.

References

1. LeVeque RJ, Li Z. The immersed interface method for elliptic equations with discontinuous coefficients and singular sources. SIAM J Numer Anal. 1994;31:1019–1044.
2. Li Z. PhD thesis. University of Washington; 1994. The immersed interface method—a numerical approach for partial differential equations with interfaces.
3. Li Z, Ito K. Frontiers in Applied Mathematics. SIAM; Philadelphia: 2006. The Immersed Interface Method: Numerical Solutions of PDEs Involving Interfaces and Irregular Domains.
4. Li Z. An overview of the immersed interface method and its applications. Taiwanese J Math. 2003;7:1–49.
5. Zhang C, LeVeque RJ. The immersed interface method for acoustic wave equations with discontinuous coefficients. Wave Motion. 1997;25:237–263.
6. Zhang C. PhD thesis. University of Washington; 1996. Immersed interface methods for hyperbolic systems of partial differential equations with discontinuous coefficients.
7. Yee KS. Numerical solution of initial boundary value problems involving Maxwell equations in isotropic media. IEEE Trans Antennas Propagat. 1966;14:302–307.
8. Xiao T, Liu QH. A staggered upwind embedded boundary (SUEB) method to eliminate the FDTD staircasing error. IEEE Trans Antennas Propagat. 2004;52:730–741.
9. Ditkowski A, Dridi K, Hesthaven JS. Convergent Cartesian grid methods for Maxwell’s equations in complex geometries. J Comput Phys. 2001;170:39–80.
10. Ditkowski A, Hesthaven JS, Teng CH. Modeling dielectric interfaces in the FDTD-method: A comparative study. 2000 Progress in Electromagnetics Research (PIERS 2000), Proceedings; Cambridge, MA. 2000.
11. Hesthaven JS. High-order accurate methods in time-domain computational electromagnetics. A review. Adv Imag Electron Phys. 2003;127:59–123.
12. Xie Z, Chan CH, Zhang B. An explicit fourth-order staggered finite-difference time-domain method for Maxwell’s equations. J Comput Appl Math. 2002;147:75–98.
13. Xie Z, Chan CH, Zhang B. An explicit fourth-order orthogonal curvilinear staggered-grid FDTD method for Maxwell’s equations. J Comput Phys. 2002;175:739–763.
14. Yefet A, Petropoulos PG. A staggered fourth-order accurate explicit finite difference scheme for the time-domain Maxwell’s equations. J Comput Phys. 2001;168:286–315.
15. Yefet A, Turkel E. Fourth order compact implicit method for the Maxwell equations with discontinuous coefficients. Appl Numer Math. 2000;33:125–134.
16. Xue C, Deng S. An upwinding boundary condition capturing method for Maxwell’s equations in media with material interfaces. J Comput Phys. 2007;225:342–362.
17. Cai W, Deng S. An upwinding embedded boundary method for Maxwell’s equations in media with material interfaces: 2D case. J Comput Phys. 2003;190:159–183.
18. Zhao S, Wei GW. High-order FDTD methods via derivative matching for Maxwell’s equations with material interfaces. J Comput Phys. 2004;200:60–103.
19. Li Z, Ito K. Maximum principle preserving schemes for interface problems with discontinuous coefficients. SIAM J Sci Comput. 2001;23:339–361.
20. Umashankar K, Taflove A. Computational Electrodynamics. Artech House; Boston: 1993.
21. Atkinson KE. An Introduction to Numerical Analysis. John Wiley & Sons; New York: 1989.