World Academy of Science, Engineering and Technology

International Journal of Computer and Information Engineering

Vol:6, No:3, 2012

Simulation of the Finite Difference Time

Domain in Two Dimension
Akram G., Jasmy Y.

∂H y ∂ ( Ezx + Ezy )
Abstract—The finite-difference time-domain (FDTD) method is µ + σ x* H y = (2)
one of the most widely used computational methods in ∂t ∂x
electromagnetic. This paper describes the design of two-dimensional
(2D) FDTD simulation software for transverse magnetic (TM)
∂Ezx ∂H y
polarization using Berenger's split-field perfectly matched layer
ε + σ x Ezx = (3)
∂t ∂x
International Science Index, Computer and Information Engineering Vol:6, No:3, 2012 waset.org/Publication/5827

(PML) formulation. The software is developed using Matlab

programming language. Numerical examples validate the software.

Keywords—Finite difference time domain (FDTD) method, ∂Ezy ∂H x

ε + σ y Ezy = − (4)
perfectly matched layer (PML), split-filed formulation, transverse
magnetic (TM) polarization.
∂t ∂y

I. INTRODUCTION In order to solve these equations numerically, they are

discretized according to finite differencing techniques. The
T HE finite-difference time-domain (FDTD) method is a
flexible and powerful solution tool for solving Maxwell’s
equations [1-5]. Further, the efficient and accurate solution of
central difference approximation is used in this case. For
example, in equations (1) and (3), central differencing results
in the following equations
electromagnetic wave interaction problems in unbounded
regions is one of the greatest challenges of the FDTD method. −σ *y ( i +1 2, j )δ t µ
For such problems, an absorbing boundary condition (ABC) H xn +1 ( i + 1 2, j ) = e H xn ( i + 1 2, j )
(5)
must be introduced at the outer grid boundary to simulate the
extension of the grid to infinity. −
(1 − e −σ *y ( i +1 2, j )δ t µ
)  E
( i + 1 2, j + 1 2 ) + Ezyn+1 2 ( i + 1 2, j + 1 2 ) 
n +1 2
zx

σ y ( i + 1 2, j ) δ x  − Ezx ( i + 1 2, j − 1 2 ) − Ezyn +1 2 ( i + 1 2, j − 1 2 ) 
* n +1 2
Recently, among a few important studies on ABC was J.P.
Berenger’s Perfectly Matched Layer (PML) method [6]. In the
E zxn +1 2 ( i + 1 2 , j + 1 2 ) = e E zxn −1 2 ( i + 1 2 , j + 1 2 )
−σ *x ( i +1 2, j +1 2 )δ t ε
PML, he created an artificial medium with magnetic

(1 − e ) H
conductivity (σ*) that will make the boundary condition to −σ *x ( i +1 2, j +1 2 )δ t ε
(6)
work as wave absorber region. This paper describes the design − n
( i, j + 1 2 ) − H yn ( i + 1, j + 1 2 ) 
of a two-dimensional (2D) FDTD simulation software for σ x* ( i + 1 2 , j + 1 2 ) δ x  y

transverse magnetic (TM) polarized incident wave using

Berenger's split-field perfectly matched layer (PML) Based on the above formulation, a 2-D FDTD code has
formulation. The software is developed using Matlab been written in Matlab programming language to simulate the
programming language. fields of a plane wave source in lossy media. The formulation
and the FDTD code are validated by checking the numerical
II. TWO-DIMENSIONAL TM POLARIZATION results for homogeneous media against the analytical solution.
The Maxwell’s curl equations as modified by Berenger are The code for 2D FDTD is shown in Appendix.
expressed in their time-dependent form as, In order to test and validate the FDTD code, an
electromagnetic wave source interacting with dielectric

∂H x ∂ ( Ezx + Ezy ) cylinder (εr = 4.0, σ = 0.12) of radius 6 cm is simulated in this

µ + σ *y H x = − (1) section, shown in Fig. 1. A 1 mW/cm2 plane wave is used as
∂t ∂y source, and the frequency is setup to 2.5 GHz. The reason for
using dielectric cylinder is there exist an analytical solution to
this problem [7]. A +y-directed plane wave source :

G. Akram is with the Faculty of Health Science & Biomedical Eng., ,

Universiti Teknologi Malaysia (UTM), 81310 Skudai, Johor, Malaysia
ezy(:,jebc+2) = ezy(:,jebc+2) + 61.4*sin(omega*n*dt)
(phone: 607-55-35785; fax: 607-55-66272; e-mail: akram@ fke.utm.my). is used to generate the incident wave at j = jebc + 2 and n
Y. Jasmy is with the Faculty of Health Science & Biomedical Eng., = 1. Where jebc is thickness of the PML region, n is time step
Universiti Teknologi Malaysia (UTM), 81310 Skudai, Johor, Malaysia.

International Scholarly and Scientific Research & Innovation 6(3) 2012 360 scholar.waset.org/1307-6892/5827
 World Academy of Science, Engineering and Technology
International Journal of Computer and Information Engineering
Vol:6, No:3, 2012

in FDTD code and dt is 5.0 psec. Fig. 2 shows a comparison of

FDTD calculated results with analytic solution along the center
axis o the cylinder. The simulation was made with a code
attached in Appendix.

Fig. 1 Schematic representation of the FDTD computational domain

Fig. 3 Test and benchmark computational domains
International Science Index, Computer and Information Engineering Vol:6, No:3, 2012 waset.org/Publication/5827

50

analytical solution
Then, the local error of the computed field in ΩT, due to
FDTD + PML ( 16 ) reflections from the PML, is obtained by subtracting the field
45
FDTD + PML ( 8 ) at any point within ΩT at a given time step from the field at the
corresponding point in ΩB. Here, Ez is used to define the error
40
n n
= Ez ,T − Ez , B i , j
n
|E|

elocal i, j i, j
(8)
35

where Ez,T and Ez,B are, respectively, the FDTD computed E-

30 fields within the test and benchmark domains. Further, the
global error was defined as
25
5 10 15 20 25 30 35 40 45 2
eglobal = ∑∑ Ez ,T
n n n
Grid Position, I
i, j
− Ez , B i , j (9)
Fig. 2 Waveform comparison of the FDTD results (with different i j
PML thicknesses) and the analytical solution for a plane wave source
In Fig. 4 and 5, different numbers of PML are used outside
III. PROCEDURE FOR PAPER SUBMISSION the FDTD simulation region to test the efficiency of the PML
In this section, various numerical simulations are performed ABC for the lossy media. The accuracy of the PML ABC is
to illustrate the effectiveness of PML code [8]. Fig. 3 shows compared with that of standard second-order Mur ABC [9].
the two FDTD computational domains used in the The numbers of PML used are 16, 8 and 4, respectively. As is
experiments: a test domain, ΩT, and the much larger evident from the result, 16 PML are good enough to absorb
benchmark domain ΩB. By calculating the difference between most of the incident waves at the boundaries.
the FDTD solutions in the two domains at each grid-point at
each time-step, a measure of the spurious reflection caused by IV. CONCLUSION
the test ABC is obtained. Initially, a 2-D test domain 100×50 The FDTD method is widely used because it is simple to
cells long was examined. A square cell was used, with δ implement numerically. It provides a flexible means for
= 0.015 m, and the time step was chosen based on the stability directly solving Maxwell’s time-dependent curl equations by
criterion δ t = δ / 2c , where c is the speed of light in vacuum. using finite differences to discretize them. This paper
In all the cases examined here, the excitation was a pulse describes the design of two-dimensional (2D) FDTD
exhibiting a very smooth transition to zero, as used in [8], and simulation software for transverse magnetic (TM) polarized
defined as follows, incident wave using Berenger's split-field perfectly matched
 1 10 − 15cos (π n / 20 ) +  layer (PML) formulation. The software is developed using
   n ≤ 40 (7)
Ez ,T 50,25 =  32 6 cos ( 2π n / 20 ) − cos ( 3π n / 20 )  Matlab programming language. The numerical results agree
n

 very well with the analytical results for homogeneous media.

0 n > 40

International Scholarly and Scientific Research & Innovation 6(3) 2012 361 scholar.waset.org/1307-6892/5827
 World Academy of Science, Engineering and Technology
International Journal of Computer and Information Engineering
Vol:6, No:3, 2012

0
10 iefbc = ie + 2*iebc;
jefbc = je + 2*jebc;
ibfbc = iefbc + 1;
-10
10 jbfbc = jefbc + 1;
%**********************************************
% Material parameters
-20
%**********************************************
Global Error

10
media = 2;
Mur ABC
eps = [1.0 4.0];
PML ABC (16 cell)
10
-30
PML ABC (8 cell)
sig = [0.0 0.12];
PML ABC (4 cell) mur = [1.0 1.0];
sim = [0.0 0.0];
10
-40
%**********************************************
% Field arrays
%**********************************************
10
-50 ez = zeros(ibfbc,jbfbc); %fields in main grid
0 100 200 300 400 500 ezx = zeros(ibfbc,jbfbc);
Time Step
ezy = zeros(ibfbc,jbfbc);
Fig. 4 Global error within a 100×50 cell two-dimensional TM test hx = zeros(iefbc,jefbc);
grid hy = zeros(iefbc,jefbc);
%**********************************************
International Science Index, Computer and Information Engineering Vol:6, No:3, 2012 waset.org/Publication/5827

-5
x 10
2 % Updating coefficients
%**********************************************
1.5 for i = 1:media
Normalized Reflection Along J = 1

eaf = dt*sig(i)/(2.0*epsz*eps(i));
1
ca(i) = (1.0-eaf)/(1.0+eaf);
cb(i) = dt/epsz/eps(i)/dx/(1.0+eaf);
haf = dt*sim(i)/(2.0*muz*mur(i));
0.5
da(i) = (1.0-haf)/(1.0+haf);
db(i) = dt/muz/mur(i)/dx/(1.0+haf);
0 end
%**********************************************
-0.5 % Geometry specification (main grid)
PML ABC ( 16 cell )
PML ABC ( 8 cell ) %**********************************************
-1 % Initialize entire main grid to free space
caezx(1:ibfbc,1:jbfbc) = ca(1);
-1.5 cbezx(1:ibfbc,1:jbfbc) = cb(1);
0 20 40 60 80 100 caezy(1:ibfbc,1:jbfbc) = ca(1);
Grid Position, I
cbezy(1:ibfbc,1:jbfbc) = cb(1);
Fig. 5 Comparative local error measures for Mur and PML ABCs dahxy(1:iefbc,1:jefbc) = da(1);
along the test grid outer boundary at time-step n = 100. dbhxy(1:iefbc,1:jefbc) = db(1);
dahyx(1:iefbc,1:jefbc) = da(1);
dbhyx(1:iefbc,1:jefbc) = db(1);
APPENDIX
%********************************************** % Add dielectric cylinder
% Fundamental constants diam = 40; % diameter of cylinder: 6 cm
%********************************************** rad = diam/2.0; % radius of cylinder: 3 cm
cc = 2.99792458e8; %speed of light in free space ic = iefbc/2; % i-coordinate of cylinder's center
muz = 4.0*pi*1.0e-7; %permeability of free space jc = jefbc/2; % j-coordinate of cylinder's center
epsz = 1.0/(cc*cc*muz); %permittivity of free space
freq = 2.5e+9; %frequency of source excitation for i = iebc+1:ie+iebc-1
omega = 2.0*pi*freq; for j = jebc+1:je+jebc-1
%********************************************** dist2 = (i-ic)^2 + (j-jc)^2;
% Grid parameters if dist2 <= rad^2
%********************************************** caezy(i,j) = ca(2);
ie = 50; %number of grid cells in x-direction cbezy(i,j) = cb(2);
je = 50; %number of grid cells in y-direction caezx(i,j) = ca(2);
cbezx(i,j) = cb(2);
ib = ie + 1; end
jb = je + 1; end
end
dx = 3.0e-3; %space increment of square lattice %**********************************************
dt = dx/(2.0*cc); %time step % Fill the PML regions
%**********************************************
nmax = 500; %total number of time steps delbc = iebc*dx;
sigmam=-log(rmax/100.0)*epsz*cc*(orderbc+1)/(2*delbc);
iebc = 8; %thickness of left and right PML region bcfactor=eps(1)*sigmam/(dx*(delbc^orderbc) * …
jebc = 8; %thickness of front and back PML region (orderbc+1));
rmax = 0.00001;
orderbc = 2; % FRONT region
ibbc = iebc+1; caezy(1:ibfbc,1) = 1.0;
jbbc = jebc+1; cbezy(1:ibfbc,1) = 0.0;

International Scholarly and Scientific Research & Innovation 6(3) 2012 362 scholar.waset.org/1307-6892/5827
 World Academy of Science, Engineering and Technology
International Journal of Computer and Information Engineering
Vol:6, No:3, 2012

for j = 2:jebc
y1 = (jebc-j+1.5)*dx; for i = 1:iebc
y2 = (jebc-j+0.5)*dx; x1 = (iebc-i+1)*dx;
sigmay = bcfactor*(y1^(orderbc+1)-y2^(orderbc+1)); x2 = (iebc-i)*dx;
ca1 = exp(-sigmay*dt/(epsz*eps(1))); sigmax = bcfactor*(x1^(orderbc+1)-x2^(orderbc+1));
cb1 = (1.0-ca1)/(sigmay*dx); sigmaxs = sigmax*(muz/(epsz*eps(1)));
caezy(1:ibfbc,j) = ca1; da1 = exp(-sigmaxs*dt/muz);
cbezy(1:ibfbc,j) = cb1; db1 = (1-da1)/(sigmaxs*dx);
end dahyx(i,1:jefbc) = da1;
sigmay = bcfactor*(0.5*dx)^(orderbc+1); dbhyx(i,1:jefbc) = db1;
ca1 = exp(-sigmay*dt/(epsz*eps(1))); end
cb1 = (1-ca1)/(sigmay*dx);
caezy(1:ibfbc,jbbc) = ca1; % RIGHT region
cbezy(1:ibfbc,jbbc) = cb1; caezx(ibfbc,1:jbfbc) = 1.0;
cbezx(ibfbc,1:jbfbc) = 0.0;
for j=1:jebc for i = 1:iebc-1
y1 = (jebc-j+1)*dx; x1 = (i+0.5)*dx;
y2 = (jebc-j)*dx; x2 = (i-0.5)*dx;
sigmay = bcfactor*(y1^(orderbc+1)-y2^(orderbc+1)); sigmax = bcfactor*(x1^(orderbc+1)-x2^(orderbc+1));
sigmays = sigmay*(muz/(epsz*eps(1))); ca1 = exp(-sigmax*dt/(epsz*eps(1)));
da1 = exp(-sigmays*dt/muz); cb1 = (1-ca1)/(sigmax*dx);
International Science Index, Computer and Information Engineering Vol:6, No:3, 2012 waset.org/Publication/5827

db1 = (1-da1)/(sigmays*dx); caezx(ibfbc-iebc+i,1:jbfbc) = ca1;

dahxy(1:iefbc,j) = da1; cbezx(ibfbc-iebc+i,1:jbfbc) = cb1;
dbhxy(1:iefbc,j) = db1; end
end sigmax = bcfactor*(0.5*dx)^(orderbc+1);
ca1 = exp(-sigmax*dt/(epsz*eps(1)));
% BACK region cb1 = (1-ca1)/(sigmax*dx);
caezy(1:ibfbc,jbfbc) = 1.0; caezx(ibfbc-iebc+1,1:jbfbc) = ca1;
cbezy(1:ibfbc,jbfbc) = 0.0; cbezx(ibfbc-iebc+1,1:jbfbc) = cb1;
for j = 1:jebc-1
y1 = (j+0.5)*dx; for i = 1:iebc
y2 = (j-0.5)*dx; x1 = i*dx;
sigmay = bcfactor*(y1^(orderbc+1)-y2^(orderbc+1)); x2 = (i-1)*dx;
ca1 = exp(-sigmay*dt/(epsz*eps(1))); sigmax = bcfactor*(x1^(orderbc+1)-x2^(orderbc+1));
cb1 = (1-ca1)/(sigmay*dx); sigmaxs = sigmax*(muz/(epsz*eps(1)));
caezy(1:ibfbc,jbfbc-jebc+j) = ca1; da1 = exp(-sigmaxs*dt/muz);
cbezy(1:ibfbc,jbfbc-jebc+j) = cb1; db1 = (1-da1)/(sigmaxs*dx);
end dahyx(iefbc-iebc+i,1:jefbc) = da1;
sigmay = bcfactor*(0.5*dx)^(orderbc+1); dbhyx(iefbc-iebc+i,1:jefbc) = db1;
ca1 = exp(-sigmay*dt/(epsz*eps(1))); end
cb1 = (1-ca1)/(sigmay*dx); %**********************************************
caezy(1:ibfbc,jefbc-jebc+1) = ca1; % BEGIN TIME-STEPPING LOOP
cbezy(1:ibfbc,jefbc-jebc+1) = cb1; %**********************************************
for n = 1:nmax
for j = 1:jebc %**********************************************
y1 = j*dx; % Update electric field EZ in main grid
y2 = (j-1)*dx; %**********************************************
sigmay = bcfactor*(y1^(orderbc+1)-y2^(orderbc+1)); ezy(:,jebc+2) = ezy(:,jebc+2) + 61.4*sin(omega*n*dt);
sigmays = sigmay*(muz/(epsz*eps(1))); ezx(2:iefbc,2:jefbc)= caezx(2:iefbc,2:jefbc).*ezx(2:iefbc,2:jefbc) + ...
da1 = exp(-sigmays*dt/muz); cbezx(2:iefbc,2:jefbc).*(hyx(2:iefbc,2:jefbc)- hyx(1:iefbc-1,2:jefbc));
db1 = (1-da1)/(sigmays*dx);
dahxy(1:iefbc,jefbc-jebc+j) = da1; ezy(2:iefbc,2:jefbc)= caezy(2:iefbc,2:jefbc).*ezy(2:iefbc,2:jefbc) + ...
dbhxy(1:iefbc,jefbc-jebc+j) = db1; cbezy(2:iefbc,2:jefbc).*(hxy(2:iefbc,1:jefbc-1)- hxy(2:iefbc,2:jefbc));
end
ez = ezx + ezy;
% LEFT region %**********************************************
caezx(1,1:jbfbc) = 1.0; % Update magnetic fields (Hx and Hy) in main grid
cbezx(1,1:jbfbc) = 0.0; %**********************************************
for i = 2:iebc hxy(:,:) = dahxy(:,:).*hxy(:,:) + dbhxy(:,:).*...
x1 = (iebc-i+1.5)*dx; (ezx(1:iefbc,1:jefbc) + ezy(1:iefbc,1:jefbc) - ...
x2 = (iebc-i+0.5)*dx; ezx(1:iefbc,2:jbfbc) - ezy(1:iefbc,2:jbfbc));
sigmax = bcfactor*(x1^(orderbc+1)-x2^(orderbc+1));
ca1 = exp(-sigmax*dt/(epsz*eps(1))); hyx(:,:) = dahyx(:,:).*hyx(:,:) + dbhyx(:,:).*...
cb1 = (1-ca1)/(sigmax*dx); (ezx(2:ibfbc,1:jefbc) + ezy(2:ibfbc,1:jefbc) - ...
caezx(i,1:jbfbc) = ca1; ezx(1:iefbc,1:jefbc) - ezy(1:iefbc,1:jefbc));
cbezx(i,1:jbfbc) = cb1;
end %************************************************
sigmax = bcfactor*(0.5*dx)^(orderbc+1); % END TIME-STEPPING LOOP
ca1 = exp(-sigmax*dt/(epsz*eps(1))); %************************************************
cb1 = (1-ca1)/(sigmax*dx); end
caezx(ibbc,1:jbfbc) = ca1;
cbezx(ibbc,1:jbfbc) = cb1;

