← Back to homepage

Engineering tutorial

Computational Electromagnetics for Engineers

FDTD and FEM for RF structures, photonics, metamaterials, radar absorbers, and open-region scattering.

This page is organized around four things: governing equations, physical quantities with units, solver choice, and extraction of defensible observables.

Page focus

  • Maxwell equations in usable engineering form
  • Field, material, and observable quantities with SI units
  • 1D to 2D progression for self-study
  • Open-region modeling with source injection and CPML

Contents

1. Governing equations & quantities 2. Engineering observables 3. FDTD and FEM selection 4. 1D FDTD 5. 1D FEM 6. 2D FDTD 7. 2D FEM 8. Workflow and validation 9. Sources and boundaries 10. Example codes

1. Governing equations and physical quantities

Computational electromagnetics starts from Maxwell’s equations. The solver changes. The physics does not.

\[ \nabla \times \mathbf{E} = -\frac{\partial \mathbf{B}}{\partial t}, \qquad \nabla \times \mathbf{H} = \mathbf{J} + \frac{\partial \mathbf{D}}{\partial t} \] \[ \nabla \cdot \mathbf{D} = \rho, \qquad \nabla \cdot \mathbf{B} = 0 \]

For linear isotropic media:

\[ \mathbf{D} = \varepsilon \mathbf{E}, \qquad \mathbf{B} = \mu \mathbf{H}, \qquad \mathbf{J} = \sigma \mathbf{E} \]
In practical RF and photonics work, \(\varepsilon\), \(\mu\), and \(\sigma\) may vary with position and may also depend on frequency. That is the point where numerical methods become necessary.

Use SI units unless a project convention states otherwise. In optics, wavelength is often reported in nm; in RF, frequency is often reported in MHz or GHz. The underlying SI units remain the reference.

Core field and material quantities

Symbol Name SI unit
\(\mathbf{E}\)Electric fieldV/m
\(\mathbf{D}\)Electric flux densityC/m²
\(\mathbf{H}\)Magnetic fieldA/m
\(\mathbf{B}\)Magnetic flux densityT
\(\mathbf{J}\)Current densityA/m²
\(\rho\)Charge densityC/m³
\(\varepsilon\)PermittivityF/m
\(\mu\)PermeabilityH/m
\(\sigma\)ConductivityS/m
\(\varepsilon_r\)Relative permittivitydimensionless
\(\mu_r\)Relative permeabilitydimensionless

Wave and observable quantities

Symbol Name SI unit
\(f\)FrequencyHz
\(\omega\)Angular frequencyrad/s
\(\lambda\)Wavelengthm
\(k_0\)Free-space wavenumberrad/m
\(\eta\)Wave impedanceΩ
\(\mathbf{S}=\mathbf{E}\times\mathbf{H}\)Poynting vectorW/m²
\(R\)Reflectancedimensionless or %
\(T\)Transmittancedimensionless or %
\(A\)Absorptancedimensionless or %
\(S_{11}, S_{12}, S_{21}, S_{22}\)S-parameterscomplex ratio; often magnitude, phase, or dB
\(Q\)Quality factordimensionless

FDTD in context. FDTD traces back to Yee’s 1966 staggered-grid formulation and became a practical mainstream tool after absorbing boundary treatments such as PML made open-region simulations reliable. Its historical strength is direct time-domain treatment of wave propagation, transients, and broadband scattering.

FEM in context. FEM grew out of variational and structural analysis methods and was established earlier in computational mechanics before becoming a standard electromagnetic tool. Its historical advantage is the weak-form framework, which naturally supports irregular geometries, boundary-value problems, and later multiphysics coupling.

2. Engineering observables

Use observables that connect directly to device performance.

Reflection / transmission

For slabs, coatings, filters, absorbers, and open-region scattering.

\[ A = 1 - R - T \]
  • \(R\): reflected power ratio. Unit: dimensionless or %.
  • \(T\): transmitted power ratio. Unit: dimensionless or %.
  • \(A\): absorbed power ratio. Unit: dimensionless or %.

S-parameters

For ports, feeds, antennas, RF packaging, and microwave components.

  • \(S_{11}\): reflection seen at port 1. Unit: dimensionless complex ratio, often shown in dB.
  • \(S_{21}\): transmission from port 1 to port 2. Unit: dimensionless complex ratio, often shown in dB.
  • \(S_{12}\): transmission from port 2 to port 1. Unit: dimensionless complex ratio, often shown in dB.
  • \(S_{22}\): reflection seen at port 2. Unit: dimensionless complex ratio, often shown in dB.

Field quantities

