FDTD Derivation

1. Introduction

This page gives the detailed mathematical derivation of the main finite-difference time-domain (FDTD) update equations used to solve Maxwell’s equations in the time domain. The purpose of this page is different from the main theory page. The aim here is not to provide a lightweight overview, but to show explicitly how the continuous curl equations are converted into the discrete algebraic expressions advanced on a computer. This is also the part of the FDTD method which is easiest to get wrong in code: sign errors, misplaced indices, or inconsistent time staggering will immediately produce unstable or unphysical results.

We therefore begin from Ampère’s law and Faraday’s law in differential form, expand them into Cartesian components, and then specialise to a two-dimensional optical problem. The spatial derivatives are replaced by centred finite differences, the time derivatives are discretised using the standard leapfrog scheme, and the resulting expressions are rearranged to give explicit update formulae for the electric and magnetic field components.

To keep the derivation readable, we reduce the equations here to a two-dimensional case in which \(\partial/\partial y = 0\) and the magnetic conductivity is set to zero, \(\sigma_m = 0\). This makes the algebra much easier to follow and also shows explicitly how the dimensional reduction is carried out in practice. In the full OghmaNano FDTD solver, however, the electromagnetic equations are solved in three dimensions. The two-dimensional derivation given here is therefore intended as a clear worked example of the method, rather than a statement that the underlying solver is limited to 2D.

Throughout this page, electric field components are evaluated at integer time steps \(n\Delta t\), while magnetic field components are evaluated at half-integer time steps \(\left(n+\tfrac{1}{2}\right)\Delta t\). This is the standard Yee-grid leapfrog arrangement used in FDTD. Unlike the shorter overview page, the indexing is stated more explicitly here so that the derivation can be read as an implementation guide as well as a mathematical note.

2. Starting equations

We start from Maxwell’s curl equations in a linear medium with electric conductivity \(\sigma\), electric permittivity \(\epsilon\), magnetic permeability \(\mu\), and, for generality, magnetic conductivity \(\sigma_m\):

\[ \nabla \times \boldsymbol{H} = \sigma \boldsymbol{E} + \epsilon \frac{\partial \boldsymbol{E}}{\partial t}, \qquad \nabla \times \boldsymbol{E} = - \sigma_m \boldsymbol{H} - \mu \frac{\partial \boldsymbol{H}}{\partial t}. \]

In determinant form, the curls are

\[ \nabla \times \boldsymbol{H} = \begin{vmatrix} \hat{\boldsymbol{x}} & \hat{\boldsymbol{y}} & \hat{\boldsymbol{z}} \\ \dfrac{\partial}{\partial x} & \dfrac{\partial}{\partial y} & \dfrac{\partial}{\partial z} \\ H_x & H_y & H_z \end{vmatrix}, \qquad \nabla \times \boldsymbol{E} = \begin{vmatrix} \hat{\boldsymbol{x}} & \hat{\boldsymbol{y}} & \hat{\boldsymbol{z}} \\ \dfrac{\partial}{\partial x} & \dfrac{\partial}{\partial y} & \dfrac{\partial}{\partial z} \\ E_x & E_y & E_z \end{vmatrix}. \]

These are the continuous-time, continuous-space equations. The purpose of the remainder of the page is to convert them into an explicit finite-difference scheme that can be stepped forward in time.

3. Component form of the curl equations

Expanding Ampère’s law into Cartesian components gives

\[ \sigma E_x + \epsilon \frac{\partial E_x}{\partial t} = \frac{\partial H_z}{\partial y} - \frac{\partial H_y}{\partial z}, \]

\[ \sigma E_y + \epsilon \frac{\partial E_y}{\partial t} = \frac{\partial H_x}{\partial z} - \frac{\partial H_z}{\partial x}, \]

\[ \sigma E_z + \epsilon \frac{\partial E_z}{\partial t} = \frac{\partial H_y}{\partial x} - \frac{\partial H_x}{\partial y}. \]

Expanding Faraday’s law gives

\[ -\sigma_m H_x - \mu \frac{\partial H_x}{\partial t} = \frac{\partial E_z}{\partial y} - \frac{\partial E_y}{\partial z}, \]

