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
CourantFriedrichsLewy (CFL) criterion which is explained here.
Courant-Friedrich-Lewy (CFL)
criterion
CourantFriedrichsLewy (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}\]
