582 IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES, VOL. 48, NO. 4, APRIL 2000 Reduction of Numerical Dispersion in FDTD Method Through Artificial Anisotropy Jaakko S. Juntunen and Theodoros D. Tsiboukis, Member, IEEE Abstract—In this paper, a simple and computationally low-cost modification of the standard finite-difference time-domain (FDTD) algorithm is presented to reduce numerical dispersion in the algo- rithm. Both two- and three-dimensional cases are considered. It is shown that the maximum error in phase velocity can be reduced by a factor of 2–7, depending on the shape of the FDTD cell. Although the reduction procedure is optimal for only single frequency, nu- merical examples show that the proposed method can also improve the accuracy significantly in wide-band inhomogeneous problems. Index Terms—FDTD method, numerical dispersion. I. INTRODUCTION N UMERICAL dispersion is an undesired nonphysical effect inherently present in the finite-difference time-do- main (FDTD) algorithm. In short, numerical dispersion means dependence of wave propagation velocity on frequency. Herein, we also include in the term the dependence of velocity on propagation direction. The latter is sometimes called numerical anisotropy. In qualitative terms, dispersion causes distortion of waveforms. Frequency dependence causes high-frequency content of a wave to lag, while direction dependence causes spherical waveforms to become slightly cubical. There are several problems associated with numerical disper- sion. First, it causes cumulative phase error. If a device is based on phase cancellation, even an apparently small error in wave propagation velocity may cumulate phase error to unacceptable amounts. Equivalently, phase error appears as mislocation of resonances in the frequency domain. Sometimes, it is possible to pre-estimate the effect of the dispersion error, and use the es- timation to choose proper spatial resolution for the problem [1]. In some cases, the numerical dispersion can be eliminated in post-processing [2]. This elimination is rarely possible since it is based on the assumption that there are waves propagating in only one direction. Another possible trouble with the numerical dispersion is nonphysical refraction [3]. If the cell shape varies over the grid, a wave experiences different numerical dispersion in different parts of the grid. This corresponds to inhomogeneous medium, and refraction takes place. In some problems, it is indeed necessary to vary the cell shape quite dramatically, e.g., Manuscript received March 17, 1999. This work was supported by the Jenny and Antti Wihuri Foundation, by the Finnish Graduate School of Electronics, Telecommunication, and Automation, and under EU Grant ERBFMBICT983462. J. S. Juntunen is with the Radio Laboratory, Helsinki University of Tech- nology, FIN-02 015 HUT Espoo, Finland (e-mail: jju@radio.hut.fi). T. D. Tsiboukis is with the Electrical and Computer Engineering Department, Aristotle University of Thessaloniki, GR-54 006 Thessaloniki, Greece. Publisher Item Identifier S 0018-9480(00)02781-2. in [4], the width-to-length ratio of a two–dimensional (2-D) cell varies a factor of 12.5. A different FDTD algorithm is proposed in [5], which is equivalent to the so-called symmetrical-condensed-node trans- mission-line matrix method (SCN–TLM). The dispersion errors in SCN–TLM are less than in the standard FDTD technique [6], [7]. A distinct disadvantage of the SCN–TLM formulation though, is the need of extensive memory requirements as opposed to the standard FDTD implementation. Another possi- bility is to use fourth-order spatial differencing in the algorithm [8]. However, associated problems are encountered as more smoothness is assumed from the field quantities, especially on the material boundaries. The present reduction method is based on carefully speeding up the wave propagation by introducing anisotropy parameters into the algorithm. A detailed Fourier-mode analysis is given for the determination of the optimal anisotropy parameters. Several simulation examples confirm the theory. Wide-band problems are also discussed. II. NUMERICAL DISPERSION RELATION IN 2-D Let us consider electrically anisotropic medium in 2-D and the TE mode. Let the relative permittivity tensor be diagonal, i.e., . The stability condition of the FDTD algorithm for this problem is (1) Here, is the speed of light in free space. The derivation of the dispersion relation is canonical. A wave is expanded into plane waves of the form (2) and the FDTD update equations are applied to these waves. The resulting numerical dispersion relation for the TE mode is (3) The TM mode does not “see” the electric anisotropy. The dual situation can be obtained by replacing the relative permittivities by relative permeabilities. In (3), we write the numerical wave vector as , where is the numerical wavenumber and 0018–9480/00$10.00 © 2000 IEEE