\[ -\sigma_m H_y - \mu \frac{\partial H_y}{\partial t} = \frac{\partial E_x}{\partial z} - \frac{\partial E_z}{\partial x}, \]

\[ -\sigma_m H_z - \mu \frac{\partial H_z}{\partial t} = \frac{\partial E_y}{\partial x} - \frac{\partial E_x}{\partial y}. \]

At this stage nothing has yet been assumed about dimensionality. The expressions above are simply the component form of the full three-dimensional curl equations.

4. Two-dimensional reduction

To make the derivation tractable, we now reduce the problem to two dimensions by assuming there is no variation in the \(y\)-direction, so \(\partial / \partial y = 0\). This step is taken purely to simplify the algebra and make the structure of the update equations easier to follow. The same procedure can be extended directly to three dimensions, and in practice OghmaNano solves the full 3D form of Maxwell’s equations.

With this assumption, the electric-field equations reduce to

\[ \sigma E_x + \epsilon \frac{\partial E_x}{\partial t} = - \frac{\partial H_y}{\partial z}, \]

\[ \sigma E_y + \epsilon \frac{\partial E_y}{\partial t} = - \frac{\partial H_z}{\partial x} + \frac{\partial H_x}{\partial z}, \]

\[ \sigma E_z + \epsilon \frac{\partial E_z}{\partial t} = \frac{\partial H_y}{\partial x}. \]

If we also set the magnetic conductivity to zero, \(\sigma_m = 0\), the magnetic-field equations become

\[ \mu \frac{\partial H_x}{\partial t} = \frac{\partial E_y}{\partial z}, \]

\[ \mu \frac{\partial H_y}{\partial t} = \frac{\partial E_z}{\partial x} - \frac{\partial E_x}{\partial z}, \]

\[ \mu \frac{\partial H_z}{\partial t} = - \frac{\partial E_y}{\partial x}. \]

These are the continuous equations that will now be discretised. It is worth stating explicitly that this is a 2D reduction of the full 3D problem. The purpose is not to claim that all FDTD is intrinsically two-dimensional, but to derive clean update equations for the common case where one direction is suppressed.

5. Yee grid, indexing, and notation

To solve the above equations numerically, both space and time are discretised. Let the computational grid spacings in the \(x\)- and \(z\)-directions be \(\Delta x\) and \(\Delta z\), respectively, and let the time step be \(\Delta t\). We use the integer indices \(i\) and \(k\) for the \(x\)- and \(z\)-directions:

\[ x_i = i \Delta x, \qquad z_k = k \Delta z. \]

Time is indexed by the integer \(n\), so that

\[ t^n = n \Delta t. \]

In the Yee scheme, the electric and magnetic fields are not stored at the same locations in either space or time. Instead, they are staggered. Electric-field components are evaluated at integer time levels \(n\Delta t\), while magnetic-field components are evaluated at half-integer time levels \(\left(n+\tfrac{1}{2}\right)\Delta t\). This is the origin of the leapfrog structure of FDTD.

Spatially, the field components are also offset so that each centred finite-difference curl is naturally defined between neighbouring grid locations. In a full code implementation, each field component carries its own pair of spatial indices, for example \(E_x^n(i,k)\), \(H_y^{n+\frac{1}{2}}(i+\tfrac{1}{2},k+\tfrac{1}{2})\), and so on. To keep the algebra readable, some of the repeated “spectator” indices are not always written explicitly in the steps below. However, the interpretation is always the same:

\(i\) indexes the \(x\)-direction,
\(k\) indexes the \(z\)-direction,
\(n\) indexes integer electric-field time levels, and
\(n+\tfrac{1}{2}\) indexes half-step magnetic-field time levels.

For example, the centred finite-difference approximation to \(\partial H_y / \partial z\) is

\[ \frac{\partial H_y}{\partial z} \approx \frac{ H_y^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_y^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z}, \]

and similarly the centred approximation to \(\partial H_z / \partial x\) is

\[ \frac{\partial H_z}{\partial x} \approx \frac{ H_z^{n+\frac{1}{2}}\!\left(i+\frac{1}{2},k\right) - H_z^{n+\frac{1}{2}}\!\left(i-\frac{1}{2},k\right) }{\Delta x}. \]

With this notation made explicit, the derivation below can be read unambiguously as a proper Yee-grid update scheme.

