# Finite Difference Time Domain (FDTD) Electromagnetics

(lightly adapted from: Understanding the Finite-Difference Time-Domain Method, John B. Schneider, http://www.eecs.wsu.edu/~schneidj/ufdtd, 2010.)
 
The FDTD method employs finite differences as approximations to both
the spatial and temporal derivatives that appear in Maxwell's
equations (specifically Ampere's and Faraday's laws). 

## Maxwell's equations

Maxwell's equations in differential form are:

\begin{equation}
\nabla \cdot \mathbf{D} = \rho \quad \text{(Gauss's Law)}
\end{equation}

\begin{equation}
\nabla \cdot \mathbf{B} = 0 \quad \text{(Gauss's Law for magentism)}
\end{equation}


\begin{equation}
\nabla \times \mathbf{E} = -\frac{\partial\mathbf{B}}{\partial t} \quad \text{(Faraday's Law)}
\end{equation}

\begin{equation}
\nabla \times \mathbf{H} = \mathbf{J} + \frac{\partial\mathbf{D}}{\partial t} \quad   \text{(Ampere's Law)}
\end{equation}

where: $\mathbf{D}$ is the electric displacement field, $\rho$ is the electric charge density, $\mathbf{B}$ is the magnetic flux density, $\mathbf{E}$ is the electric field, $\mathbf{H}$ is the magnetic field, and $\mathbf{J}$ is the electric current density.

In a region with no sources and no currents, they become:

\begin{equation}
\nabla \cdot \mathbf{D} = 0
\end{equation}

\begin{equation}
\nabla \cdot \mathbf{B} = 0
\end{equation}


\begin{equation}
\nabla \times \mathbf{E} = -\frac{\partial\mathbf{B}}{\partial t}
\end{equation}

\begin{equation}
\nabla \times \mathbf{H} = \frac{\partial\mathbf{D}}{\partial t}
\end{equation}

For materials without polarization and magnetization:

\begin{equation}
\mathbf{D} = \epsilon \mathbf{E}
\end{equation}

\begin{equation}
\mathbf{B} = \mu \mathbf{H}
\end{equation}


And the equations become

\begin{equation}
\nabla \cdot \mathbf{E} = 0
\end{equation}

\begin{equation}
\nabla \cdot \mathbf{H} = 0
\end{equation}


\begin{equation}
\nabla \times \mathbf{E} = -\mu \frac{\partial\mathbf{H}}{\partial t}
\end{equation}

\begin{equation}
\nabla \times \mathbf{H} = \epsilon \frac{\partial\mathbf{E}}{\partial t}
\end{equation}

(Note that these equations are written in terms of $\mathbf{E}$ and $\mathbf{H}$, the common practice in computational electromagnetics, rather than $\mathbf{E}$ and $\mathbf{B}$, which is common in other fields.)

## Maxwell's equations in 1D

Consider a one-dimensional space where there are only variations in
the $x$ direction.  Assume that the electric field only has a $z$
component.  In this case Faraday's law can be written

\begin{equation}
  -\mu \frac{\partial \mathbf{H}}{\partial t} =
  \nabla \times \mathbf{E} =
  \left|
  \begin{array}{ccc}
     \hat{\mathbf{a}}_{x} & \hat{\mathbf{a}}_{y} & \hat{\mathbf{a}}_{z} \\
     \frac{\partial}{\partial x} & 0 & 0 \\
     0 & 0 & E_z
  \end{array}
  \right|
  = -\hat{\mathbf{a}}_{y}\frac{\partial E_z}{\partial x}.
 \label{eq:faradayVec}
\end{equation}

Thus $H_y$ must be the only non-zero component of the magnetic field
which is time varying.  (Since the right-hand side of this equation
has only a $y$ component, the magnetic field may have non-zero
components in the $x$ and $z$ directions, but they must be static.  We
will not be concerned with static fields here.)  Knowing this,
Ampere's law can be written

\begin{equation}
  \epsilon \frac{\partial \mathbf{E}}{\partial t} =
  \nabla \times \mathbf{H} =
  \left|
  \begin{array}{ccc}
     \hat{\mathbf{a}}_{z} & \hat{\mathbf{a}}_{y} & \hat{\mathbf{a}}_{z} \\
     \frac{\partial}{\partial x} & 0 & 0 \\
     0 & H_y & 0
  \end{array}
  \right|
  = \hat{\mathbf{a}}_{z}\frac{\partial H_y}{\partial x}.
 \label{eq:ampereVec}
\end{equation}


The two scalar equations obtained from the previous equations are

\begin{eqnarray}
  \mu \frac{\partial H_y}{\partial t} &=&
    \frac{\partial E_z}{\partial x},  \label{eq:faradayScalar}
  \\
  \epsilon \frac{\partial E_z}{\partial t} &=&
    \frac{\partial H_y}{\partial x}. \label{eq:ampereScalar}
\end{eqnarray}

The first equation gives the temporal derivative of the magnetic field
in terms of the spatial derivative of the electric field.  Conversely,
the second equation gives the temporal derivative of the electric
field in terms of the spatial derivative of the magnetic field.  As
will be shown, the first equation will be used to advance the magnetic
field in time while the second will be used to advance the electric
field.  A method in which one field is advanced and then the other,
and then the process is repeated, is known as a leap-frog method.


## Central differencing

Consider the
Taylor series expansions of the function $f(x)$ expanded about the
point $x_0$ with an offset of $\pm\delta/2$:

\begin{eqnarray}
  f\!\left(x_0+\frac{\delta}{2}\right) &=&
    f(x_0) + \frac{\delta}{2} f'(x_0) +
        \frac{1}{2!}\left(\frac{\delta}{2}\right)^2 f''(x_0) + 
        \frac{1}{3!}\left(\frac{\delta}{2}\right)^3 f'''(x_0) + \ldots,
  \\
  f\!\left(x_0-\frac{\delta}{2}\right) &=&
    f(x_0) - \frac{\delta}{2} f'(x_0) +
        \frac{1}{2!}\left(\frac{\delta}{2}\right)^2 f''(x_0) -
        \frac{1}{3!}\left(\frac{\delta}{2}\right)^3 f'''(x_0) + \ldots
\end{eqnarray}

where the primes indicate differentiation.

Subtracting the second
equation from the first yields

\begin{equation}
  f\!\left(x_0+\frac{\delta}{2}\right) -
  f\!\left(x_0-\frac{\delta}{2}\right) = 
   \delta f'(x_0) +
        \frac{2}{3!}\left(\frac{\delta}{2}\right)^3 f'''(x_0) + \ldots
\end{equation}


Dividing by $\delta$ produces

\begin{equation}
  \frac{f\!\left(x_0+\frac{\delta}{2}\right) -
        f\!\left(x_0-\frac{\delta}{2}\right)}{\delta} = 
   f'(x_0) +
        \frac{1}{3!}\frac{\delta^2}{2^2} f'''(x_0) + \ldots
\end{equation}


Thus the term on the left is equal to the derivative of the function
at the point $x_0$ plus a term which depends on $\delta^2$ plus an
infinite number of other terms which are not shown.  For the terms
which are not shown, the next would depend on $\delta^4$ and all
subsequent terms would depend on even higher powers of $\delta$.


Rearranging slightly, this relationship is often stated as

\begin{equation}
  \left.\frac{df(x)}{dx}\right|_{x=x_0} = 
  \frac{f\!\left(x_0+\frac{\delta}{2}\right) -
        f\!\left(x_0-\frac{\delta}{2}\right)}{\delta} + O(\delta^2).
\end{equation}

The "big-Oh" term represents all the terms that are not explicitly
shown and the value in parentheses, i.e., $\delta^2$, indicates the
lowest order of $\delta$ in these hidden terms.  If $\delta$ is
sufficiently small, a reasonable approximation to the derivative may
be obtained by simply neglecting all the terms represented by the
"big-Oh" term.


Thus, the central-difference approximation is given
by

\begin{equation}
  \left.\frac{df(x)}{dx}\right|_{x=x_0} \approx
  \frac{f\!\left(x_0+\frac{\delta}{2}\right) -
        f\!\left(x_0-\frac{\delta}{2}\right)}{\delta}.
\end{equation}

Note that the central difference provides an approximation of the
derivative of the function at $x_0$, but the function is not actually
sampled there.  Instead, the function is sampled at the neighboring
points $x_0+\delta/2$ and $x_0-\delta/2$.  Since the lowest power of
$\delta$ being ignored is second order, the central difference is said
to have second-order accuracy or second-order behavior.  This implies
that if $\delta$ is reduced by a factor of $10$, the error in the
approximation should be reduced by a factor of $100$ (at least
approximately).  In the limit as $\delta$ goes to zero, the
approximation becomes exact.


## Applying central differing
The next step is to replace the derivatives in
the previous equations with finite
differences.  To do this, space and time need to be discretized.  The
following notation will be used to indicate the location where the
fields are sampled in space and time

\begin{eqnarray}
  E_z(x,t) \!&=&\! E_z(m\Delta_x, q\Delta_t) = E_z^{q}\!\left[m\right], \\
  H_y(x,t) \!&=&\! H_y(m\Delta_x, q\Delta_t) = H_y^{q}\!\left[m\right],
\end{eqnarray}

where $\Delta_x$ is the spatial offset between sample points and $\Delta_t$
is the temporal offset.  The index $m$ corresponds to the spatial
step, effectively the spatial location, while the index $q$
corresponds to the temporal step.  When written as a superscript $q$
still represents the temporal step---it is not an exponent.  When
implementing FDTD algorithms we will see that the spatial indices are
used as array indices while the temporal index, which is essentially a
global parameter, is not explicitly specified for each field location.
Hence, it is reasonable to keep the spatial indices as an explicit
argument while indicating the temporal index separately.

## Gridding
Although we only have one spatial dimension, time can be thought of as another dimension. Thus this is effectively a form of two-dimensional problem. The question now is: How should the electric and magnetic field sample points, also known as nodes, be arranged in space and time? The answer is shown in the figure below. The electric-field nodes are shown as circles and the magnetic- field nodes as triangles. Assume that all the fields below the dashed line are known—they are considered to be in the past—while the fields above the dashed line are future fields and hence unknown. The FDTD algorithm provides a way to obtain the future fields from the past fields.

<figure>
  
  <figcaption>The arrangement of electric- and magnetic-field nodes in space and time. The electric- field nodes are shown as circles and the magnetic-field nodes as triangles. The indicated point is where the difference equation is expanded to obtain an update equation for $H_y$.</figcaption>
</figure>


Consider Faraday's
law at the space-time point $((m+1/2)\Delta_x,q\Delta_t)$

\begin{equation}
  \left.\mu
  \frac{\partial H_y}{\partial t}\right|_{(m+1/2)\Delta_x,q\Delta_t}
  =
  \left.\frac{\partial E_z}{\partial x}\right|_{(m+1/2)\Delta_x,q\Delta_t}.
  \label{eq:faradayDiscrete}
\end{equation}


The temporal derivative is replaced by a finite difference involving
$H_y^{q+\frac{1}{2}}\!\left[m+\frac{1}{2}\right]$ and $H_y^{q-\frac{1}{2}}\!\left[m+\frac{1}{2}\right]$
(i.e., the magnetic field at a fixed location but two different times)
while the spatial derivative is replaced by a finite difference
involving $E_z^{q}\!\left[m+1\right]$ and $E_z^{q}\!\left[m\right]$ (i.e., the
electric field at two different locations but one time).  This yields

\begin{equation}
  \mu\frac{H_y^{q+\frac{1}{2}}\!\!\left[m+\frac{1}{2}\right] -
           H_y^{q-\frac{1}{2}}\!\!\left[m+\frac{1}{2}\right]} {\Delta_t} = 
  \frac{E_z^{q}\!\left[m+1\right] - E_z^{q}\!\left[m\right]}{\Delta_x}.
  \label{eq:faradayFdtd1D}
\end{equation}


Solving this for $H_y^{q+\frac{1}{2}}\!\!\left[m+\frac{1}{2}\right]$ yields:

\begin{equation}
  H_y^{q+\frac{1}{2}}\!\!\left[m+\frac{1}{2}\right] = H_y^{q-\frac{1}{2}}\!\!\left[m+\frac{1}{2}\right] +
  \frac{\Delta_t}{\mu\Delta_x}
  \left(E_z^{q}\!\left[m+1\right] - E_z^{q}\!\left[m\right]\right).
  \label{eq:updateHy}
\end{equation}

This is known as an update equation, specifically the update equation
for the $H_y$ field.  It is a generic equation which can be applied to
any magnetic-field node.  It shows that the future value of $H_y$
depends on only its previous value and the neighboring electric
fields.

After applying this equation to all the magnetic-field
nodes, the dividing line between future and past values has advanced a
half time-step.  The space-time grid thus appears as shown in the figure below which is identical to the figure above  except for the advancement of the past/future
dividing line.

<figure>
  
  <figcaption>Space-time after updating the magnetic field. The dividing line between future and past values has moved forward a half temporal step. The indicated point is where the difference equation is written to obtain an update equation for $E_z$.</figcaption>
</figure>



Starting with Ampere's law, a smilar update equation for the $E_z$ field can be derived:

\begin{equation}
  E_z^{q+1}\!\left[m\right] = E_z^{q}\!\left[m\right] + 
  \frac{\Delta_t}{\epsilon\Delta_x}
   \left(H_y^{q+\frac{1}{2}}\!\!\left[m+\frac{1}{2}\right] - H_y^{q+\frac{1}{2}}\!\!\left[m-\frac{1}{2}\right]\right).
\end{equation}

The indices in this equation are generic so that the same
equation holds for every $E_z$ node.  Similar to the update equation for
the magnetic field, here we see that the future value of $E_z$ depends
on only its past value and the value of the neighboring magnetic
fields.

After applying this equation to every electric-field node in the
grid, the dividing line between what is known and what is unknown
moves forward another one-half temporal step.  One is essentially back
to the situation depicted in the first figure: the future
fields closest to the dividing line between the future and past are
magnetics fields.  They would be updated again, then the electric
fields would be updated, and so on.

## Turning this into a computer program

You will have noticed indices are either an integer ($m$) or an integer plus one half ($m+\frac{1}{2}$). While we coud use these indices in a list in some languages, common practice is to simply to use integers.  If we take the figure above, and compress the time axis so that we just view the space axis, we can pick a standard way to express the indices of the $E_z$ and $H_y$ fields as follows (note that we also need to choose variable names that don't involve subscripts).

<figure>
  
  <figcaption>A one-dimensional FDTD space showing the assumed spatial arrangement of the electric- and magnetic-field nodes in the arrays ez and hy. Note that an electric-field node is assumed to exist to the left of the magnetic-field node with the same index.</figcaption>
</figure>

And the update equations then are:

```  hy[m] = hy[m] + (ez[m + 1] - ez[m]) / imp0```

and

```  ez[m] = ez[m] + (hy[m] - hy[m - 1]) * imp0```

```imp0``` is the impedance of free space:
$imp_0 = \frac{E_z}{H_y} = \mu_0 c_0 = \sqrt{\frac{\mu_0}{\epsilon_0}} = \frac{1}{\epsilon_0 c_0}$


where 
$\mu_0$ is the magnetic constant,
$\epsilon_0$ is the electric constant,
and $c_0$ is the speed of light in free space.


```python

```

529cb1ff7c3718997ae8663a038381bfa51b1abf
ipynb
Jupyter Notebook
fdtd_intro.ipynb
danielskatz/repro-fdtd1d
e960dcac74147de1fd9cd1ffb7ccb98e4e5e0f00
MIT
2020-08-13T19:12:01.000Z
2022-01-12T21:41:33.000Z
fdtd_intro.ipynb
danielskatz/repro-fdtd1d
e960dcac74147de1fd9cd1ffb7ccb98e4e5e0f00
MIT
fdtd_intro.ipynb
danielskatz/repro-fdtd1d
e960dcac74147de1fd9cd1ffb7ccb98e4e5e0f00
MIT
2018-07-29T12:04:40.000Z
2020-11-25T11:11:24.000Z
Qwen/Qwen-72B
1. YES
2. YES