The FDTD method can be classed as a special case of the finite difference method , a standard method for solving PDEs. It is a time-stepping method, wherein starting at time t=0, the electromagnetic properties are updated at each time level to give the evolving solution.

The FDTD method can be used in 1, 2 and 3 dimensions, but it is the more practical 3D case that is the most important and which the Yee algorithm addresses. It is the 3 dimensional FDTD solution of electromagnetic problems that is considered here.

1  Modelling

For a linear, isotropic nondispersive materials (i.e. materials having field-independent, direction-independent, and frequency-independent electric and magnetic properties), then Maxwell's curl equations take the following form:

¶H ¶t = - 1 m Ñ×E - r m H ,

(1)

¶E ¶t = 1 e Ñ×H - s e E ,

(2)

where H=[H_x,H_y,H_z] is the magnetic field intensity, E=[E_x,E_y,E_z] is the electric field intensity, m is the permeability, e is the permittivity, r is the magnetic resistivity and s is the conductivity. (Note that the magnetic resistivity (r) is generally taken to be zero. Its inclusion above is simply to demonstrate the similar form of the two equations.) In general m, r, e and s can vary throughout the medium. For the background to these equations, click here .

The individual components of the equations (1) and (2) are as follows

¶Hx ¶t = - 1 m ( ¶Ey ¶z - ¶Ez ¶y - rHx ) ,

(3)

¶Hy ¶t = - 1 m ( ¶Ez ¶x - ¶Ex ¶z - rHy) ,

(4)

¶Hz ¶t = - 1 m ( ¶Ex ¶y - ¶Ey ¶x - rHz ) ,

(5)

from equation (1) and

¶Ex ¶t = - 1 e ( ¶Hz ¶y - ¶Hy ¶z - sEx ) ,

(6)

¶Ey ¶t = - 1 e ( ¶Hx ¶z - ¶Hz ¶x - sEy) ,

(7)

¶Ez ¶t = - 1 e ( ¶Hy ¶x - ¶Hx ¶y - sEz ) .

(8)

The system of six coupled PDEs forms the basis of the FD-TD numerical algorithm.

2  The Yee Algorithm

In 1966, Kane Yee [2] introduced a set of finite-difference equations for the time-dependent Maxwell's curl equations system for a lossless material (r = 0 and s = 0). In order to derive the method it is assumed that the three-dimensional region of interest is divided into a mesh of cubic cells. Equation (2) is used to update E at each time level (n), equation (1) is used to update H at each mid time level (n+[1/2]). The mesh used for updating E and H is different; the edge of the cubes in one mesh lie at the centres of the cubes in the other mesh; the cubes interlace each other.

The components of the electromagnetic field E_x, E_y, E_z, H_x, H_y, H_z are computed at different points in the domain. The above figure shows the mesh for computing H from E with a typical origin. The coordinates show the positioning of the electromagnetic components with respect to that origin. The finite difference equivalent of equation (1) for computing H at the [1/2] time level from the original time level is as follows:

Hx|[1/2],0,0[1/2] = A[1/2],0,0Hx|[1/2],0,0-[1/2]+ C[1/2],0,0 é ê ë - Ez|[1/2],[1/2],00 - Ez|[1/2],-[1/2],00 Dy + Ey|[1/2],0,[1/2][1/2] - Ey|[1/2],0,-[1/2]0 Dz ù ú û

Hy|0,[1/2],0[1/2] = A0,[1/2],0 Hy|0,[1/2],0-[1/2]+ C0,[1/2],0 é ê ë - Ex|0,[1/2],[1/2]0 - Ex|0,[1/2],-[1/2]0 Dz + Ez|[1/2],[1/2],0[1/2] - Ez|-[1/2],[1/2],00 Dx ù ú û

Hz|0,0,[1/2][1/2] = A0,0,[1/2] Hz|0,0,[1/2]-[1/2]+C0,0,[1/2] é ê ë - Ey|[1/2],0,[1/2]0 - Ey|-[1/2],0,[1/2]0 Dx + Ex|0,[1/2],[1/2][1/2] - Ey|0,-[1/2],[1/2]0 Dy ù ú û

where

A*,*,* = 1.0 ,C*,*,* = Dt m*,*,*  ,

and parameter r is presumed zero.

For computing E from H the following figure shows a cube from the mesh in the region of the same origin. The following finite difference equivalents of equation (2) may be derived

Ex|1,[1/2],[1/2]1 = B1,[1/2],[1/2] Ex|1,[1/2],[1/2]0+ D1,[1/2],[1/2] é ê ë Hz|1,1,[1/2][1/2] - Hz|1,0,[1/2][1/2] Dy - Hy|1,[1/2],1[1/2] - Hy|1,[1/2],0[1/2] Dz ù ú û