6. Derivation of the electric-field update equations

The electric field is evaluated at integer time steps \(n\Delta t\), while the magnetic field is evaluated at half-integer time steps \(\left(n+\tfrac{1}{2}\right)\Delta t\). To derive the update equations, the time derivative is approximated by a centred finite difference, and the conductive term is evaluated using the average of the old and new electric fields. This gives a numerically stable trapezoidal treatment of conductivity.

This averaging step is important. If the conductivity term were treated inconsistently, the update coefficients would no longer have the standard damped form, and the correspondence with the usual FDTD equations would be lost.

6.1 Update equation for \(E_x\)

Start from

\[ \sigma E_x + \epsilon \frac{\partial E_x}{\partial t} = - \frac{\partial H_y}{\partial z}. \]

Evaluate the conductivity term using the average of \(E_x^{n+1}\) and \(E_x^n\), and evaluate the time derivative using a centred finite difference:

\[ \sigma \frac{E_x^{n+1}(i,k) + E_x^n(i,k)}{2} + \epsilon \frac{E_x^{n+1}(i,k) - E_x^n(i,k)}{\Delta t} = - \frac{ H_y^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_y^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z}. \]

Collect the \(E_x^{n+1}\) terms on the left:

\[ \left( \frac{\sigma}{2} + \frac{\epsilon}{\Delta t} \right) E_x^{n+1}(i,k) = - \frac{ H_y^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_y^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z} + \left( \frac{\epsilon}{\Delta t} - \frac{\sigma}{2} \right) E_x^n(i,k). \]

Therefore

\[ E_x^{n+1}(i,k) = \frac{ \left( \dfrac{\epsilon}{\Delta t} - \dfrac{\sigma}{2} \right) E_x^n(i,k) - \dfrac{ H_y^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_y^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z} }{ \dfrac{\epsilon}{\Delta t} + \dfrac{\sigma}{2} }. \]

Multiplying numerator and denominator by \(2\Delta t\) gives the familiar coefficient form:

\[ E_x^{n+1}(i,k) = \frac{ 2\epsilon - \sigma \Delta t }{ 2\epsilon + \sigma \Delta t } E_x^n(i,k) - \frac{ 2 \Delta t }{ 2\epsilon + \sigma \Delta t } \frac{ H_y^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_y^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z}. \]

6.2 Update equation for \(E_y\)

Start from

\[ \sigma E_y + \epsilon \frac{\partial E_y}{\partial t} = - \frac{\partial H_z}{\partial x} + \frac{\partial H_x}{\partial z}. \]

Discretise both derivative terms on the Yee grid:

\[ \sigma \frac{E_y^{n+1}(i,k) + E_y^n(i,k)}{2} + \epsilon \frac{E_y^{n+1}(i,k) - E_y^n(i,k)}{\Delta t} = - \frac{ H_z^{n+\frac{1}{2}}\!\left(i+\frac{1}{2},k\right) - H_z^{n+\frac{1}{2}}\!\left(i-\frac{1}{2},k\right) }{\Delta x} + \frac{ H_x^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_x^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z}. \]

Rearranging gives

\[ E_y^{n+1}(i,k) = \frac{ 2\epsilon - \sigma \Delta t }{ 2\epsilon + \sigma \Delta t } E_y^n(i,k) + \frac{ 2 \Delta t }{ 2\epsilon + \sigma \Delta t } \left[ - \frac{ H_z^{n+\frac{1}{2}}\!\left(i+\frac{1}{2},k\right) - H_z^{n+\frac{1}{2}}\!\left(i-\frac{1}{2},k\right) }{\Delta x} + \frac{ H_x^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_x^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z} \right]. \]

This is the only electric-field component in the 2D reduction which contains contributions from two separate magnetic-field derivatives, and it is therefore the easiest one to get wrong in code. In particular, the relative sign between the \(x\)- and \(z\)-derivative terms follows directly from the curl expansion and should not be altered.

6.3 Update equation for \(E_z\)

Start from

\[ \sigma E_z + \epsilon \frac{\partial E_z}{\partial t} = \frac{\partial H_y}{\partial x}. \]

Discretisation gives

