Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2886521

Formats

Article sections

- Abstract
- 1. Introduction
- 2. IIM for two-dimensional Maxwell’s equations
- 3. Numerical results
- 4. Conclusions
- References

Authors

Related links

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.001PMCID: PMC2886521

NIHMSID: NIHMS121792

Shaozhong Deng^{*}

Shaozhong Deng, Department of Mathematics and Statistics, University of North Carolina at Charlotte, Charlotte, NC 28223-0001, USA;

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.

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 [8–10]. 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 [11–15]. 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.

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

$${U}_{t}+A{U}_{x}+B{U}_{y}=0,$$

(1)

where

$$\begin{array}{l}U=\left(\begin{array}{l}u(x,y,t)\hfill \\ v(x,y,t)\hfill \\ p(x,y,t)\hfill \end{array}\right),\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}A=\left(\begin{array}{ccc}0& 0& 1/\rho \\ 0& 0& 0\\ K& 0& 0\end{array}\right),\\ B=\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& 1/\rho \\ 0& K& 0\end{array}\right).\end{array}$$

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:

$$[\mathbf{n}\xb7\mathbf{u}]=0,\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}[\mathbf{t}\xb7\rho \mathbf{u}]=0,\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}[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 *u _{y}* −

In this note, we consider the two-dimensional *z*-transverse magnetic (TM* _{z}*) set of Maxwell’s equations

$$\begin{array}{l}\frac{\partial {H}^{x}}{\partial t}+\frac{1}{\mu}\frac{\partial {E}^{z}}{\partial y}=0,\\ \frac{\partial {H}^{y}}{\partial t}-\frac{1}{\mu}\frac{\partial {E}^{z}}{\partial x}=0,\\ \frac{\partial {E}^{z}}{\partial t}+\frac{1}{\epsilon}\left(\frac{\partial {H}^{x}}{\partial y}-\frac{\partial {H}^{y}}{\partial x}\right)=0,\end{array}$$

(3)

where *E ^{z}* and

$$[\mathbf{t}\xb7\mathbf{H}]=0,\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}[\mathbf{n}\xb7\mu \mathbf{H}]=0,\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}[{E}^{z}]=0.$$

(4)

Now note that if we let

$$u=-{H}^{y},\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}v={H}^{x},\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}p={E}^{z},\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\rho =\mu ,\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}K=1/\epsilon ,$$

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 · **H** = 0, we have *u _{y}* −

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

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.

$$\begin{array}{l}{U}_{ij}^{n+1}={U}_{ij}^{n}-\frac{k}{2h}(A({U}_{i+1,j}^{n}-{U}_{i-1,j}^{n})+B({U}_{i,j+1}^{n}-{U}_{i,j-1}^{n}))\\ +\frac{1}{2}{\left(\frac{k}{h}\right)}^{2}({A}^{2}({U}_{i+1,j}^{n}-2{U}_{ij}^{n}+{U}_{i-1,j}^{n})\\ +{B}^{2}({U}_{i,j+1}^{n}-2{U}_{ij}^{n}+{U}_{i,j-1}^{n})\\ +\frac{1}{4}(AB+BA)({U}_{i+1,j+1}^{n}-{U}_{i+1,j-1}^{n}-{U}_{i-1,j+1}^{n}\\ +{U}_{i-1,j-1}^{n})).\end{array}$$

