## Abstract

Inverse method has wide application on medical diagnosis where non-destructive evaluation is the key factor .Back scattered waves or echoes generated from the forward moving waves has information about its geometry, size and location. In this paper we have investigated how well different geometries of the object is determined from the back scattered waves by a high accuracy Non-Standard Finite Difference Time Inverse (NSFD-TI) Maxwell’s algorithm and how the refractive index of the object plays a deterministic role on its size.

## 1. Introduction

In medical diagnosis, detection of malignant tumors of various size and shape from back scattered waves [1] with the application of the Inverse theorem is widespread. The challenge lies how well this method can detect the tumors when the refractive index of the malignant cells and the normal cells differs by a small quantity. In this paper we have investigated three different reconstructive spectrums of three different shapes of scatterer (circle, rectangle and triangle) using NSFD-TI algorithm based on Non-Standard Finite Difference [2,4] and also we have observed the effect of refractive index upon the size of the circular scatterer.

## 2. Model and Calculation Methods

Maxwell’s equations for a non-conducting dielectric are given by:

$∂tB(x,t)=−∇×E(x,t)$
$∂tD(x,t)=∇×H(x,t)−J(x,t)$
$∇·B(x,t)=0$
$∇·D(x,t)=ρ(x)$
where, H is the magnetic field, B is the magnetic induction, E the electric field, D the electric displacement, J the electric current density, and ρ is the density of free charge. Materials, such as iron or glass, modify the electric and magnetic fields. B describes how the material affects the magnetic field, and D describes how it affects electric field. The vector position is x = $(x,y,z)$, $∂t=∂/∂t$, and $∇=(∂x,∂y,∂z)$ is a vector differential operator. The curl of a vector is defined by$∇×V=(∂yVz−∂zVy,∂zVx−∂xVz,∂xVy−∂yVz)$, and the divergence of V is$∇·V=∂xVx+∂yVy+∂zVz$, where$V=(Vx,Vy,Vz)$.

In most optical materials there is no free charge and$ρ=0$. In a non-conducting dielectric, such as glass, no electric current flow so J = 0. Also to an excellent approximation

$B(x,t)=μH(x,t)$
where, μ, the magnetic permeability, is a constant. In general, the electrical permittivity, ε depends upon the frequency, but when this dependence is weak, it is a good approximation to take
$D(x,t)=ε(x)E(x,t)$
The permittivity is different for different materials, however so it is a function of position. For example, in glass$ε≅1.2$. Assuming that Eq. (5) and Eq. (6) hold, Maxwell’s equations for a non-magnetic, non-conducting medium with no free charges or electrical currents is
$μ∂tH(x,t)=−∇×E(x,t)$
$ε(x)∂tE(x,t)=∇×H(x,t)$
Now, let us introduce what is called a Finite Difference (FD) approximation. The first derivative is approximately given by
$f′(x)≅f(x+h2)−f(x−h2)h$
For convenience let us define the difference operator
$dxf(x)=f(x+h2)−f(x−h2)$
with this definition the FD approximation to $f′$ is
$f′(x)≅dxhf(x)$
Defining the vector difference operator:
$d=(dx,dy,dz)$
and,
$∇f(x)≅dhf(x)$
where, $Δx=Δy=Δz=h$. The second derivative$f″$ can be approximated by applying Eq. (11) twice:
$f″(x)≅dx2h2f(x).$
and from the definition of $dx$ it is easy to show that $dx2=dxdx$ is given by:
$dx2f(x)=f(x+h)+f(x−h)−2f(x)$
Defining,
$d2=dx2+dy2+dz2$
We now have,
$∇2f(x)≅d2h2f(x)$
and also replacing the derivatives of Maxwell’s Eq. (7) by FD approximations we obtain
$dtH(x,t)=−1μΔthd×E(x,t)$
Here the electromagnetic field components are arranged on the grid in such a way those central FDs can be taken. Expanding$dtH(x,t)$, and solving for $H(x,t+Δt2)$ gives
$H(x,t+Δt2)=H(x,t−Δt2)−1μΔthd×E(x,t)$
Doing the same for Eq. (8) we have
$dtE(x,t+Δt2)=−1ε(x)Δthd×H(x,t+Δt2)$
and we compute$dtE(x,t+Δt2)$, rather than $dtE(x,t)$ in order to use the updated value of H. Expanding$dtE(x,t+Δt2)$, and solving for $E(x,t+Δt)$ gives
$E(x,t+Δt)=E(x,t)−1ε(x)Δthd×H(x,t+Δt2)$
Equation (19) and Eq. (21) are called the Yee algorithm [3]. The notation does not show it, but each component of H and E is evaluated at a different position. For example to compute$dtHz∼dxEy−dyEx$, if $Hz$ lies on $(x,y,z)$ then $Ey$ must be known at$(x±h2,y,z)$, and $Ex$ at$(x,y±h2,z)$.