For hot spots, mode shape, localization, and power flow.

  • \(\mathbf{E}\): electric field. Unit: V/m.
  • \(\mathbf{H}\): magnetic field. Unit: A/m.
  • \(|\mathbf{E}|^2\): electric-field intensity proxy. Unit: V²/m².
  • \(|\mathbf{H}|^2\): magnetic-field intensity proxy. Unit: A²/m².
  • \(\langle \mathbf{S} \rangle\): time-averaged power-flow density. Unit: W/m².
Field plots help diagnose behavior. Final design decisions should be made from normalized observables.

3. FDTD and FEM selection

Method Formulation Main strength Main limitation Typical use
FDTD Time domain, staggered grid Broadband response from one pulse; good physical intuition Time stepping cost; CFL limit; structured grid Transient scattering, wideband spectra, source propagation
FEM Frequency domain, weak form Handles irregular geometry and boundary-value problems cleanly Broadband response requires repeated frequency solves Resonators, bounded domains, complex geometry, multiphysics coupling

4. 1D finite-difference time-domain (FDTD)

Method idea
Electric and magnetic fields are staggered in space and time. The update is explicit and physically transparent.
Update equations
\[ H_{i+\frac{1}{2}}^{n+\frac{1}{2}} = H_{i+\frac{1}{2}}^{n-\frac{1}{2}} + \frac{\Delta t}{\mu \Delta x}\left(E_{i+1}^{n} - E_i^{n}\right) \] \[ E_i^{n+1} = E_i^n + \frac{\Delta t}{\varepsilon_i \Delta x}\left(H_{i+\frac{1}{2}}^{n+\frac{1}{2}} - H_{i-\frac{1}{2}}^{n+\frac{1}{2}}\right) \] \[ \Delta t \le \frac{\Delta x}{c} \]
Minimal example
Gaussian pulse incident on a dielectric slab. Extract \(R(\omega)\) and \(T(\omega)\) from monitor signals.

5. 1D finite element method (FEM)

Method idea
The differential equation is converted into a weak form and assembled into a sparse linear system.
Governing form
\[ \frac{d}{dx}\left(\frac{1}{\mu}\frac{dE}{dx}\right) + \omega^2 \varepsilon(x) E = f(x) \] \[ \int_{\Omega} \frac{d v}{d x}\frac{1}{\mu}\frac{dE}{d x}\,dx - \omega^2 \int_{\Omega} v\,\varepsilon(x)\,E\,dx = \int_{\Omega} v\,f\,dx \]
Minimal example
Layered dielectric interval solved at a selected frequency. Useful for understanding weak form, assembly, and boundary conditions.

6. 2D finite-difference time-domain (FDTD)

Method idea
In 2D, Maxwell's equations split into TM and TE modes. The grid is staggered in both \(x\) and \(y\) directions (Yee Cell).
Field Set
For **TEz** mode: primary components are \(E_x\), \(E_y\), and \(H_z\).
Use case
Scattering from pillars, waveguides, and 2D photonic crystals.
\[ E_x^{n+1}(i+\frac{1}{2}, j) = E_x^n(i+\frac{1}{2}, j) + \frac{\Delta t}{\varepsilon \Delta y} \left[ H_z^{n+\frac{1}{2}}(i+\frac{1}{2}, j+\frac{1}{2}) - H_z^{n+\frac{1}{2}}(i+\frac{1}{2}, j-\frac{1}{2}) \right] \] \[ E_y^{n+1}(i, j+\frac{1}{2}) = E_y^n(i, j+\frac{1}{2}) - \frac{\Delta t}{\varepsilon \Delta x} \left[ H_z^{n+\frac{1}{2}}(i+\frac{1}{2}, j+\frac{1}{2}) - H_z^{n+\frac{1}{2}}(i-\frac{1}{2}, j+\frac{1}{2}) \right] \] \[ H_z^{n+\frac{1}{2}}(i+\frac{1}{2}, j+\frac{1}{2}) = H_z^{n-\frac{1}{2}}(i+\frac{1}{2}, j+\frac{1}{2}) + \frac{\Delta t}{\mu \Delta y} \left[ E_x^n(i+\frac{1}{2}, j+1) - E_x^n(i+\frac{1}{2}, j) \right] - \frac{\Delta t}{\mu \Delta x} \left[ E_y^n(i+1, j+\frac{1}{2}) - E_y^n(i, j+\frac{1}{2}) \right] \]
2D Yee Cell Staggered Grid
Diagram showing the staggered placement of E and H field components in a 2D FDTD grid
Spatial discretization in 2D FDTD. Field components are offset by half a grid cell, enabling second-order accurate central difference updates.

7. 2D finite element method (FEM)