International Scholarly and Scientific Research & Innovation 6(3) 2012 363 scholar.waset.org/1307-6892/5827
 World Academy of Science, Engineering and Technology
International Journal of Computer and Information Engineering
Vol:6, No:3, 2012

REFERENCES
[1] K.S. Yee, “Numerical solution of initial boundary value problems
involving Maxwell’s equations in isotropic media,” IEEE Trans.
Antennas & Propag., vol. 14, 1966, pp: 302-307.
[2] A. Taflove, and M.E. Brodwin, "Numerical solution of steady-state
electromagnetic scattering problems using the time-dependent
Maxwell's equations," IEEE Trans. on Microwave Theo. & Tech., vol.
23, no. 8, 1975, pp.:623–630.
[3] A. Taflove, ” Application of the finite-difference time-domain method
to sinusoidal steady state electromagnetic penetration problems,” IEEE
Trans. Electromag. Compatibil., vol. 22, 1980, pp.: 191-202.
[4] A. Taflove, S. Hagness, “ Computational Electrodynamics: The Finite-
Difference Time-Domain Method,” 3rd Ed.. Artech House Publishers.
2005.
[5] K. Umashankar, and A. Taflove, "A novel method to analyze
electromagnetic scattering of complex objects," IEEE Trans.
Electromag. Compatibil., vol. 24, 1982, pp.:397–405.
[6] J.P. Berenger, “A perfectly matched layer for the absorption of
electromagnetic waves,” J. Comput. Phys., vol. 114, 1994, pp.: 185-200.
[7] R. Harrington, “Time Harmonic Electromagnetic Fields,” New York:
International Science Index, Computer and Information Engineering Vol:6, No:3, 2012 waset.org/Publication/5827

McGraw-Hill. 1961.
[8] T. Moore, J. Blaschak, A. Taflove, and G. Kriegsmann, “Theory and
application of radiation boundary operators,” IEEE Trans. Antennas &
Propagat., vol. 36, 1988, pp. 1797-1812.
[9] G. Mur, “Absorbing boundary conditions for the finite-difference
approximation of the time-domain electromagnetic-field equations,”
IEEE Trans. Electromag. Compatibil., vol. 23, no. 4, 1981, pp. 377-382.

G. Akram (M’06) Gasmelseed Akram received his B.Sc. degree in

Electrical Engineering and Informatics – major in Computer Engineering –
and M.Sc degree in Electrical Engineering and Informatics – major in
Biomedical Engineering – from Budapest University of Technology and
Economics (BME), Budapest, Hungary, in 1993 and 1999, respectively.He
received the PhD degree in Electrical Engineering – major in Biomedical
Engineering – from Universiti Teknologi Malasysia (UTM), Malaysia, in
2009. He is the founder and chair of IEEE EMBS Scoiety. His biography has
listed in the 25th Silver Anniversary Edition of Who’s Who in the World
2009 and 2011.
Y. Jasmy is a Professor at the Faculty of Health Science & Biomedical
Engineering, Universiti Teknologi Malaysia (UTM). He obtained his MSc
and PhD in Electronics from the University of Kent at Canterbury, UK. He
has been involved in developing the Medical Electronics area in the faculty,
heading the computer aided rehabilitation engineering group. His research is
in speech therapy, electromagnetic effects in hospital environment, and
computers as aids for the handicapped.

International Scholarly and Scientific Research & Innovation 6(3) 2012 364 scholar.waset.org/1307-6892/5827