(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 (*x _{i}, y_{j}*) is said to be regular if all nine grid points in the standard 9-point stencil (5) centered at (

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

$${U}_{tt}={A}^{2}{U}_{xx}+(AB+BA){U}_{xy}+{B}^{2}{U}_{yy}={c}^{2}({U}_{xx}+{U}_{yy}),$$

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

$$\begin{array}{l}{U}_{ij}^{n+1}={U}_{ij}^{n}-\frac{k}{2h}(A({U}_{i+1,j}^{n}-{U}_{i-1,j}^{n})+B({U}_{i,j+1}^{n}-{U}_{i,j-1}^{n}))\\ +\frac{1}{2}{\left(\frac{k}{h}\right)}^{2}{c}^{2}({U}_{i+1,j}^{n}+{U}_{i-1,j}^{n}+{U}_{i,j+1}^{n}+{U}_{i,j-1}^{n}-4{U}_{ij}^{n})\end{array}$$

(6)

at regular grid points. Note that now a grid point (*x _{i}, y_{j}*) is said to be regular if all five grid points in the standard 5-point stencil (6) centered at (

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 (*x _{i}, y_{j}*), the following 6-point difference scheme

$${U}_{ij}^{n+1}={U}_{ij}^{n}+\frac{k}{h}\sum _{l=1}^{M}{\mathrm{\Gamma}}_{ij,l}{U}_{ij,l}^{n}$$

(7)

with *M* = 6 is used, where the coefficients *Γ _{ij,l}* are 3 × 3 matrices, and
${U}_{ij,l}^{n}$ for

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=\frac{1}{h}\sum _{l=1}^{M}{\mathrm{\Gamma}}_{ij,l}{U}_{ij,l}-\left({U}_{t}({x}_{i},{y}_{j})+\frac{1}{2}k{U}_{tt}({x}_{i},{y}_{j})\right)+O({k}^{2}).$$

(8)

In more detail, to calculate the coefficient matrices *Γ _{ij,l}* for an irregular grid point (

$$\left(\begin{array}{l}\xi \hfill \\ \eta \hfill \end{array}\right)=\left(\begin{array}{cc}\hfill cos\theta \hfill & \hfill sin\theta \hfill \\ \hfill -sin\theta \hfill & \hfill cos\theta \hfill \end{array}\right)\left(\begin{array}{l}x-{x}^{\ast}\hfill \\ y-{y}^{\ast}\hfill \end{array}\right),$$

(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

$${\overline{U}}_{t}+A{\overline{U}}_{\xi}+B{\overline{U}}_{\eta}=0,$$

(10)

Where

$$\begin{array}{l}\overline{U}(\xi ,\eta ,t)=\left(\begin{array}{l}\overline{u}(\xi ,\eta ,t)\hfill \\ \overline{v}(\xi ,\eta ,t)\hfill \\ p(\xi ,\eta ,t)\hfill \end{array}\right)={Q}_{0}^{-1}U(x,y,t),\\ {Q}_{0}=\left(\begin{array}{ccc}cos\theta & -sin\theta & 0\\ sin\theta & cos\theta & 0\\ 0& 0& 1\end{array}\right).\end{array}$$

Here *ū* = **n** · **u** and = **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

$$[\overline{u}]=0,\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}[\rho \overline{v}]=0,\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}[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]:

$$\begin{array}{l}{\overline{U}}^{+}={Q}_{1}{\overline{U}}^{-},\\ {\overline{U}}_{\eta}^{+}={Q}_{1}{\overline{U}}_{\eta}^{-},\\ {\overline{U}}_{\xi}^{+}={Q}_{2}{\overline{U}}_{\xi}^{-}+{Q}_{3}{\overline{U}}_{\eta}^{-},\\ {\overline{U}}_{\eta \eta}^{+}={Q}_{1}{\overline{U}}_{\eta \eta}^{-}-{\chi}^{\u2033}(0)(({Q}_{2}-{Q}_{1}){\overline{U}}_{\xi}^{-}+{Q}_{3}{\overline{U}}_{\eta}^{-}),\\ {\overline{U}}_{\xi \eta}^{+}={Q}_{2}{\overline{U}}_{\xi \eta}^{-}+{Q}_{3}{\overline{U}}_{\eta \eta}^{-}+{\chi}^{\u2033}(0){Q}_{4}(({Q}_{2}-{Q}_{1}){\overline{U}}_{\xi}^{-}+{Q}_{3}{\overline{U}}_{\eta}^{-}),\\ {\overline{U}}_{\xi \xi}^{+}={\left(\frac{{c}^{-}}{{c}^{+}}\right)}^{2}{Q}_{1}{\overline{U}}_{\xi \xi}^{-}+\left({\left(\frac{{c}^{-}}{{c}^{+}}\right)}^{2}-1\right){Q}_{1}{\overline{U}}_{\eta \eta}^{-}\\ +{\chi}^{\u2033}(0)(({Q}_{2}-{Q}_{1}){\overline{U}}_{\xi}^{-}+{Q}_{3}{\overline{U}}_{\eta}^{-}),\end{array}$$

(12)

where *Q _{i}*,

$$\begin{array}{l}{Q}_{1}=\text{diag}(1,{\rho}^{-}/{\rho}^{+},1),\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}{Q}_{2}=\text{diag}({K}^{-}/{K}^{+},1,{\rho}^{+}/{\rho}^{-}),\\ {Q}_{3}=\frac{{\rho}^{-}}{{\rho}^{+}}\left({\left(\frac{{c}^{-}}{{c}^{+}}\right)}^{2}-1\right)\left(\begin{array}{lll}0\hfill & 1\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill \end{array}\right),\\ {Q}_{4}=\left(\begin{array}{lll}0\hfill & 1\hfill & 0\hfill \\ -1\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill \end{array}\right).\end{array}$$

(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

$${\mathrm{\Gamma}}_{ij,l}={Q}_{0}{\widehat{\mathrm{\Gamma}}}_{ij,l}{Q}_{0}^{-1},$$

where * _{ij,l}* are obtained by solving the following system of six matrix equations

$$\begin{array}{ll}\sum _{l=1}^{M}{\widehat{\mathrm{\Gamma}}}_{ij,l}{\widehat{Q}}_{1l}=0,\hfill & \sum _{l=1}^{M}{\widehat{\mathrm{\Gamma}}}_{ij,l}{\widehat{Q}}_{4l}=\frac{k}{h}{c}^{2}I-2\frac{{\xi}_{1}}{h}A,\hfill \\ \sum _{l=1}^{M}{\widehat{\mathrm{\Gamma}}}_{ij,l}{\widehat{Q}}_{2l}=-A,\hfill & \sum _{l=1}^{M}{\widehat{\mathrm{\Gamma}}}_{ij,l}{\widehat{Q}}_{5l}=-\frac{2}{h}({\eta}_{1}A+{\xi}_{1}B),\hfill \\ \sum _{l=1}^{M}{\widehat{\mathrm{\Gamma}}}_{ij,l}{\widehat{Q}}_{3l}=-B,\hfill & \sum _{l=1}^{M}{\widehat{\mathrm{\Gamma}}}_{ij,l}{\widehat{Q}}_{6l}=\frac{k}{h}{c}^{2}I-2\frac{{\eta}_{1}}{h}B.\hfill \end{array}$$

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

$$\begin{array}{l}{\widehat{Q}}_{1l}=I,\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}{\widehat{Q}}_{2l}=\frac{{\xi}_{l}}{h}I,\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}{\widehat{Q}}_{3l}=\frac{{\eta}_{l}}{h}I,\\ {\widehat{Q}}_{4l}={\left(\frac{{\xi}_{l}}{h}\right)}^{2}I,\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}{\widehat{Q}}_{5l}=2\frac{{\xi}_{l}{\eta}_{l}}{{h}^{2}}I,\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}{\widehat{Q}}_{6l}={\left(\frac{{\eta}_{l}}{h}\right)}^{2}I.\end{array}$$

