IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES,VOL. 49, NO. 2, FEBRUARY 2001 377 Short Papers_______________________________________________________________________________ Analyzing the Stability of the FDTD Technique by Combining the von Neumann Method with the Routh–Hurwitz Criterion José A. Pereda, Luis A. Vielva, Ángel Vegas, and Andrés Prieto Abstract—This paper addresses the problem of stability analysis of fi- nite-difference time-domain (FDTD) approximations for Maxwell’s equa- tions. The combination of the von Neumann method with the Routh–Hur- witz criterion is proposed as an algebraic procedure for obtaining analytical closed-form stability expressions. This technique is applied to the problem of determining the stability conditions of an extension of the FDTD method to incorporate dispersive media previously reported in the literature. Both Debye and Lorentz dispersive media are considered. It is shown that, for the former case, the stability limit of the conventional FDTD method is preserved. However, for the latter case, a more restrictive stability limit is obtained. To overcome this drawback, a new scheme is presented, which al- lows the stability limit of the conventional FDTD method to be maintained. Index Terms—Dispersive media, FDTD methods, Routh’s methods, sta- bility analysis, von Neumann method. I. INTRODUCTION The finite-difference time-domain (FDTD) method is a powerful numerical technique for the solution of electromagnetic problems. Based on the finite-difference approximation of the time-dependent Maxwell’s curl equations, this method leads to a conditionally stable scheme. The stability of a difference scheme refers to the unstable growth or stable decay of errors in the arithmetic operations required to solve the finite-difference equations. Due to the conditional stability of the FDTD method, it is essential to choose a number of parameters—time step, size of the spatial mesh, etc.—in such a way that the difference scheme remains stable. Thus, the derivation of closed-form analytical stability conditions is of great theoretical and practical interest. The establishment of the analytical stability condition for the con- ventional FDTD method (Yee’s scheme for isotropic, nondispersive lossless media) was an early development [1]. More recently, stability conditions were determined for extensions of the conventional FDTD method to incorporate dispersive media [2], [3], anisotropic media, and lossy media [5], [6]. The von Neumann method is probably the most popular procedure for determining the stability of finite-difference approximations for space–time problems [7]. For a given difference scheme, this technique provides an associated polynomial. The condition for stability is that all the roots of this polynomial have modulus less than or equal to unity. When applying this method, the root locus is obtained as a function of the parameters of interest. This technique usually requires exhaustive numerical root searching, which makes the derivation of closed-form stability conditions difficult. The numerical searching may also intro- Manuscript received June 18, 1999. This work was supported by the Comisión Interministerial de Ciencia y Tecnología, Spain, under Research Project TIC2000-1612-C03-01. The authors are with the Departamento de Ingeniería de Comunicaciones, Universidad de Cantabria, 39005 Santander, Cantabria, Spain (e-mail: pereda@dicom.unican.es). Publisher Item Identifier S 0018-9480(01)01080-8. duce inaccuracies in the location of the roots and, therefore, in the ranges of stability for the parameters of interest. The Routh–Hurwitz (R–H) criterion is a procedure widely used to study the stability of continuous- and discrete-time linear systems [8]. This method is able to establish the location of the zeros of a real-coef- ficient polynomial with respect to the imaginary axis, without actually solving for the zeros. This paper proposes combining the von Neumann method with the R–H criterion as an algebraic procedure for deriving closed-form sta- bility expressions for the FDTD method. To illustrate this technique, it is applied to the stability analysis of an extension of the FDTD method, previously published, which is able to incorporate dispersive media [9]. It is shown that the scheme introduced in [9] to treat Debye media pre- serves the stability limit of the conventional FDTD method. However, a more restrictive stability limit is obtained for the scheme given in [9] to treat Lorentz media. To overcome this limitation, a new scheme is presented, which recovers the stability limit of the conventional FDTD method. II. THEORY A. von Neumann Method Each difference scheme has a theoretical exact solution. However, when explicit calculations are carried out in a computer, errors are committed due to the finite precision of the arithmetic operations. The study of the stability of a finite-difference scheme consists of finding the conditions under which the error—difference between the theoret- ical and numerical solutions of the finite-difference equation—remains bounded as the number of time iterations tends to infinity. The von Neumann method basically consists of considering a Fourier series expansion of the error at the mesh nodes at a given time instant . Due to linearity, only a single term of this expansion needs to be considered, i.e., (1) where is a complex amplitude, indexes , , denote the position of the nodes in the mesh, are the sizes of the discretiza- tion cell, and are the wavenumbers of the discrete modes in the -direction. In (1), is a complex variable, often called the amplification factor, which gives the growth of the error in a time iteration, i.e., . To ensure that a finite-difference scheme will be stable, the error must not grow as the time increases, thus, the condition must be satisfied. To derive the stability condition for a particular scheme, solutions of the form (1) are substituted into the difference equations. This leads to a characteristic polynomial in (2) The condition for stability can now be written as where are the roots of . That is, the roots of must lie inside or on the unit circle in the -plane. To determine the range of the parameters of interest for which the scheme is stable, it is usually 0018–9480/01$10.00 © 2001 IEEE