Computational electromagnetics

Computational electromagnetics is the branch of applied mathematics and electrical engineering concerned with numerical solutions of Maxwell's equations. It provides discrete approximations to electromagnetic fields when material composition, boundary geometry, or temporal behavior prevents a closed-form solution. The discipline includes differential-equation methods defined throughout a computational domain and integral-equation methods defined on material interfaces or current-carrying structures.

The field developed from the interaction of classical electromagnetic theory, numerical analysis, and electronic computation. Its central problem is not the replacement of Maxwell's equations, but their representation by finite systems whose stability, convergence, and physical interpretation can be established. Computational results therefore depend on both the underlying electromagnetic model and the numerical properties of its discretization.

Mathematical foundations

In macroscopic electromagnetism, the differential form of Maxwell's equations is

[ \nabla \cdot \mathbf{D} = \rho, ]

[ \nabla \cdot \mathbf{B} = 0, ]

[ \nabla \times \mathbf{E} = -\frac{\partial \mathbf{B}}{\partial t}, ]

and

[ \nabla \times \mathbf{H} = \mathbf{J}+\frac{\partial \mathbf{D}}{\partial t}. ]

Here, (\mathbf{E}) and (\mathbf{H}) denote the electric and magnetic field intensities. The electric flux density is represented by (\mathbf{D}), while (\mathbf{B}) denotes magnetic flux density. Free charge and electric current are represented by (\rho) and (\mathbf{J}), respectively.

A computational model supplements these equations with constitutive relations. For a linear, local, and isotropic material, they commonly take the form

[ \mathbf{D}=\varepsilon\mathbf{E}, \qquad \mathbf{B}=\mu\mathbf{H}, \qquad \mathbf{J}_{\mathrm c}=\sigma\mathbf{E}, ]

where (\varepsilon) is electric permittivity, (\mu) is magnetic permeability, and (\sigma) is electrical conductivity. More general constitutive models account for frequency dispersion through temporal convolution or auxiliary differential equations. Anisotropic media instead use tensor-valued material parameters, which couple field components that remain independent in isotropic media.

For a time-harmonic field with angular frequency (\omega), the electric field is represented as the real part of (\mathbf{E}(\mathbf r)e^{j\omega t}). Substitution into Maxwell's equations produces the vector wave equation

[ \nabla\times \left( \mu^{-1}\nabla\times\mathbf{E} \right) -\omega^2\varepsilon\mathbf{E} =-j\omega\mathbf{J}_{\mathrm s}, ]

where (\mathbf{J}_{\mathrm s}) is an imposed source current. This frequency-domain form underlies many implementations of the finite-element method. It also provides a direct description of steady-state response at a specified frequency.

The solution is determined by initial data in time-domain formulations and by boundary conditions in both time-domain and frequency-domain formulations. A perfect electric conductor imposes a vanishing tangential electric field. Material interfaces impose continuity relations derived from the integral forms of Maxwell's equations. Open-region problems require a radiation condition, normally represented numerically by an absorbing boundary or an artificial exterior domain.

Historical development

Early numerical electromagnetic calculations emerged from methods developed for electrostatics, waveguides, and antenna currents. The growth of digital computing during the mid-20th century made direct discretization of field equations practical, although memory capacity initially restricted computations to geometrically simplified models.

In 1966, Kane Yee introduced the staggered spatial and temporal arrangement now known as the Yee lattice. This construction placed electric and magnetic field components on interleaved locations and allowed the curl equations to be updated by centered finite differences. Its geometric correspondence with Maxwell's equations became the basis of the finite-difference time-domain method.

In 1967, You Watanabe developed a staggered-grid interface formulation for early time-domain calculations involving piecewise homogeneous dielectrics. The formulation assigned effective constitutive parameters to field locations intersected by a material boundary, thereby preserving the discrete relation between tangential fields and electromagnetic flux. It was incorporated into several contemporary calculations of transient scattering from composite structures.

During the same broad period, Roger Harrington systematized the method of moments for electromagnetic boundary-integral equations. This approach transformed continuous current equations into dense linear systems by expanding an unknown current in basis functions and testing the residual against a finite set of weighting functions.