\[ \sigma \frac{E_z^{n+1}(i,k) + E_z^n(i,k)}{2} + \epsilon \frac{E_z^{n+1}(i,k) - E_z^n(i,k)}{\Delta t} = \frac{ H_y^{n+\frac{1}{2}}\!\left(i+\frac{1}{2},k\right) - H_y^{n+\frac{1}{2}}\!\left(i-\frac{1}{2},k\right) }{\Delta x}. \]

Hence

\[ E_z^{n+1}(i,k) = \frac{ 2\epsilon - \sigma \Delta t }{ 2\epsilon + \sigma \Delta t } E_z^n(i,k) + \frac{ 2 \Delta t }{ 2\epsilon + \sigma \Delta t } \frac{ H_y^{n+\frac{1}{2}}\!\left(i+\frac{1}{2},k\right) - H_y^{n+\frac{1}{2}}\!\left(i-\frac{1}{2},k\right) }{\Delta x}. \]

At this point the discrete electric-field update equations have been obtained in their standard damped form. The next step is to derive the corresponding magnetic-field equations advanced at half-integer time levels.

7. Derivation of the magnetic-field update equations

The magnetic field is updated at half-integer time steps. With \(\sigma_m = 0\), the magnetic-field equations contain no damping term, so the update formulae take the simpler forward form. This part of the derivation is algebraically shorter than the electric-field case because there is no trapezoidal conductivity term to carry through.

7.1 Update equation for \(H_x\)

Start from

\[ \mu \frac{\partial H_x}{\partial t} = \frac{\partial E_y}{\partial z}. \]

Discretising in time between the half steps \(n-\tfrac{1}{2}\) and \(n+\tfrac{1}{2}\), and in space using a centred difference in the \(z\)-direction, gives

\[ \mu \frac{ H_x^{n+\frac{1}{2}}(i,k+\frac{1}{2}) - H_x^{n-\frac{1}{2}}(i,k+\frac{1}{2}) }{\Delta t} = \frac{ E_y^n(i,k+1) - E_y^n(i,k) }{\Delta z}. \]

In the compact half-index notation often used in derivations, the same result can be written as

\[ \mu \frac{ H_x^{n+\frac{1}{2}} - H_x^{n-\frac{1}{2}} }{\Delta t} = \frac{ E_y^n\!\left(k+\frac{1}{2}\right) - E_y^n\!\left(k-\frac{1}{2}\right) }{\Delta z}. \]

Therefore

\[ H_x^{n+\frac{1}{2}} = H_x^{n-\frac{1}{2}} + \frac{\Delta t}{\mu} \frac{ E_y^n\!\left(k+\frac{1}{2}\right) - E_y^n\!\left(k-\frac{1}{2}\right) }{\Delta z}. \]

7.2 Update equation for \(H_y\)

Start from

\[ \mu \frac{\partial H_y}{\partial t} = \frac{\partial E_z}{\partial x} - \frac{\partial E_x}{\partial z}. \]

Discretisation gives

\[ \mu \frac{ H_y^{n+\frac{1}{2}}(i+\frac{1}{2},k+\frac{1}{2}) - H_y^{n-\frac{1}{2}}(i+\frac{1}{2},k+\frac{1}{2}) }{\Delta t} = \frac{ E_z^n(i+1,k+\frac{1}{2}) - E_z^n(i,k+\frac{1}{2}) }{\Delta x} - \frac{ E_x^n(i+\frac{1}{2},k+1) - E_x^n(i+\frac{1}{2},k) }{\Delta z}. \]

Again, in the shortened half-index notation,

\[ \mu \frac{ H_y^{n+\frac{1}{2}} - H_y^{n-\frac{1}{2}} }{\Delta t} = \frac{ E_z^n\!\left(i+\frac{1}{2}\right) - E_z^n\!\left(i-\frac{1}{2}\right) }{\Delta x} - \frac{ E_x^n\!\left(k+\frac{1}{2}\right) - E_x^n\!\left(k-\frac{1}{2}\right) }{\Delta z}. \]

Hence