(15)

Otherwise, we have

$$\begin{array}{l}{\widehat{Q}}_{1l}={Q}_{1},\\ {\widehat{Q}}_{2l}=\frac{{\xi}_{l}}{h}{Q}_{2}+\frac{1}{2h}{\chi}^{\u2033}(0)({\xi}_{l}^{2}-{\eta}_{l}^{2}+2{\xi}_{l}{\eta}_{l}{Q}_{4})({Q}_{2}-{Q}_{1}),\\ {\widehat{Q}}_{3l}=\frac{{\eta}_{l}}{h}{Q}_{1}+\frac{{\xi}_{l}}{h}{Q}_{3}+\frac{1}{2h}{\chi}^{\u2033}(0)({\xi}_{l}^{2}-{\eta}_{l}^{2}+2{\xi}_{l}{\eta}_{l}{Q}_{4}){Q}_{3},\\ {\widehat{Q}}_{4l}={\left(\frac{{\xi}_{l}{c}^{-}}{h{c}^{+}}\right)}^{2}{Q}_{1},\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}{\widehat{Q}}_{5l}=2\frac{{\xi}_{l}{\eta}_{l}}{{h}^{2}}{Q}_{2},\\ {\widehat{Q}}_{6l}={\left(\frac{{\eta}_{l}}{h}\right)}^{2}{Q}_{1}+2\frac{{\xi}_{l}{\eta}_{l}}{{h}^{2}}{Q}_{3}+{\left(\frac{{\xi}_{l}}{h}\right)}^{2}\left({\left(\frac{{c}^{-}}{{c}^{+}}\right)}^{2}-1\right){Q}_{1}.\end{array}$$

