Knowledge base

Numerical Stability

1 Numerical stability

In finite-difference time-domain (FDTD), the temporal resolution, \(\Delta t\), affects the stability of numerical solutions for Maxwell\('\)s equations. To ensure that there are bounded (or stable) solutions for discretized Maxwell\('\)s equations, the time step \(\Delta t\) should be selected fine enough relative to the spatial resolution. The dependency of the time step size on the spatial resolution in FDTD is expressed through the Courant–Friedrichs–Lewy (CFL) criterion which is explained here.

1.2 Courant-Friedrich-Lewy (CFL) criterion

Courant–Friedrichs–Lewy (CFL) criterion determines the maximum size of the time step in terms of the spatial discretization i.e. \(\Delta x\), \(\Delta y\), and \(\Delta z\) to avoid instability (unboundedness) of numerical solutions. As mentioned in article on numerical dispersion, the 2D numerical dispersion equation for a traveling plane wave in free space can be found as [1] \[\begin{equation} \left[\left(\frac{1}{c\Delta t}\right)\sin\left(\frac{\omega\Delta t}{2}\right)\right]^{2}=\left[\frac{1}{\Delta x}\sin\left(\frac{\tilde{k_{x}}\Delta x}{2}\right)\right]^{2}+\left[\frac{1}{\Delta z}\sin\left(\frac{\tilde{k_{z}}\Delta z}{2}\right)\right]^{2} \end{equation}\] where \(\tilde{k}_{x}\) and \(\tilde{k}_{z}\) are the numerical wavevectors in x and z direction. The \(\omega\) and \(c\) represent the angular frequency and light velocity in free space, respectively. From the above equation, the angular frequency can be obtained as \[\begin{equation} \omega=\frac{2}{\Delta t}\sin^{-1}\left(c\Delta t\sqrt{\left[\frac{1}{\Delta x}\sin\left(\frac{\tilde{k_{x}}\Delta x}{2}\right)\right]^{2}+\left[\frac{1}{\Delta z}\sin\left(\frac{\tilde{k_{z}}\Delta z}{2}\right)\right]^{2}}\right) \end{equation}\] Since the plane wave is traveling in free space, it does not experience any variation in the amplitude, implying a real-value angular frequency (note that an angular frequency with a complex value is equivalent to dissipation in the medium). To have a real-value angular frequency from (2), the following condition should be met [1] \[\begin{equation} c\Delta t\sqrt{\left[\frac{1}{\Delta x}\sin\left(\frac{\tilde{k_{x}}\Delta x}{2}\right)\right]^{2}+\left[\frac{1}{\Delta z}\sin\left(\frac{\tilde{k_{z}}\Delta z}{2}\right)\right]^{2}}\leq1 \end{equation}\] this gives rise to the following inequality (known as CFL criterion) \[\begin{equation} \Delta t\leq\frac{1}{c}\sqrt{\frac{1}{\frac{1}{\Delta x^{2}}+\frac{1}{\Delta z^{2}}}} \end{equation}\] For the three dimension (3D) case, the CFL criterion can be derived as \[\begin{equation} \Delta t\leq\frac{1}{c}\sqrt{\frac{1}{\frac{1}{\Delta x^{2}}+\frac{1}{\Delta y^{2}}+\frac{1}{\Delta z^{2}}}} \end{equation}\] The CFL criterion emphasizes that within a simulation grid the distance traveled by a wave in each time step cannot exceed the cell size. For a uniform mesh (\(\Delta x = \Delta z = \Delta\)) this condition can be re-expressed in terms of the “Courant factor” S \[\begin{eqnarray} S=\frac{c\Delta t}{\Delta} \end{eqnarray}\] where \(S\leq\frac{1}{\sqrt{2}}\) or \(S\leq\frac{1}{\sqrt{3}}\) for 2D and 3D, respectively. The factors of \(\frac{1}{\sqrt{2}}\) and \(\frac{1}{\sqrt{3}}\) are the limits of stability and therefore in to keep the simulation stable it is necessary to adhere to just under the stability criterion, with values typically set to \(0.99 S\). This allows (5) to be rewritten as \[\begin{equation} \Delta t=\frac{0.99}{c}\sqrt{\frac{1}{\frac{1}{\Delta x^{2}}+\frac{1}{\Delta y^{2}}+\frac{1}{\Delta z^{2}}}} \end{equation}\]


  1. A. Taflove, S. C. Hagness, Computational Electrodynamics: The
    Finite-Difference Time-Domain Method 3rd edition. Artech House,


Learn More About Our 30-Day Evaluations