To develop the Nonstandard-FDTD (NS-FDTD) model or in short NSFD we first construct an alternative Standard-FDTD (S-FDTD) model or SFD of Maxwell’s equations. Defining$H′=μhH/Δt$, we obtain

$dtH′(x,t)=−d1×E(x,t),$
$dtE(x,t+Δt2)=1εμΔt2h2d1×H′(x,t+Δt2)$
We introduced in earlier publication [4,7], the vector difference operator $d0$ which is also highly isotropic in nature with the property that
$d1·d0=d0·d1=d02.$
whereas, $d1·d1=d02$,$d0·d0≠d02$. In two dimensions$d0=(dx(0),dy(0))$, where
$dx(0)=dx+(1−γ4)dxdy2$
$dy(0)=dy+(1−γ4)dydx2$
$andd0=d1(1+1−γ4dxdy)$
Here, γ which is also a function of $k=k0ε$ (assuming uniform μ), given by the expression $γ=23−190(kh)2−115120(kh)4(11−52)−⋯$ is chosen to minimize the error of the NSFD approximation [4].

Making the substitutions $d1→d0$ and $Δt2/(h2εiμ)→u0(ε)2$ in Eq. (23) we obtain a new NSFD model of Maxwell’s equations,

$dtH′(x,t)=−d1×E(x,t)$
$dtE(x,t+Δt2)=u02εd0(ε)×H′(x,t+Δt2)$
Solving Eq. (28) and Eq. (29) for $H′(x,t+Δt2)$ and$E(x,t+Δt)$, we obtain the NS-Yee algorithm
$H′(x,t+Δt2)=H′(x,t−Δt2)−d1×E(x,t)$
$E(x,t+Δt)=E(x,t)+u02εd0(ε)×H′(x,t+Δt2)$
This NSFD algorithm is somewhat different from the one introduced in [4]. We find that it gives better results when μ is constant. To reverse the time direction we need to reverse the sequence of the above computations; doing so we get,
$H′(x,t−Δt2)=H′(x,t+Δt2)+d1×E(x,t+Δt)$
and
$E(x,t)=E(x,t+Δt)−u02εd0(ε)×H′(x,t−Δt2)$
Here, we observed that we basically get the same equations as we already derived for forward propagating wave. This depicts that even if we iterate the wave in time reverse direction the nature of the wave remain same. For this purpose we have applied the periodic boundary conditions instead of Mur boundary conditions because the later causes field explosion when we iterate in time reverse condition [5]. Let, the computational domain be$x=0,h,2h,⋯Nxh$. In periodic boundary we define, $ψ(−h,t)=ψ(Nxh,t)$where, waves that exit the computational domain at $x=0$ re-enter it from$x=Nxh$.

## Algorithm Stability

When any algorithm is iterated is can be become unstable. Let A be an algorithm to solve for an unknown vector ψ. Let $φ0$ be the initial solution vector. Then$φ1=Aφ0$, $φ2=Aφ1$, ⋯, and thus $φn=Anφ0$. If $lim|φn|n→∞→∞$ when $|ψ|$ is finite, the algorithm is said to be unstable.

It can be shown that the D-dimensional FDTD algorithm is stable if$vΔth≤DD$. Taking finite-difference algorithm of the form:

$ψ(t+Δt)=pψ(t−Δt)+2qψ(t)$
Defining $t=nΔt$and $ψn=ψ(nΔt)$we get,
$ψn+1=pψn−1+2qψn$
Now we postulate a solution:Inserting $ψn=λn$ and dividing by$λn−1$ we get,
$λ2−2qλ−p=0$
Solving for λ we get,
$λ±=q±q2+p$
Thus the most general solution is:
$ψn=c+λ+n+c−λ−n$
The condition for the stability will be$|ψn|$ being finite as $n→∞$ is$|λ±|≤1$.

We want finite oscillatory solutions of the difference Eq. (36), for the oscillation the solutions must have an imaginary part. If $p<0⇒p=−|p|$ and $p+q2<0$ and $λ±=q±q2−|p|=q±i|p|−q2⇒|λ±|2=q2−|p|+q2=|p|$ . Thus $|p|≤1$ and $q2<|p|$ are the conditions for stability for an oscillating solution.

Applying these results to the wave equation given below:

$(∂t2−v2∇2)ψ(x,t)=0,$
Generalized Finite Difference model is:
$(dt2−u02d2)ψ(x,t)=0$
and the FDTD algorithm is
$ψ(x,t+Δt)=−ψ(x,t−Δt)+(2+u02d2)ψ(x,t)$
While in S-FDTD algorithm, $d2=d12$ and$u02=v2Δt2h2$, whereas in NS-FDTD $d2=d02$ and$u02=sin2(ωΔt/2)sin2(kh/2)$. To direct the FDTD algorithm in the form of (34) we need to compute$d2ψ$. Substituting $ψ=φk=ei(k.x−ωt)$ we can write, $d2φk=−2d2φk(x,t)$ where the values of $d2$ vary case by case [4]:

For SFD:

Max$(d12)=2D$where, D = Dimension

For NSFD:

also,
Till now we have not specified the value ofγ. For the stability we optimize the value of γ by $γ0$ (also termed as weight coefficient for$d12$and$d22$)Where
$γ0=1−sin2(k′xh/2)+sin2(k′yh/2)−sin2(kh/2)2sin2(k′xh/2)sin2(k′yh/2)$
and$(k′x,k′y)=k(2−1/4,1−2−1/2)$.Thus for the NSFD algorithm we have,
$d02=γ0d12+(1−γ0)d22$
From the expanded value of γ given above, we get Max ($γ0$) = 2/3. The above equation can be rewritten as:
$d02(k)=γ0d12(k)+(1−γ0)d22(k)$
where, $k=k(cosθ,sinθ)$and Max$(d02)=8/3$ in 2 Dimensions.If we write, $ψ(x,t)=ψxt$ then the FDTD algorithm (41) becomes:
$φxt+1=−φxt−1+2(1−d2u02)φxt$
This equation is in the form of (36) and the solutions are:
$λ±=1−u02d2±(1−u02d2)2−1$
To accomplish the stability condition$|λ±|≤1$, we have to have$(1−u02d2)2≤1$.

If $d2>0$and$u02>0$, for non-vanishing monochromatic plane waves, the stability condition turns out to be [4]:

$u02≤2d2$
While$d2$is a function of the magnitude and direction of k, stability is guaranteed by inserting the maximum value of$d2$, Max$(d2)$ instead of$d2$.
$u02≤2Max(d2)$
Considering NSFD algorithm, we have $d02=d12=d2$and Max$(d02)=2,8/3$for dimensions, D = 1, 2 respectively. Defining earlier, $u02=sin2(ωΔt/2)sin2(kh/2)$the stability condition (48) becomes,
$sin(ωΔt/2)sin(kh/2)≤2Max(d02)$
Replacing, $ϖ=ωΔt,k¯=kh$and$ϖ=ck¯$, where the value of c is to be found out so that (49) does not violate. In D dimensions expressing $Max(d02)=mdnsfd$ the Eq. (49) turns out to be:
$sin(ck¯/2)≤2mdnsfdsin(k¯/2)$
In the range$0≤θ≤π/2$, $sin(θ)$rises monotonically with increasing θ, so it is obvious that$sinθ1≤sinθ2⇒θ1≤θ2$. The Nyquist limit can be articulated in the form$k<π/D$. Hence we shall never figure out a FDTD algorithm which disobeys $λh>2D$ condition, and presuming that, $c≤1$ we conclude from (50) that
$c<2dπarcsin[2mdnsfdsin(π2D)]$
Putting the numerical values of $mdnsfd$from above, we have

## 3. Verifications and Practical Tests

In Fig. 1 we compare FDTD calculations of Mie scattering [6] for an infinite dielectric cylinder of radius $0.65λ0$ and refractive index$n=1.8$. The cylinder is parallel to the z-axis. An infinite plane wave propagating in the $+x$ direction impinges upon the cylinder. The electric field is polarized parallel to the cylinder axis, and the scattered intensity is visualized in shades of red for the analytic solution (left) the NS-FDTD algorithm (middle) and the S-FDTD algorithm (right) with a discretization of $λ0/h=8$ outside the cylinder, and$λ0/nh≅4.4$ inside. Similar result was published by one of the author in his previous publication [7].

Fig. 1 Mie scattering from an infinite dielectric cylinder. Scattered intensity is visualized in shades of red. E is polarized parallel cylinder axis. Infinite plane wave impinges from the left. Cylinder radius = $0.65λ0$, $λ0$ = vacuum wavelength, refractive index $n=1.8$.

In Fig. 2(a) the scattered fields shown in Fig. 1 are plotted as a function of angle the (θ) from the $+x$-axis on a circular contour of radius $λ0$ centered on the axis of the cylinder. At discretization $λ0/h=8$ the root mean square error for the NS-FDTD and S-FDTD algorithms is $εNS-8=0.04$, and $εS-8=0.20$, respectively. In Fig. 2(b) at discretization $λ0/h=24$ we compare the S-FDTD calculation with the analytic solution. The root mean square error is $εS-24=0.04$ the same as the NS-FDTD algorithm at$λ0/h=8$. At $λ0/h=8$ the NS-FDTD algorithm delivers the same accuracy as the S-FDTD one does at $λ0/h=24$. Because of the lower coarser discretization, the computational cost of NSFDTD algorithm is only $1/27$th that of the Standard FDTD one.

Fig. 2 Angular distribution of scattered intensity about a circular contour of radius $λ0$ centered on the axis of the cylinder. (a) Analytic solution (black) compared with the NSFD calculation (NSFD-8, red) and the SFD one (SFD-8, blue) at $λ0/h=8$. (b) Analytic solution (black) compared with the SFD (SFD-24, blue) one at $λ0/h=24$.

Understanding the fact that NSFD algorithm is better than S-FDTD one, we extended our study for the Time Inverse algorithm on the detection of size, shape of different scatterers. for this purpose we have generated a weighted source array continuous wave Gaussian beam that turns on at the time t = 1. The beam width of the Gaussian pulse taken for the simulation purpose is 1.25$λ0$, where $λ0$ is the wavelength of the central frequency of the pulse. For the simulation we have chosen three types of the scatterer; Circular, Rectangular and triangular. Now we have allowed the light to be passed through the scatterer and it further moves to a distance which is very close to the periodic boundary of the computational space ($Nx=$140h), then we remove the scatterer and iterate the forward moving wave in the time reverse in order to trace the location and shape of the scatterer and examine how well it will match with the actual location and size. For this experiment, we have confined our observation on TE mode only, i.e. we will compute x-intensity and y-intensity as depicted in Fig. 3 .

Fig. 3 Time Inverse Computed Intensities using NS-FDTD for: (a)Circular scatterer (diameter = 2$λ0$)(b) Rectangular scatterer (side = 2$λ0$) (c)Triangular scatterer(base = 2$λ0$,height = 1$λ0$),where, $λ0/h=10$

In the next experiment, our objective is to investigate the effect of refractive index on the size of the scatterer. The minimum detectable size in case is called distinguishable size. The refractive index of the scatterer is fixed and in all the cases is taken to 1.8. In this experiment, we have considered circular scatterer because it is homogeneous in shape and no steep edges and corners. Now, gradually we have changed the refractive index of the circular scatterer from 1.1 to 1.8 and compute the distinguishable size of the circular scatterer for a particular refractive index and plot the graph shown in Fig. 4 for TM Mode of polarization.

Fig. 4 Distinguishable Size vs. Refractive Index

## 4. Results and discussions

From the Fig. 3 we noticed that the Time Inverse algorithm (TE mode) has two intensities, x-intensity and y intensity in which left one is the computed x-intensity and right one is the computed y-intensity. It is clearly seen from the figure that x-intensities gives the clear indication of the shape, which means the first one is circular, the middle one is the rectangular and the last one is the triangular scatterer. The y- intensity is perpendicularly oriented to the x-intensity and it does not differentiate the edges well as it does for the x-intensity. But considering both the intensities we can distinguish the exact shape, size and location of the scatterer by Non standard Time Domain Time Inverse algorithm. It is very helpful in medical imaging because a small cancerous tumor of arbitrary shape is perfectly detected and chance of misdetection or false detection is eliminated.

In Fig. 4, shown below applying the TM Mode of polarization of incoming light we have noticed that as the refractive index of the circular scatterer gradually increases the distinguishable size of the scatterer gradually decreases. One of the probable reasons is that, the scatterer having the refractive index close to vacuum, it size must be large enough in order to get some considerable contrast. We have also noticed that our Non-Standard Time Domain Inverse algorithm is good enough to detect the shape, size and location of the object in comparison to Standard Time Domain Inverse algorithm even if the refractive index is very close to unity. In the later case more time steps are required and the picture contrast is not good as that of NSFD algorithm.

## Acknowledgements

We would like to acknowledge the Ministry of Education, Government of Japan for granting one of the authors (KC) with a scholarship during the period of this work.

1. L. W. Schmerr, Jr., in Fundamentals of Ultrasonic Nondestructive Evaluation: A Modeling Approach, (Plenum Press, New York, 1998).

2. R. Mickens, E, in Nonstandard finite difference models of differential equations, (World Scientific Publishing Co., Inc., River Edge, NJ, 1994).

3. K. S. Kunz, and R. J. Luebbers, in The Finite Difference Time Domain Method for Electromagnetics, (CRC Press, New York, 1993).

4. J. B. Cole, S. Banerjee, and M. Haftel, in Advances in the Applications of Nonstandard Finite Difference Schemes, ed. R. E. Mickens (World Scientific Singapore, 2007), Chap.4.

5. R. Sorrentino, L. Roselli, and P. Mezzanotte, “Time reversal in finite difference time domain method,” IEEE Microw. Guid. Wave Lett. 3(11), 402–404 (1993). [CrossRef]

6. P. W. Barber, and S. C. Hill, in Light Scattering by Particles: Computational Methods, (World Scientific, Singapore, 1990).

7. J. B. Cole, “High-Accuracy Yee algorithm Based on Nonstandard Finite Difference: New Developments and Verifications,” IEEE Trans. Antenn. Propag. 50(9), 1185–1191 (2002). [CrossRef]

8. H. Kudo, et al., “Numerical Dispersion and Stability Condition of the Nonstandard FDTD Method” Electronics and Communications in Japan, 85, 22–30(2002), http://www3.interscience.wiley.com/cgi-bin/fulltext/93514073/PDFSTART.

### References

• View by:
• |
• |
• |

1. L. W. Schmerr, Jr., in Fundamentals of Ultrasonic Nondestructive Evaluation: A Modeling Approach, (Plenum Press, New York, 1998).
2. R. Mickens, E, in Nonstandard finite difference models of differential equations, (World Scientific Publishing Co., Inc., River Edge, NJ, 1994).
3. K. S. Kunz, and R. J. Luebbers, in The Finite Difference Time Domain Method for Electromagnetics, (CRC Press, New York, 1993).
4. J. B. Cole, S. Banerjee, and M. Haftel, in Advances in the Applications of Nonstandard Finite Difference Schemes, ed. R. E. Mickens (World Scientific Singapore, 2007), Chap.4.
5. R. Sorrentino, L. Roselli, and P. Mezzanotte, “Time reversal in finite difference time domain method,” IEEE Microw. Guid. Wave Lett. 3(11), 402–404 (1993).
[CrossRef]
6. P. W. Barber, and S. C. Hill, in Light Scattering by Particles: Computational Methods, (World Scientific, Singapore, 1990).
7. J. B. Cole, “High-Accuracy Yee algorithm Based on Nonstandard Finite Difference: New Developments and Verifications,” IEEE Trans. Antenn. Propag. 50(9), 1185–1191 (2002).
[CrossRef]
8. H. Kudo, et al., “Numerical Dispersion and Stability Condition of the Nonstandard FDTD Method” Electronics and Communications in Japan, 85, 22–30(2002), http://www3.interscience.wiley.com/cgi-bin/fulltext/93514073/PDFSTART .

#### 2002 (1)

J. B. Cole, “High-Accuracy Yee algorithm Based on Nonstandard Finite Difference: New Developments and Verifications,” IEEE Trans. Antenn. Propag. 50(9), 1185–1191 (2002).
[CrossRef]

#### 1993 (1)

R. Sorrentino, L. Roselli, and P. Mezzanotte, “Time reversal in finite difference time domain method,” IEEE Microw. Guid. Wave Lett. 3(11), 402–404 (1993).
[CrossRef]

#### Cole, J. B.

J. B. Cole, “High-Accuracy Yee algorithm Based on Nonstandard Finite Difference: New Developments and Verifications,” IEEE Trans. Antenn. Propag. 50(9), 1185–1191 (2002).
[CrossRef]

#### Mezzanotte, P.

R. Sorrentino, L. Roselli, and P. Mezzanotte, “Time reversal in finite difference time domain method,” IEEE Microw. Guid. Wave Lett. 3(11), 402–404 (1993).
[CrossRef]

#### Roselli, L.

R. Sorrentino, L. Roselli, and P. Mezzanotte, “Time reversal in finite difference time domain method,” IEEE Microw. Guid. Wave Lett. 3(11), 402–404 (1993).
[CrossRef]

#### Sorrentino, R.

R. Sorrentino, L. Roselli, and P. Mezzanotte, “Time reversal in finite difference time domain method,” IEEE Microw. Guid. Wave Lett. 3(11), 402–404 (1993).
[CrossRef]

#### IEEE Microw. Guid. Wave Lett. (1)

R. Sorrentino, L. Roselli, and P. Mezzanotte, “Time reversal in finite difference time domain method,” IEEE Microw. Guid. Wave Lett. 3(11), 402–404 (1993).
[CrossRef]

#### IEEE Trans. Antenn. Propag. (1)

J. B. Cole, “High-Accuracy Yee algorithm Based on Nonstandard Finite Difference: New Developments and Verifications,” IEEE Trans. Antenn. Propag. 50(9), 1185–1191 (2002).
[CrossRef]

#### Other (6)

H. Kudo, et al., “Numerical Dispersion and Stability Condition of the Nonstandard FDTD Method” Electronics and Communications in Japan, 85, 22–30(2002), http://www3.interscience.wiley.com/cgi-bin/fulltext/93514073/PDFSTART .

P. W. Barber, and S. C. Hill, in Light Scattering by Particles: Computational Methods, (World Scientific, Singapore, 1990).

L. W. Schmerr, Jr., in Fundamentals of Ultrasonic Nondestructive Evaluation: A Modeling Approach, (Plenum Press, New York, 1998).

R. Mickens, E, in Nonstandard finite difference models of differential equations, (World Scientific Publishing Co., Inc., River Edge, NJ, 1994).

K. S. Kunz, and R. J. Luebbers, in The Finite Difference Time Domain Method for Electromagnetics, (CRC Press, New York, 1993).

J. B. Cole, S. Banerjee, and M. Haftel, in Advances in the Applications of Nonstandard Finite Difference Schemes, ed. R. E. Mickens (World Scientific Singapore, 2007), Chap.4.

### Cited By

OSA participates in CrossRef's Cited-By Linking service. Citing articles from OSA journals and other participating publishers are listed here.

### Figures (4)

Fig. 1

Mie scattering from an infinite dielectric cylinder. Scattered intensity is visualized in shades of red. E is polarized parallel cylinder axis. Infinite plane wave impinges from the left. Cylinder radius = $0.65 λ 0$ , $λ 0$ = vacuum wavelength, refractive index $n = 1.8$ .

Fig. 2

Angular distribution of scattered intensity about a circular contour of radius $λ 0$ centered on the axis of the cylinder. (a) Analytic solution (black) compared with the NSFD calculation (NSFD-8, red) and the SFD one (SFD-8, blue) at $λ 0 / h = 8$ . (b) Analytic solution (black) compared with the SFD (SFD-24, blue) one at $λ 0 / h = 24$ .

Fig. 3

Time Inverse Computed Intensities using NS-FDTD for: (a)Circular scatterer (diameter = 2 $λ 0$ )(b) Rectangular scatterer (side = 2 $λ 0$ ) (c)Triangular scatterer(base = 2 $λ 0$ ,height = 1 $λ 0$ ),where, $λ 0 / h = 10$

Fig. 4

Distinguishable Size vs. Refractive Index

### Equations (57)

$∂ t B ( x , t ) = − ∇ × E ( x , t )$
$∂ t D ( x , t ) = ∇ × H ( x , t ) − J ( x , t )$
$∇ · B ( x , t ) = 0$
$∇ · D ( x , t ) = ρ ( x )$
$B ( x , t ) = μ H ( x , t )$
$D ( x , t ) = ε ( x ) E ( x , t )$
$μ ∂ t H ( x , t ) = − ∇ × E ( x , t )$
$ε ( x ) ∂ t E ( x , t ) = ∇ × H ( x , t )$
$f ′ ( x ) ≅ f ( x + h 2 ) − f ( x − h 2 ) h$
$d x f ( x ) = f ( x + h 2 ) − f ( x − h 2 )$
$f ′ ( x ) ≅ d x h f ( x )$
$d = ( d x , d y , d z )$
$∇ f ( x ) ≅ d h f ( x )$
$f ″ ( x ) ≅ d x 2 h 2 f ( x ) .$
$d x 2 f ( x ) = f ( x + h ) + f ( x − h ) − 2 f ( x )$
$d 2 = d x 2 + d y 2 + d z 2$
$∇ 2 f ( x ) ≅ d 2 h 2 f ( x )$
$d t H ( x , t ) = − 1 μ Δ t h d × E ( x , t )$
$H ( x , t + Δ t 2 ) = H ( x , t − Δ t 2 ) − 1 μ Δ t h d × E ( x , t )$
$d t E ( x , t + Δ t 2 ) = − 1 ε ( x ) Δ t h d × H ( x , t + Δ t 2 )$
$E ( x , t + Δ t ) = E ( x , t ) − 1 ε ( x ) Δ t h d × H ( x , t + Δ t 2 )$
$d t H ′ ( x , t ) = − d 1 × E ( x , t ) ,$
$d t E ( x , t + Δ t 2 ) = 1 ε μ Δ t 2 h 2 d 1 × H ′ ( x , t + Δ t 2 )$
$d 1 · d 0 = d 0 · d 1 = d 0 2 .$
$d x ( 0 ) = d x + ( 1 − γ 4 ) d x d y 2$
$d y ( 0 ) = d y + ( 1 − γ 4 ) d y d x 2$
$and d 0 = d 1 ( 1 + 1 − γ 4 d x d y )$
$d t H ′ ( x , t ) = − d 1 × E ( x , t )$
$d t E ( x , t + Δ t 2 ) = u 0 2 ε d 0 ( ε ) × H ′ ( x , t + Δ t 2 )$
$H ′ ( x , t + Δ t 2 ) = H ′ ( x , t − Δ t 2 ) − d 1 × E ( x , t )$
$E ( x , t + Δ t ) = E ( x , t ) + u 0 2 ε d 0 ( ε ) × H ′ ( x , t + Δ t 2 )$
$H ′ ( x , t − Δ t 2 ) = H ′ ( x , t + Δ t 2 ) + d 1 × E ( x , t + Δ t )$
$E ( x , t ) = E ( x , t + Δ t ) − u 0 2 ε d 0 ( ε ) × H ′ ( x , t − Δ t 2 )$
$ψ ( t + Δ t ) = p ψ ( t − Δ t ) + 2 q ψ ( t )$
$ψ n + 1 = p ψ n − 1 + 2 q ψ n$
$λ 2 − 2 q λ − p = 0$
$λ ± = q ± q 2 + p$
$ψ n = c + λ + n + c − λ − n$
$( ∂ t 2 − v 2 ∇ 2 ) ψ ( x , t ) = 0 ,$
$( d t 2 − u 0 2 d 2 ) ψ ( x , t ) = 0$
$ψ ( x , t + Δ t ) = − ψ ( x , t − Δ t ) + ( 2 + u 0 2 d 2 ) ψ ( x , t )$
$γ 0 = 1 − sin 2 ( k ′ x h / 2 ) + sin 2 ( k ′ y h / 2 ) − sin 2 ( k h / 2 ) 2 sin 2 ( k ′ x h / 2 ) sin 2 ( k ′ y h / 2 )$
$d 0 2 = γ 0 d 1 2 + ( 1 − γ 0 ) d 2 2$
$d 0 2 ( k ) = γ 0 d 1 2 ( k ) + ( 1 − γ 0 ) d 2 2 ( k )$
$φ x t + 1 = − φ x t − 1 + 2 ( 1 − d 2 u 0 2 ) φ x t$
$λ ± = 1 − u 0 2 d 2 ± ( 1 − u 0 2 d 2 ) 2 − 1$
$u 0 2 ≤ 2 d 2$
$u 0 2 ≤ 2 Max ( d 2 )$
$sin ( ω Δ t / 2 ) sin ( k h / 2 ) ≤ 2 Max( d 0 2 )$
$sin ( c k ¯ / 2 ) ≤ 2 m d n s f d sin ( k ¯ / 2 )$
$c < 2 d π arcsin [ 2 m d n s f d sin ( π 2 D ) ]$