(16)

Here, (*ξ _{1}, η_{1}*) are the local coordinates of the

Finally, for completeness, let us also mention the two-dimensional *z*-transverse electric (TE* _{z}*) set of Maxwell’s equations,

$$\begin{array}{l}\frac{\partial {E}^{x}}{\partial t}-\frac{1}{\epsilon}\frac{\partial {H}^{z}}{\partial y}=0,\\ \frac{\partial {E}^{y}}{\partial t}+\frac{1}{\epsilon}\frac{\partial {H}^{z}}{\partial x}=0,\\ \frac{\partial {H}^{z}}{\partial t}-\frac{1}{\mu}\left(\frac{\partial {E}^{x}}{\partial y}-\frac{\partial {E}^{y}}{\partial x}\right)=0,\end{array}$$

(17)

where *H ^{z}* and

$$[\mathbf{t}\xb7\mathbf{E}]=0,\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}[\mathbf{n}\xb7\epsilon \mathbf{E}]=0,\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}[{H}^{z}]=0.$$

(18)

Likely, if we let

$$u={E}^{y},\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}v=-{E}^{x},\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}p={H}^{z},\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\rho =\epsilon ,\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}K=1/\mu ,$$

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, *u _{y}* −

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)

$$min\frac{1}{2}\sum _{l=1}^{M}{({\mathrm{\Gamma}}_{ij,l}-{D}_{l})}^{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

$$\begin{array}{l}{U}_{ij}^{n+1}={U}_{ij}^{n}+\frac{k}{h}[\left(\frac{k}{2h}{c}^{2}I+\frac{1}{2}A\right){U}_{i-1,j}^{n}\\ +\left(\frac{k}{2h}{c}^{2}I-\frac{1}{2}A\right){U}_{i+1,j}^{n}+\left(\frac{k}{2h}{c}^{2}I+\frac{1}{2}B\right){U}_{i,j-1}^{n}\\ +\left(\frac{k}{2h}{c}^{2}I-\frac{1}{2}B\right){U}_{i,j+1}^{n}-2\frac{k}{h}{c}^{2}{U}_{i,j}^{n}].\end{array}$$

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 *D _{l}* in the optimization problem (19) as

$$\begin{array}{l}{D}_{1}=-2\frac{k}{h}{c}^{2}I,\\ {D}_{2}=\frac{k}{2h}{c}^{2}I+\frac{1}{2}A,\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}{D}_{3}=\frac{k}{2h}{c}^{2}I-\frac{1}{2}A,\\ {D}_{4}=\frac{k}{2h}{c}^{2}I+\frac{1}{2}B,\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}{D}_{5}=\frac{k}{2h}{c}^{2}I-\frac{1}{2}B,\\ {D}_{l}=0,\phantom{\rule{0.38889em}{0ex}}l=6,\dots ,M.\end{array}$$

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

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

$$(A{U}_{x}+B{U}_{y}){\mid}_{ij}=A\left(\frac{{U}_{i+1,j}-{U}_{i-1,j}}{2h}\right)+B\left(\frac{{U}_{i,j+1}-{U}_{i,j-1}}{2h}\right).$$