\[ H_y^{n+\frac{1}{2}} = H_y^{n-\frac{1}{2}} + \frac{\Delta t}{\mu} \left[ \frac{ E_z^n\!\left(i+\frac{1}{2}\right) - E_z^n\!\left(i-\frac{1}{2}\right) }{\Delta x} - \frac{ E_x^n\!\left(k+\frac{1}{2}\right) - E_x^n\!\left(k-\frac{1}{2}\right) }{\Delta z} \right]. \]

7.3 Update equation for \(H_z\)

Start from

\[ \mu \frac{\partial H_z}{\partial t} = - \frac{\partial E_y}{\partial x}. \]

Discretising gives

\[ \mu \frac{ H_z^{n+\frac{1}{2}}(i+\frac{1}{2},k) - H_z^{n-\frac{1}{2}}(i+\frac{1}{2},k) }{\Delta t} = - \frac{ E_y^n(i+1,k) - E_y^n(i,k) }{\Delta x}. \]

Equivalently, in compact notation,

\[ \mu \frac{ H_z^{n+\frac{1}{2}} - H_z^{n-\frac{1}{2}} }{\Delta t} = - \frac{ E_y^n\!\left(i+\frac{1}{2}\right) - E_y^n\!\left(i-\frac{1}{2}\right) }{\Delta x}. \]

Therefore

\[ H_z^{n+\frac{1}{2}} = H_z^{n-\frac{1}{2}} - \frac{\Delta t}{\mu} \frac{ E_y^n\!\left(i+\frac{1}{2}\right) - E_y^n\!\left(i-\frac{1}{2}\right) }{\Delta x}. \]

With these substitutions complete, the standard leapfrog FDTD update equations follow directly.

8. Final update equations

Collecting the results, the two-dimensional FDTD update equations are

\[ H_x^{n+\frac{1}{2}}(i,k+\tfrac{1}{2}) = H_x^{n-\frac{1}{2}}(i,k+\tfrac{1}{2}) + \frac{\Delta t}{\mu} \frac{ E_y^n(i,k+1) - E_y^n(i,k) }{\Delta z}, \]

\[ H_y^{n+\frac{1}{2}}(i+\tfrac{1}{2},k+\tfrac{1}{2}) = H_y^{n-\frac{1}{2}}(i+\tfrac{1}{2},k+\tfrac{1}{2}) + \frac{\Delta t}{\mu} \left[ \frac{ E_z^n(i+1,k+\tfrac{1}{2}) - E_z^n(i,k+\tfrac{1}{2}) }{\Delta x} - \frac{ E_x^n(i+\tfrac{1}{2},k+1) - E_x^n(i+\tfrac{1}{2},k) }{\Delta z} \right], \]

\[ H_z^{n+\frac{1}{2}}(i+\tfrac{1}{2},k) = H_z^{n-\frac{1}{2}}(i+\tfrac{1}{2},k) - \frac{\Delta t}{\mu} \frac{ E_y^n(i+1,k) - E_y^n(i,k) }{\Delta x}. \]

Written in this form, the meaning of the leapfrog scheme is explicit. The electric field at time level \(n+1\) is obtained from the magnetic field at time level \(n+\tfrac{1}{2}\), while the magnetic field at time level \(n+\tfrac{1}{2}\) is obtained from the electric field at time level \(n\). The fields therefore advance in an alternating sequence in time.

9. Remarks on implementation

In a practical code, the coefficients multiplying the old field values and the curl terms are often precomputed and stored, especially when \(\epsilon\), \(\mu\), or \(\sigma\) vary from cell to cell. For example, the electric-field update commonly takes the form

\[ E^{n+1} = C_a E^n + C_b (\nabla \times H), \]

where the local coefficients \(C_a\) and \(C_b\) encode the material parameters and the time step. This makes the implementation cleaner and avoids recalculating the same factors every time step.

It is also worth emphasising that the derivation above gives only the core FDTD update scheme. A usable simulation additionally requires sources, boundary conditions, absorbing layers such as perfectly matched layers (PML), and a time step satisfying the Courant stability condition. None of those additions changes the fundamental leapfrog structure derived here, but each of them modifies how the update equations are applied at particular grid locations.

Finally, although the notation above has been written carefully enough to be implementation-relevant, the exact spatial placement of each component depends on the grid convention used in the code. The crucial requirement is not the typography itself, but consistency: the curl operator, the staggering, and the signs must all correspond to the same Yee-cell arrangement throughout.