Method idea
Triangular elements make geometry part of the discretization. This is the main reason FEM is attractive for irregular domains.
Helmholtz form
\[ \nabla \cdot (\nabla u) + k_0^2\varepsilon_r(x,y)u = f(x,y) \] \[ \int_{\Omega} \nabla v \cdot \nabla u\, d\Omega - k_0^2 \int_{\Omega} \varepsilon_r(x,y)vu\, d\Omega = \int_{\Omega} vf\, d\Omega \]
Minimal example
2D dielectric inclusion in a bounded domain. Useful for mesh generation, assembly, and geometry-dependent error control.

8. Workflow and validation

Before writing code or running a model, it is crucial to understand the overarching simulation pipeline. A rigorous workflow is the only way to prevent "garbage in, garbage out" and ensure your results represent reality.

The Standard Simulation Pipeline

1

Define Geometry & Materials

Set up the physical structure and critical length scales. Assign correct material properties (\(\varepsilon, \mu, \sigma\)) and verify unit conversions.

2

Select Discretization (Mesh)

Choose the appropriate solver (FDTD/FEM) and generate a grid or mesh that adequately resolves all structural features and wavelengths.

3

Apply Boundaries & Sources

Define how energy enters the system (ports, plane waves) and how the simulation edges behave. (This is heavily detailed in Section 9).

4

Solve & Extract Observables

Run the simulation and extract engineering observables (S-parameters, Reflectance) rather than relying solely on visual field plots.

5

Validation & Convergence

Ensure the extracted observables are independent of the mesh density and boundary distances.

Critical Checks

  • Mesh refinement sensitivity.
  • Time-step (CFL) limits for FDTD.
  • Domain-size and boundary-distance sensitivity.
  • Comparison against analytical reference models.

Common Failure Modes

  • Under-resolved interfaces or curved geometry.
  • Boundary reflection contaminating the answer.
  • Wrong normalization of \(R\), \(T\), or S-parameters.
  • Incorrect material data or unit scales.
\[ \delta_Q^{(k)} = \frac{|Q_{k+1}-Q_k|}{|Q_{k+1}| + 10^{-16}} \]

Use \(Q\) for any meaningful observable: reflected power, resonant frequency, transmission peak, or field energy. Convergence is defined on the observable, not on visual comfort.

As shown in the convergence curve, the observable value approaches a stable "converged value" as the mesh density increases. A proper mesh convergence check ensures that your computational results are grid-independent and physically defensible rather than an artifact of the discretization.

Mesh convergence check
Mesh convergence check graph showing observable value stabilizing as mesh density increases
Plotting an observable against mesh density. The simulation is considered converged when further mesh refinement produces negligible changes in the engineering observable of interest.
Transitioning to Boundaries: Now that the pipeline is clear, it's vital to recognize that Step 3—Sources and Boundaries—is often where simulations fail. The mathematical edges of your simulation must correctly mimic the physical setup. Section 9 explores exactly how to do this.

9. Sources and boundaries

Boundaries define the edges of your mathematical model and drastically influence accuracy. They are not a cosmetic implementation detail. The choice of boundary depends entirely on the physics of your problem.

Open Space / PML (Perfectly Matched Layer)

What it is: PML is an artificial absorbing sponge placed at the edges of the simulation domain to absorb outgoing waves without causing reflections.

When to use: Primarily used for modeling infinite free space, such as in scattering problems (e.g., radar cross-section, Mie scattering) or antenna propagation analysis.

How to use: PMLs are often paired with a Total-Field Scattered-Field (TFSF) boundary. A plane wave is injected at the connecting surface. The innermost area (Region 1) calculates the total field (incident + scattered), while the outer area (Region 2) tracks only the scattered field before it is safely absorbed by the PML (Region 3).

TF/SF source injection with PML
Simulation domain showing Total Field, Scattered Field, and PML regions
The interacting structure is enclosed by a connecting surface where the source is injected. The outer PML layer ensures no unphysical reflections return to the Scattered Field region.

Periodic Boundary Conditions (PBC)

What it is: PBCs wrap the fields from one side of the simulation domain to the opposite side, enforcing a periodic phase relationship.

When to use: Indispensable for modeling infinitely repeating structures such as metamaterials, frequency selective surfaces, and optical periodic gratings.

How to use: The simulation domain should cover exactly one unit cell of your structure. Apply a soft-source excitation to inject the wave. The transverse boundaries are set to PBC to maintain the infinite lateral extent, while the propagation axis is terminated with PML to absorb transmission and reflection.

Periodic Boundary Condition
Simulation domain showing Periodic Boundary Conditions
A single unit cell representing an infinite array. The incident wave propagates longitudinally, while PBCs explicitly map the fields horizontally.

Symmetry Boundaries (PEC / PMC)

Symmetry boundary conditions are optimization tools that can significantly reduce the required simulation volume and computation time by factors of 2, 4, or 8.

What it is: Perfect Electric Conductor (Symmetric) or Perfect Magnetic Conductor (Anti-Symmetric) boundaries force specific field components to zero along the plane of symmetry.

When to use: You can only use these when both the electromagnetic fields and the physical structure have a true plane of symmetry directly through the middle of the simulation region.

How to use: The correct option depends on your source polarization:

  • If the source polarization is tangential to the plane of symmetry, select the symmetry option (PEC / Symmetric).
  • If the source polarization is normal to the plane of symmetry, select the anti-symmetry option (PMC / Anti-Symmetric).
Important Limitations and Rules:
  • No warnings for wrong settings: It is easy to accidentally select the wrong symmetry option. The software will not display a warning, but the simulation results will be completely incorrect. Always verify by running a full non-symmetric model first.
  • Do not manually shrink the domain: Keep the original dimensions of the simulation region. If selected correctly, the solver will automatically shade the portion that isn't directly simulated and will automatically unfold the field data during analysis.

10. Example codes

1D FDTD slab

Gaussian pulse, dielectric slab, FFT-based \(R/T\) extraction.

Download fdtd_1d_slab.py

1D FEM layered medium

Linear elements, piecewise \(\varepsilon\), sparse solve.

Download fem_1d_layered.py

2D FDTD without ABC

Shows boundary reflection contamination in an open-region transient problem.

Download fdtd_2d_no_abc.py

2D FDTD with Mur ABC

Simple absorbing boundary demonstration.

Download fdtd_2d_mur_abc.py

2D FEM Helmholtz

Triangular mesh assembly for a dielectric inclusion.

Download fem_2d_helmholtz.py

Suggested self-study path

  • Reproduce \(R/T\) for a dielectric slab
  • Add loss and compare \(A\)
  • Move from ABC to CPML
  • Compare FDTD and FEM on one shared structure

Selected engineering examples and validation notes

Finite-Difference Time-Domain (FDTD) Method

What this example demonstrates

Below I provide a two-dimensional TE FDTD implementation that includes:

  • Yee-grid time stepping for TE polarization in 2D
  • Total-field / scattered-field (TF/SF) formulation to inject an incident wave while separating scattered fields
  • PML absorbing boundaries to approximate an open domain

The demonstration problem is electromagnetic scattering from a core–shell nanocylinder. The relative permittivities of the shell and core are set to 4 and 1, respectively. The simulation domain is divided into a total-field region, a scattered-field region, and surrounding PML layers. The animation visualizes the time evolution of the field (Ey) and provides an intuitive view of how the incident wave interacts with the core–shell structure and generates scattered waves.

The MATLAB implementation can be accessed here.

Notes on interpretation
  • Why time-domain? In frequency-domain solvers you typically solve one frequency per run; in FDTD, a broadband pulse can reveal a wide spectral response from a single transient simulation.
  • Why TF/SF? TF/SF allows clean separation between incident and scattered fields, which is particularly useful for scattering and extinction calculations.
  • Why PML matters? Boundary reflections can contaminate scattering signatures; a properly configured PML reduces this artifact and improves the reliability of near- and far-field observables.
  1. Kane Yee (1966). “Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media”. IEEE Transactions on Antennas and Propagation, 14(3), 302–307.
  2. Allen Taflove and Susan C. Hagness (2005). Computational Electrodynamics: The Finite-Difference Time-Domain Method, 3rd ed. Artech House Publishers. ISBN 978-1-58053-832-9.

Plasmonic Metal Nanoparticles

Over the last two decades, plasmon resonances in metallic nanoparticles have been studied extensively because metallic nanostructures can strongly absorb and scatter light through localized surface plasmon resonances. This behavior underpins a broad range of applications, including optical sensing, photothermal heating and solar vapor generation, and bio-imaging.

Here I provide a 2D-TE FDTD code to compute extinction, scattering, and absorption efficiencies for a gold-shell dielectric-core nanocylinder. The Fortran implementation can be accessed here.

Analytical benchmark for verification

For 2D scattering problems with cylindrical symmetry, it is valuable to have an analytical solution as a reference for verifying numerical results. For this reason, I also provide an analytical code to calculate the extinction spectrum of a core–shell nanocylinder. As an example, results for a gold-shell dielectric-core nanocylinder with outer radius 50 nm and inner radius 40 nm are compared against the corresponding FDTD simulations in the figure below. The analytical implementation can be accessed here.

Comparison between FDTD and analytical results for a core-shell nanocylinder
Comparison between FDTD and analytical results for a core–shell nanocylinder.
  • Why this comparison is useful: matching peak positions and spectral trends provides a practical sanity check on discretization, TF/SF injection, and boundary absorption settings.
  • What can differ: small deviations may arise from finite grid resolution, numerical dispersion, PML configuration, and the time-window / Fourier-transform settings used to extract spectra from transient fields.

Back to Teaching