(21)

On the other hand, at an irregular grid point (*x _{i}, y_{j}*) the difference equation

$$(A{U}_{x}+B{U}_{y}){\mid}_{ij}=\frac{1}{h}\sum _{l=1}^{M}{\mathrm{\Gamma}}_{ij,l}{U}_{ij,l}$$

(22)

is used, where the coefficients *Γ _{ij,l}* are 3 × 3 matrices, and

$$T=\frac{1}{h}\sum _{l=1}^{M}{\mathrm{\Gamma}}_{ij,l}{U}_{ij,l}-(A{U}_{x}({x}_{i},{y}_{j})+B{U}_{y}({x}_{i},{y}_{j}))+O(h).$$

(23)

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

$${\mathrm{\Gamma}}_{ij,l}={Q}_{0}{\widehat{\mathrm{\Gamma}}}_{ij,l}{Q}_{0}^{-1},$$

where in general the * _{ij,l}* matrices are obtained by solving the constrained quadratic optimization problem (19) subject to

$$\begin{array}{ll}\sum _{l=1}^{M}{\widehat{\mathrm{\Gamma}}}_{ij,l}{\widehat{Q}}_{1l}=0,\hfill & \sum _{l=1}^{M}{\widehat{\mathrm{\Gamma}}}_{ij,l}{\widehat{Q}}_{4l}=-2\frac{{\xi}_{1}}{h}A,\hfill \\ \sum _{l=1}^{M}{\widehat{\mathrm{\Gamma}}}_{ij,l}{\widehat{Q}}_{2l}=-A,\hfill & \sum _{l=1}^{M}{\widehat{\mathrm{\Gamma}}}_{ij,l}{\widehat{Q}}_{5l}=-\frac{2}{h}({\eta}_{1}A+{\xi}_{1}B),\hfill \\ \sum _{l=1}^{M}{\widehat{\mathrm{\Gamma}}}_{ij,l}{\widehat{Q}}_{3l}=-B,\hfill & \sum _{l=1}^{M}{\widehat{\mathrm{\Gamma}}}_{ij,l}{\widehat{Q}}_{6l}=-2\frac{{\eta}_{1}}{h}B.\hfill \end{array}$$

(24)

Here the 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 *D _{l}* matrices in the optimization problem (19) shall be chosen as

$$\begin{array}{l}{D}_{2}=\frac{1}{2}A,\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}{D}_{3}=-\frac{1}{2}A,\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}{D}_{4}=\frac{1}{2}B,\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}{D}_{5}=-\frac{1}{2}B,\\ {D}_{l}=0,\phantom{\rule{0.38889em}{0ex}}\text{otherwise}.\end{array}$$

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.

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

$${E}_{\text{inc}}^{z}={e}^{-\text{i}({k}_{1}x-\omega t)},\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}{H}_{\text{inc}}^{y}={e}^{-\text{i}({k}_{1}x-\omega t)},$$

where *k*_{1} = *ω*(*μ*_{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/(2\sqrt{2}{\lambda}_{max})$ and *k* = *h*/(10λ_{max}) for the Lax–Wendroff scheme-based IIM and the semi-discrete scheme-based one, respectively, where λ_{max} = max{*c*_{1}*, c*_{2}} and *c _{i}* = 1/(

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 *H ^{x}* at the time

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 *H*^{x} 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 (*x _{i}, y_{j}*), the sixth grid point in the 6-point stencil is chosen as the closest corner point to (

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 *H ^{x}* at the time

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 *H*^{x} 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 than10^{300} 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.

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
$\sqrt{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. 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

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

$${u}_{n+1}={u}_{n}+\mathrm{\Delta}tf\left({t}_{n}+\frac{\mathrm{\Delta}t}{2},{u}_{n}+\frac{\mathrm{\Delta}t}{2}f({t}_{n},{u}_{n})\right)$$

when solving the following initial value problem

$$\frac{\text{d}u}{\text{d}t}=f(t,u),\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}u({t}_{0})={u}_{0}.$$

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.

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.

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.

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.

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.

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |