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.

© 2010 OSA

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=(yVzzVy,zVxxVz,xVyyVz), and the divergence of V is·V=xVx+yVy+zVz, whereV=(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(xh2)h
For convenience let us define the difference operator
dxf(x)=f(x+h2)f(xh2)
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 derivativef 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(xh)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. ExpandingdtH(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 computedtE(x,t+Δt2), rather than dtE(x,t) in order to use the updated value of H. ExpandingdtE(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 computedtHzdxEydyEx, 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. DefiningH=μ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·d0d02. In two dimensionsd0=(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 γ=23190(kh)2115120(kh)4(1152) is chosen to minimize the error of the NSFD approximation [4].

Making the substitutions d1d0 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) andE(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 bex=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 fromx=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 ifvΔthDD. Taking finite-difference algorithm of the form:

ψ(t+Δt)=pψ(tΔt)+2qψ(t)
Defining t=nΔtand ψn=ψ(nΔt)we get,
ψn+1=pψn1+2qψn
Now we postulate a solution:Inserting ψn=λn and dividing byλn1 we get,
λ22qλ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<0p=|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:

(t2v22)ψ(x,t)=0,
Generalized Finite Difference model is:
(dt2u02d2)ψ(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 andu02=v2Δt2h2, whereas in NS-FDTD d2=d02 andu02=sin2(ωΔt/2)sin2(kh/2). To direct the FDTD algorithm in the form of (34) we need to computed2ψ. 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:

d2=d12=2sin2(kh/2): 1Dimension
d2=d12=2[sin2(kxh/2)+sin2(kyh/2)]: 2Dimensions
Max(d12)=2Dwhere, D = Dimension

For NSFD:

d02=d12=dx2and Max(d02)=2: 1Dimension
also,
d22=1cos(kxh/2)cos(kxh/2): 2Dimensions
Till now we have not specified the value ofγ. For the stability we optimize the value of γ by γ0 (also termed as weight coefficient ford12andd22)Where
γ0=1sin2(kxh/2)+sin2(kyh/2)sin2(kh/2)2sin2(kxh/2)sin2(kyh/2)
and(kx,ky)=k(21/4,121/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=φxt1+2(1d2u02)φxt
This equation is in the form of (36) and the solutions are:
λ±=1u02d2±(1u02d2)21
To accomplish the stability condition|λ±|1, we have to have(1u02d2)21.

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

u022d2
Whiled2is a function of the magnitude and direction of k, stability is guaranteed by inserting the maximum value ofd2, Max(d2) instead ofd2.
u022Max(d2)
Considering NSFD algorithm, we have d02=d12=d2and Max(d02)=2,8/3for 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¯=khandϖ=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 range0θπ/2, sin(θ)rises monotonically with increasing θ, so it is obvious thatsinθ1sinθ2θ1θ2. The Nyquist limit can be articulated in the formk<π/D. Hence we shall never figure out a FDTD algorithm which disobeys λh>2D condition, and presuming that, c1 we conclude from (50) that
c<2dπarcsin[2mdnsfdsin(π2D)]
Putting the numerical values of mdnsfdfrom above, we have

c1nsfd=1: In One Dimension
c1nsfd22πarcsin[32sin(π22)]0.84: In Two Dimensions

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 indexn=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/nh4.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.

Download Full Size | PPT Slide | PDF

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/27th 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.

Download Full Size | PPT Slide | PDF

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

Download Full Size | PPT Slide | PDF

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

Download Full Size | PPT Slide | PDF

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.

References and links

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.

Alert me when this article is cited.


Figures (4)

Fig. 1
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
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
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
Fig. 4

Distinguishable Size vs. Refractive Index

Equations (57)

Equations on this page are rendered with MathJax. Learn more.

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 )
d 2 = d 1 2 = 2 sin 2 ( k h / 2 ) : 1Dimension
d 2 = d 1 2 = 2 [ sin 2 ( k x h / 2 ) + sin 2 ( k y h / 2 ) ] : 2Dimensions
d 0 2 = d 1 2 = d x 2 and Max ( d 0 2 ) = 2 : 1Dimension
d 2 2 = 1 cos ( k x h / 2 ) cos ( k x h / 2 ) : 2Dimensions
γ 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 ) ]
c 1 nsfd = 1 : In One Dimension
c 1 nsfd 2 2 π arcsin [ 3 2 sin ( π 2 2 ) ] 0.84 : In Two Dimensions

Metrics