Computational Electromagnetics: Finite Difference vs. Finite Element Method

Closed-form analytical solutions to Maxwell’s equations typically exist for highly idealized geometries and canonical boundary conditions. For realistic electromagnetic structures—inhomogeneous media, curved boundaries, or complex scatterers—numerical PDE solvers are essential.

Every deterministic PDE solver in computational electromagnetics follows the same three-stage pipeline:

Governing Equation: Define the continuous differential operator and associated boundary conditions.
Discretization Scheme: Project the continuous operator onto a discrete computational grid or basis space to construct a sparse linear system $\mathbf{A}\mathbf{x} = \mathbf{b}$.
Linear Solve: Solve the matrix equation via direct factorization (e.g., sparse LU) or iterative Krylov methods (e.g., BiCGSTAB, GMRES).

While time-domain approaches like FDTD evolve Maxwell’s equations step-by-step in time, frequency-domain boundary value problems convert steady-state wave behavior directly into linear systems. Below, we compare how the Finite Difference (FD) and Finite Element Method (FEM) approach the same 1D scalar Helmholtz equation, culminating in an asymptotic comparison of their resulting system matrices.

1.1 The Continuous Model: 1D Scalar Helmholtz Equation

Consider time-harmonic fields with an implied $e^{j\omega t}$ time dependence. For a transverse electric (TE) mode propagating along $x$ with fields polarized as $\mathbf{E} = E_y(x)\hat{y}$ and $\mathbf{H} = H_z(x)\hat{z}$, Maxwell’s curl equations reduce to:

$\frac{dE_y}{dx} = -j\omega\mu_0 H_z \quad (1.1)$ $\frac{dH_z}{dx} = -j\omega\epsilon_0 \epsilon_r(x) E_y \quad (1.2)$

Differentiating (1.1) with respect to $x$ and substituting into (1.2) yields the 1D scalar Helmholtz equation on the domain $\Omega = (a, b)$:

\[\frac{d^2 E_y}{dx^2} + k_0^2 \epsilon_r(x) E_y = 0 \quad (1.3)\]

where $k_0 = \omega\sqrt{\mu_0 \epsilon_0}$ is the free-space wavenumber and $\epsilon_r(x)$ is the relative permittivity profile.

Boundary Conditions: Dirichlet vs. Neumann

To ensure equation (1.3) is well-posed and the resulting linear system is full rank, boundary conditions must be enforced at $x = a$ and $x = b$:

Essential (Dirichlet) Boundary Conditions: Prescribe the primary field value directly: $E_y(a) = K_a, \quad E_y(b) = K_b \quad (1.4a)$ In electromagnetics, this represents an ideal electric conductor (PEC) when set to zero, or an applied voltage/source constraint when nonzero.
Natural (Neumann) Boundary Conditions: Prescribe the normal derivative of the field: $\frac{dE_y}{dx}\Big|_{x=a} = C_a, \quad \frac{dE_y}{dx}\Big|_{x=b} = C_b \quad (1.4b)$ From (1.1), $\frac{dE_y}{dx} \propto H_z$. Specifying a Neumann condition directly constrains the tangential magnetic field (e.g., modeling a perfect magnetic conductor, PMC, or an injected surface current).

1.2 The Finite Difference Approach (Strong Form)

The Finite Difference (FD) method works directly on the strong form of the PDE (1.3). It discretizes the continuous domain into $N$ uniformly spaced grid nodes $x_m = a + (m-1)\Delta$ with pitch $\Delta = (b - a)/(N - 1)$.

Continuous differential operators are replaced by local difference quotients derived via Taylor series expansions:

\[E(x_m \pm \Delta) = E(x_m) \pm \Delta \frac{dE}{dx}\Big|_{x_m} + \frac{\Delta^2}{2}\frac{d^2 E}{dx^2}\Big|_{x_m} \pm \frac{\Delta^3}{6}\frac{d^3 E}{dx^3}\Big|_{x_m} + \mathcal{O}(\Delta^4)\]

Summing the forward and backward expansions cancels the odd-order terms, yielding the canonical second-order central difference stencil:

\[\frac{d^2 E_y}{dx^2}\Big|_{x_m} = \frac{E_{m+1} - 2E_m + E_{m-1}}{\Delta^2} + \mathcal{O}(\Delta^2) \quad (1.5)\]

Substituting stencil (1.5) into the Helmholtz equation (1.3) at interior node $m$ gives the algebraic balance:

\[\frac{1}{\Delta^2} E_{m-1} + \left(-\frac{2}{\Delta^2} + k_0^2 \epsilon_{r,m}\right) E_m + \frac{1}{\Delta^2} E_{m+1} = 0 \quad (1.6)\]

Assembling the Linear System

For an $N$-node domain with Dirichlet boundary conditions $E_1 = K_a$ and $E_N = K_b$, the interior stencils along with the boundary equations yield a sparse tridiagonal matrix equation $\mathbf{A}_{\text{FD}} \mathbf{E} = \mathbf{b}$:

\[\begin{bmatrix} 1 & 0 & 0 & \dots & 0 \\ \frac{1}{\Delta^2} & \left(-\frac{2}{\Delta^2} + k_0^2 \epsilon_2\right) & \frac{1}{\Delta^2} & \dots & 0 \\ 0 & \frac{1}{\Delta^2} & \left(-\frac{2}{\Delta^2} + k_0^2 \epsilon_3\right) & \frac{1}{\Delta^2} & \dots \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & \dots & \frac{1}{\Delta^2} & \left(-\frac{2}{\Delta^2} + k_0^2 \epsilon_{N-1}\right) & \frac{1}{\Delta^2} \\ 0 & \dots & \dots & 0 & 1 \end{bmatrix} \begin{bmatrix} E_1 \\ E_2 \\ E_3 \\ \vdots \\ E_{N-1} \\ E_N \end{bmatrix} = \begin{bmatrix} K_a \\ 0 \\ 0 \\ \vdots \\ 0 \\ K_b \end{bmatrix}\]

Note on Boundary Enforcement in FD:

Dirichlet: Implemented explicitly by replacing the boundary rows with identity stencils ($A_{1,1}=1, b_1=K_a$).

Neumann: Because FD evaluates derivatives pointwise, enforcing $\frac{dE}{dx}\big _{x_1} = C_a$ requires either a one-sided finite difference or an imaginary “ghost node” $x_0$ located at $a - \Delta$. Eliminating $E_0$ through central difference substitution couples the boundary equation into the second-derivative stencil.

1.3 The Finite Element Approach (Weak Form)

Rather than forcing the PDE to evaluate to zero point-by-point, the Finite Element Method (FEM) satisfies the equation in an integral. This requires reformulating (1.3) into its weak (variational) form.

Weak Formulation and Variational Principles

Multiply (1.3) by an arbitrary smooth test function $T(x)$ and integrate over domain $\Omega$:

\[\int_a^b T(x) \left( \frac{d^2 E_y}{dx^2} + k_0^2 \epsilon_r(x) E_y \right) dx = 0 \quad (1.7)\]

Applying integration by parts to the second-order derivative shifts one differential operator onto the test function:

\[\int_a^b T \frac{d^2 E_y}{dx^2}\,dx = \left[ T \frac{dE_y}{dx} \right]_a^b - \int_a^b \frac{dT}{dx}\frac{dE_y}{dx}\,dx \quad (1.8)\]

Substituting (1.8) into (1.7) yields the continuous weak form:

\[\int_a^b \left( \frac{dT}{dx}\frac{dE_y}{dx} - k_0^2 \epsilon_r(x) T E_y \right) dx = \left. T \frac{dE_y}{dx} \right|_b - \left. T \frac{dE_y}{dx} \right|_a \quad (1.9)\]

This transformation produces two key structural advantages:

Reduced Regularity: The solution $E_y(x)$ no longer needs $\mathcal{C}^2$ continuity (twice differentiable). It only requires first derivatives that are square-integrable ($E_y \in H^1(\Omega)$), permitting piecewise polynomial approximations that naturally handle material interfaces.
Natural Boundary Handling: The boundary terms on the right-hand side of (1.9) contain the normal derivatives $\frac{dE_y}{dx}$. Neumann boundary conditions enter the weak form directly through these surface integrals without requiring ghost cells.

Galerkin Discretization and Matrix Assembly

We approximate the field by expanding it into a set of $N$ localized basis functions $B_n(x)$:

\[E_y(x) \approx \sum_{n=1}^N e_n B_n(x) \quad (1.10)\]

For a 1D mesh with element size $\Delta$, the canonical choice is the piecewise-linear hat function:

\[B_n(x) = \begin{cases} \frac{x - x_{n-1}}{\Delta}, & x_{n-1} \le x \le x_n \\ \frac{x_{n+1} - x}{\Delta}, & x_n \le x \le x_{n+1} \\ 0, & \text{otherwise} \end{cases} \quad (1.11)\]

Using the Galerkin method, we choose test functions identical to our trial bases: $T_m(x) = B_m(x)$. Assuming Dirichlet conditions at the boundaries, the test functions vanish at $x=a$ and $x=b$, dropping the boundary terms in (1.9). Substituting (1.10) into (1.9) yields the system:

\[\sum_{n=1}^N e_n \int_a^b \left( \frac{dB_m}{dx} \frac{dB_n}{dx} - k_0^2 \epsilon_r(x) B_m B_n \right) dx = 0 \quad (1.12)\]

Assuming piecewise-constant permittivity within each element ($\epsilon_r(x) = \epsilon_m$ for $x \in [x_m, x_{m+1}]$), the integrals over overlapping basis triangles evaluate element-by-element to:

Stiffness matrix term: $\int \frac{dB_m}{dx}\frac{dB_n}{dx}\,dx$
Mass matrix term: $\int B_m B_n\,dx$

For row $m$, the non-zero couplings evaluate to:

$Z_{m,m-1} = -\frac{1}{\Delta} - \frac{k_0^2 \Delta}{6}\epsilon_{m-1} \quad (1.13)$ $Z_{m,m} = \frac{2}{\Delta} - \frac{k_0^2 \Delta}{3}\left(\frac{\epsilon_{m-1} + \epsilon_m}{2}\right) \quad (1.14)$ $Z_{m,m+1} = -\frac{1}{\Delta} - \frac{k_0^2 \Delta}{6}\epsilon_m \quad (1.15)$

Dividing through by $-\Delta$ scales the system to match the dimensions of the finite difference operator:

\[\mathbf{A}_{\text{FEM}} = \begin{bmatrix} 1 & 0 & 0 & \dots & 0 \\ \left(\frac{1}{\Delta^2} + \frac{k_0^2 \epsilon_1}{6}\right) & \left(-\frac{2}{\Delta^2} + \frac{k_0^2(\epsilon_1 + \epsilon_2)}{3}\right) & \left(\frac{1}{\Delta^2} + \frac{k_0^2 \epsilon_2}{6}\right) & \dots & 0 \\ 0 & \left(\frac{1}{\Delta^2} + \frac{k_0^2 \epsilon_2}{6}\right) & \left(-\frac{2}{\Delta^2} + \frac{k_0^2(\epsilon_2 + \epsilon_3)}{3}\right) & \dots & 0 \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & \dots & \dots & 0 & 1 \end{bmatrix}\]

1.4 Method Comparison: Asymptotic Equivalence and Structural Differences

Comparing $\mathbf{A}{\text{FD}}$ directly with $\mathbf{A}{\text{FEM}}$ highlights the fundamental bridge between the two numerical philosophies:

The Stiffness Operator: The spatial derivative term $\frac{1}{\Delta^2}[1, -2, 1]$ is identical in both systems.
The Mass Operator: In FD, the permittivity term $k_0^2 \epsilon_{r,m}$ resides entirely on the diagonal (a “lumped” mass representation). In FEM, the basis overlap spreads the permittivity across adjacent nodes through the off-diagonal terms $\frac{k_0^2 \epsilon}{6}$ (a “consistent” mass matrix).
Asymptotic Convergence: In the continuum limit as the grid is refined ($\Delta \to 0$): $\frac{1}{\Delta^2} \gg \frac{k_0^2 \epsilon}{6}$ The mass matrix perturbations become negligible relative to the singular differential terms. Both schemes converge to the same leading-order continuous operator with $\mathcal{O}(\Delta^2)$ truncation error.

Choosing Between FD and FEM

Metric / Feature	Finite Difference (FD / FDFD)	Finite Element Method (FEM)
Primary Formulation	Strong form (pointwise differential quotients)	Weak form (integral variational statements)
Mesh Geometry	Structured Cartesian grids (staircasing on curves)	Unstructured conformal meshes (tetrahedra/triangles)
Material Interfaces	Requires effective index averaging at boundaries	Naturally supported across conformal element faces
Implementation Complexity	Minimal; intuitive stencil arrays	High; requires numerical quadrature, Jacobians, and assembly
Matrix Sparsity	Highly regular, banded, uniform stencils	Sparse, variable bandwidth depending on mesh topology