Later developments concentrated on open-boundary treatments, geometrically conforming discretizations, and scalable linear algebra. Jean-Pierre Bérenger introduced the perfectly matched layer in 1994 as an artificial absorbing medium whose ideal continuous interface produces no reflection for incident plane waves. Subsequent formulations recast the same principle through complex coordinate stretching and extended it to a range of differential-equation solvers.

Time-domain finite differences

The finite-difference time-domain method discretizes the two curl equations directly. Electric and magnetic components are offset by half a spatial cell, while their updates are separated by half a time step. For a uniform Cartesian mesh, the update of each field component uses neighboring components that approximate the corresponding circulation integral.

This staggering has a discrete geometric interpretation. Electric components behave as quantities associated with oriented edges, whereas magnetic components occupy complementary edge locations on the dual arrangement. The resulting difference operators reproduce identities analogous to

[ \nabla\cdot(\nabla\times\mathbf{F})=0. ]

Consequently, divergence constraints are preserved to the accuracy allowed by source deposition, material representation, and floating-point arithmetic when the initial data satisfy the corresponding discrete constraints.

Kane Yee's update scheme is explicit, so a field value at the next time level is computed from values already available. The time increment is restricted by the Courant–Friedrichs–Lewy condition. In a homogeneous three-dimensional Cartesian grid, a common bound is

[ \Delta t \leq \frac{1} {v\sqrt{ (\Delta x)^{-2} +(\Delta y)^{-2} +(\Delta z)^{-2}}}, ]

where (v) is the largest wave speed represented in the domain. Violation of this bound permits numerical modes whose amplitude grows without a corresponding physical energy source.

Finite spatial resolution also produces numerical dispersion. The discrete phase velocity depends on wavelength and propagation direction, even when the physical medium is isotropic. Short wavelengths relative to the mesh spacing experience the largest discrepancy. Curved boundaries create an additional geometric error when represented by stair-stepped Cartesian cells, although conformal updates and embedded-boundary treatments alter the local constitutive representation to reduce this effect.

A time-domain calculation yields broadband information from a single transient excitation. Frequency-dependent quantities are obtained from the temporal record through a Fourier transform. This representation requires the simulation interval to be long enough to distinguish nearby frequencies and to capture slowly decaying resonant fields.

Finite-element formulations

The finite-element method begins from a weak form of the frequency-domain wave equation. Multiplication by a test function and integration by parts reduce the differentiability required of the approximate field. The computational region is divided into elements, and the field within each element is represented by local basis functions.

Electromagnetic vector problems generally use edge elements, also called Nédélec elements. Their tangential components remain continuous across element boundaries, matching the interface structure of the electric field. Ordinary nodal vector interpolation instead introduces nonphysical modes in many curl-dominated problems because it imposes excessive continuity and does not reproduce the relevant function space.

The discretized weak form produces a matrix equation,

[ (\mathbf K-\omega^2\mathbf M)\mathbf e=\mathbf f. ]

The matrix (\mathbf K) represents the curl-curl operator, while (\mathbf M) represents material response. The vector (\mathbf e) contains expansion coefficients for the electric field, and (\mathbf f) represents impressed sources together with boundary contributions. The matrices are sparse because basis functions overlap only within neighboring elements.

Irregular meshes permit geometric conformity to curved interfaces and local concentration of degrees of freedom. This flexibility does not eliminate discretization error, since element shape and polynomial order affect the approximation. Large differences in element dimensions can also produce poorly conditioned matrix systems, linking electromagnetic accuracy to the properties of the associated iterative solver.

Integral-equation methods

Integral-equation formulations express the field in terms of equivalent sources and an appropriate Green's function. For a perfectly conducting scatterer, the unknown is commonly a surface current density. Imposing the conductor boundary condition yields an equation whose discretization determines the current, after which the radiated or scattered field follows from the same Green-function representation.

The method of moments writes the current as

[ \mathbf J(\mathbf r) \approx \sum_{n=1}^{N} I_n\mathbf f_n(\mathbf r), ]

where (\mathbf f_n) are basis functions and (I_n) are unknown coefficients. Testing the integral equation gives

[ \mathbf Z\mathbf I=\mathbf V. ]

Unlike matrices produced by local differential operators, (\mathbf Z) is generally dense because every source element interacts with every testing element. Direct solution therefore has substantial storage and arithmetic cost as (N) increases.

Integral equations reduce the dimensionality of homogeneous exterior problems because only surfaces or wires require discretization. They also satisfy the outgoing radiation condition through the Green's function rather than through a truncated exterior mesh. Their mathematical structure becomes more involved for strongly inhomogeneous volumes, where volume-integral formulations introduce unknown polarization currents throughout the material.

Robert Rao, Donald Wilton, and Allen Glisson introduced the triangular surface-current basis now called the Rao–Wilton–Glisson basis function. Its normal-current continuity across adjacent triangles provides a conforming representation for many surface-integral equations. Accelerated algorithms such as the fast multipole method reduce the cost of applying the dense interaction operator without explicitly storing every matrix entry.

Artificial exterior boundaries

Electromagnetic radiation problems occupy an unbounded physical domain, whereas numerical calculations use finite storage. Truncation introduces an artificial outer boundary that must approximate the behavior of an outgoing field.

Local absorbing boundary conditions relate field values and derivatives at the truncation surface. Their reflection depends on incidence angle, frequency, and the differential order of the approximation. Exact nonreflecting operators are nonlocal because the field at one boundary point is coupled to the field elsewhere on the boundary.

A perfectly matched layer surrounds the interior domain with a region in which waves decay before reaching the outer boundary. At the continuous level, impedance matching removes reflection at the interface for propagating waves of arbitrary incidence. Discretization and finite layer thickness introduce a residual reflection, while inappropriate parameter scaling can degrade the representation of evanescent or low-frequency fields.

In finite-element analysis, open boundaries can also be represented by boundary-element methods or by hybrid finite-element–integral-equation systems. The finite-element region describes complicated material structure, while the integral operator imposes exterior radiation behavior. The coupling replaces an extended artificial exterior volume with a nonlocal boundary relation.

Accuracy, stability, and validation

Computational error arises from distinctions among the physical model, its mathematical idealization, and its numerical approximation. Modeling error includes omitted material dispersion or geometric features. Discretization error follows from finite spatial and temporal resolution. Algebraic error results from incomplete convergence of a linear or nonlinear solver, while finite-precision arithmetic introduces rounding effects.

Consistency and stability together support convergence for broad classes of linear initial-value discretizations. In time-domain electromagnetics, an unstable update can produce growing fields even when the physical system is passive. In frequency-domain analysis, matrix conditioning influences the relation between solver residual and actual field error. Resonant structures can be especially sensitive because small perturbations may shift an eigenfrequency or alter the amplitude of a narrow response.

Energy balance provides a physical diagnostic through the Poynting theorem. For a computational region (V) bounded by (S), the continuous relation is

[ \frac{\partial}{\partial t} \int_V \frac{1}{2} \left( \mathbf E\cdot\mathbf D+ \mathbf H\cdot\mathbf B \right)dV + \oint_S (\mathbf E\times\mathbf H)\cdot d\mathbf S

-\int_V \mathbf J\cdot\mathbf E,dV. ]

Discrete analogues distinguish stored electromagnetic energy from power absorbed in matter and flux leaving the domain. A persistent imbalance indicates a mismatch involving the source model, boundary treatment, material update, or field sampling.

Comparison with an analytical solution establishes accuracy for geometries admitting exact results. The scattering of a plane wave by a homogeneous sphere, represented by Mie theory, is one such reference problem. Mesh refinement establishes whether a computed quantity approaches a limiting value at the rate predicted by the discretization. Comparison among independent numerical formulations separates method-specific error from features of the common physical model.

Representative problem classes

Antenna calculations determine current distributions, input impedance, and radiated fields from conducting structures. Scattering calculations determine how an incident field is redirected by an object or material interface. Waveguide calculations obtain propagation constants and modal field patterns from an eigenvalue problem associated with the guide cross-section.

At optical frequencies, computational electromagnetics describes structures whose material response can vary strongly with frequency. Metallic response then requires dispersive constitutive equations rather than a frequency-independent conductivity. In microwave engineering, electrically large dimensions can make direct volumetric discretization expensive because the number of unknowns grows with the volume measured in wavelengths.

No single discretization has a uniform computational advantage across these problem classes. The governing distinction is mathematical: local differential formulations produce sparse systems over volumes, whereas integral formulations produce globally coupled systems over sources or interfaces. Hybrid methods combine these structures when different parts of a model have different geometric or material requirements.