Ey|[1/2],1,[1/2]1 = B[1/2],1,[1/2] Ey|[1/2],1,[1/2]0+ D[1/2],1,[1/2] é ê ë Hx|[1/2],1,1[1/2] - Hx|[1/2],1,0[1/2] Dz - Hz|1,1,[1/2][1/2] - Hz|0,1,[1/2][1/2] Dx ù ú û

Ez|[1/2],[1/2],11 = B[1/2],[1/2],1 Ez|[1/2],[1/2],10+ D[1/2],[1/2],1 é ê ë Hy|1,[1/2],1[1/2] - Hy|0,[1/2],1[1/2] Dx - Hx|[1/2],1,1[1/2] - Hx|[1/2],0,1[1/2] Dy ù ú û

where

B*,*,* = 1- s*,*,* Dt 2e*,*,* 1+ s*,*,* Dt 2e*,*,* ,         D*,*,* = Dt e*,*,* 1+ s*,*,* Dt 2e*,*,*

Most simply, the complete electromagnetic domain of interest is taken to be a cuboid with N_x, N_y, and N_z cuboids in the x,y and z directions. From the information given by the components in this cell, H_x,H_y,H_z at the n-[1/2]^th time level and E_x,E_y,E_z at the n^th time level, H_x,H_y,H_z at the n+[1/2]^th time level can be computed using the following formulae:

For i=1..N_x, j=1..N_y, k=1..N_z

Hx|i+[1/2],j,kn+[1/2] = Ai+[1/2],j,kHx|i+[1/2],j,kn-[1/2]+ Ci+[1/2],j,k é ê ë - Ez|i+[1/2],j+[1/2],kn - Ez|i+[1/2],j-[1/2],kn Dy + Ey|i+[1/2],j,k+[1/2]n+[1/2] - Ey|i+[1/2],j,k-[1/2]n Dz ù ú û

Hy|i,j+[1/2],kn+[1/2] = Ai,j+[1/2],k Hy|i,j+[1/2],kn-[1/2]+ Ci,j+[1/2],k é ê ë - Ex|i,j+[1/2],k+[1/2]n - Ex|i,j+[1/2],k-[1/2]n Dz + Ez|i+[1/2],j+[1/2],kn+[1/2] - Ez|i-[1/2],j+[1/2],kn Dx ù ú û

Hz|i,j,k+[1/2]n+[1/2] = Ai,j,k+[1/2] Hz|i,j,k+[1/2]n-[1/2]+Ci,j,k+[1/2] é ê ë - Ey|i+[1/2],j,k+[1/2]n - Ey|i-[1/2],j,k+[1/2]n Dx + Ex|i,j+[1/2],k+[1/2]n+[1/2] - Ey|i,j-[1/2],k+[1/2]n Dy ù ú û

where

A*,*,* = 1,         C*,*,* = Dt m*,*,*

Ex|i+1,j+[1/2],k+[1/2]n+1 = Bi+1,j+[1/2],k+[1/2] Ex|i+1,j+[1/2],k+[1/2]n +

Di+1,j+[1/2],k+[1/2] é ê ë Hz|i+1,j+1,k+[1/2]n+[1/2] - Hz|i+1,j,k+[1/2]n+[1/2] Dy - Hy|i+1,j+[1/2],k+1n+[1/2] - Hy|i+1,j+[1/2],kn+[1/2] Dz ù ú û

Ey|i+[1/2],j+1,k+[1/2]n+1 = Bi+[1/2],j+1,k+[1/2] Ey|i+[1/2],j+1,k+[1/2]n+

Di+[1/2],j+1,k+[1/2] é ê ë Hx|i+[1/2],j+1,k+1n+[1/2] - Hx|i+[1/2],j+1,kn+[1/2] Dz - Hz|i+1,j+1,k+[1/2]n+[1/2] - Hz|i,j+1,k+[1/2]n+[1/2] Dx ù ú û

Ez|i+[1/2],j+[1/2],k+1n+1 = Bi+[1/2],j+[1/2],k+1 Ez|i+[1/2],j+[1/2],k+1n+

Di+[1/2],j+[1/2],k+1 é ê ë Hy|i+1,j+[1/2],k+1n+[1/2] - Hy|i,j+[1/2],k+1n+[1/2] Dx - Hx|i+[1/2],j+1,k+1n+[1/2] - Hx|i+[1/2],j,k+1n+[1/2] Dy ù ú û

where

B*,*,* = 1- s*,*,* Dt 2e*,*,* 1+ s*,*,* Dt 2e*,*,* ,         D*,*,* = Dt e*,*,* 1+ s*,*,* Dt 2e*,*,*

Practicalities