Yee FDTD and transient Maxwell finite elements

Let $F=\{(x_0,y,z):y_0\le y\le y_1,z_0\le z\le z_1\}$ be an oriented rectangle with normal $e_x$, and let $E$ be a $C^1$ vector field near $F$. Orient $\mathrm{\partial }F$ by the outward-normal-first convention. Then ${\int }_{\mathrm{\partial }F}E\cdot tds={\int }_F(\mathrm{curl}E)\cdot e_{x}dS={\int }_{y_0}^{y_1}{\int }_{z_0}^{z_1}({\mathrm{\partial }}_y{E}_z-{\mathrm{\partial }}_z{E}_y)(x_0,y,z)dzdy\mathrm{.}$ The identities for normals $e_y$ and $e_z$ follow by cyclic permutation of coordinates. Reversing a face orientation reverses both integrals.

Let $\mathrm{\Omega }\subset {\mathbb{R}}^3$ be a bounded Lipschitz polyhedron, or a periodic rectangular cell. The real electric field is $E(t,x)\in {\mathbb{R}}^3$ and the magnetic flux density is $B(t,x)\in {\mathbb{R}}^3$; the magnetic intensity is $H={\mu }^{-1}B$. Let $ϵ(x)$ and $\mu (x)$ be measurable, time-independent, symmetric matrix fields for which there are constants $0<{ϵ}_0\le {ϵ}_1$ and $0<{\mu }_0\le {\mu }_1$ with ${ϵ}_0\mathrm{\mid }v{\mathrm{\mid }}^2\le v^Tϵ(x)v\le {ϵ}_1\mathrm{\mid }v{\mathrm{\mid }}^2,{\mu }_0\mathrm{\mid }v{\mathrm{\mid }}^2\le v^T\mu (x)v\le {\mu }_1\mathrm{\mid }v{\mathrm{\mid }}^2$ for almost every $x$ and every $v\in {\mathbb{R}}^3$. A conductivity tensor $\sigma (x)$, when present, is measurable, symmetric, and positive semidefinite. The coefficients may jump across finitely many Lipschitz material interfaces. No derivative of a coefficient is part of the model.

In each material region the constitutive fields satisfy ${\mathrm{\partial }}_t(ϵE)=\mathrm{curl}({\mu }^{-1}B)-J_{\mathrm{s}\mathrm{r}\mathrm{c}}-\sigma E,{\mathrm{\partial }}_{t}B+\mathrm{curl}E=0.$ Here $J_{\mathrm{s}\mathrm{r}\mathrm{c}}$ is the impressed free-current density, while $\sigma E$ is the conductive current. Electric charge $\rho$ is constrained by $\mathrm{div}(ϵE)=\rho$ and the total-current continuity equation ${\mathrm{\partial }}_t\rho +\mathrm{div}(J_{\mathrm{s}\mathrm{r}\mathrm{c}}+\sigma E)=0$; the magnetic constraint is $\mathrm{div}B=0$. These equalities are distributional. In the lossless branch set $\sigma =0$.

The boundary choice is one of the following. A periodic boundary identifies opposite traces of the rectangular cell. A perfect electric conductor (PEC) boundary imposes zero tangential trace $n\times E=0$ on $\mathrm{\partial }\mathrm{\Omega }$. At a material interface with no singular surface charge or current, the weak equations are understood with the transmission conditions $[n\times E]=0$, $[n\cdot B]=0$, $[n\times ({\mu }^{-1}B)]=0$, and $[n\cdot (ϵE)]=0$, where $[q]$ is the trace from one side minus the trace from the other. The interface conditions are part of the model convention; no smooth Taylor expansion across a coefficient jump is asserted.

For PEC, the weak Ampere equation is $(ϵE_t,v{)}_{L^2}-({\mu }^{-1}B,\mathrm{curl}v{)}_{L^2}=-(J_{\mathrm{s}\mathrm{r}\mathrm{c}}+\sigma E,v{)}_{L^2}$ for every $v\in H_0(\mathrm{curl};\mathrm{\Omega })$ as defined in Zero tangential trace H(curl) space; the periodic branch uses periodic $H(\mathrm{curl})$ test fields. Faraday’s equation remains $B_t+\mathrm{curl}E=0$ in distributions. These weak forms define the test-space and boundary convention used by the edge finite element method.

The weighted field norms are $\mathrm{\parallel }E{\mathrm{\parallel }}_{ϵ}^2={\int }_{\mathrm{\Omega }}{E}^TϵEdx$ and $\mathrm{\parallel }B{\mathrm{\parallel }}_{{\mu }^{-1}}^2={\int }_{\mathrm{\Omega }}{B}^T{\mu }^{-1}Bdx$. The uniform coefficient bounds make them equivalent to the ordinary $L^2$ norms. All later smooth-error claims require smoothness separately on each mesh-aligned material region.

Let $K$ be a finite oriented three-dimensional cell complex with vertex, edge, face, and volume cochains. An edge cochain records oriented electric line integrals, a face cochain records oriented magnetic fluxes, and a vertex cochain records electric charge. Write $N_v,N_e,N_f,N_c$ for the retained vertex, edge, face, and cell counts. On a periodic cell, identify opposite cells with orientation preserved and let $G,C,D$ be the full vertex-to-edge gradient, edge-to-face curl, and face-to-cell divergence incidence matrices, each the transpose of its signed oriented boundary matrix. For a PEC closure, “interior” means not in the boundary subcomplex. Keep only interior electric edges as unknowns and only interior vertices as volume-charge locations, but keep all face fluxes, including boundary-normal magnetic fluxes, and all cells. This closure has no separate boundary surface-charge cochain. Define $G$ and $C$ by restricting the absolute vertex-edge and edge-face incidence matrices to those rows and columns. This mixed closure imposes zero tangential electric edge cochains on the boundary while allowing a fixed normal magnetic flux on boundary faces. In both cases, $G\in {\mathbb{R}}^{N_e\times N_v}$, $C\in {\mathbb{R}}^{N_f\times N_e}$, and $D\in {\mathbb{R}}^{N_c\times N_f}$ satisfy $CG=0$ and $DC=0$.

The electric, magnetic, and conductive material matrices are real symmetric matrices $M_e\in {\mathbb{R}}^{N_e\times N_e}$, $M_b\in {\mathbb{R}}^{N_f\times N_f}$, and $S\in {\mathbb{R}}^{N_e\times N_e}$, respectively. The first two are positive definite and represent the discrete weighted inner products induced by $ϵ$ and ${\mu }^{-1}$; $S$ is positive semidefinite and represents conductivity. Write $⟨e_1,e_2{⟩}_{M_e}=e_1^T{M}_e{e}_2$, $⟨b_1,b_2{⟩}_{M_b}=b_1^T{M}_b{b}_2$, and $\mathrm{\parallel }e{\mathrm{\parallel }}_{M_e}^2=⟨e,e{⟩}_{M_e}$, $\mathrm{\parallel }b{\mathrm{\parallel }}_{M_b}^2=⟨b,b{⟩}_{M_b}$. In a homogeneous isotropic periodic Cartesian grid with spacings $h_x,h_y,h_z$ and cell volume $V=h_x{h}_y{h}_z$, the diagonal entries for an $x$-edge and an $x$-normal face are $ϵV\mathrm{/}{h}_x^2=ϵh_y{h}_z\mathrm{/}{h}_x$ and ${\mu }^{-1}{h}_x\mathrm{/}(h_y{h}_z)$; the $y$ and $z$ entries are obtained by cyclic permutation. These values are dual-volume Hodge weights for line-integral and face-flux cochains. A compatible PEC closure uses the corresponding truncated dual volumes for free electric edges, while the magnetic mass includes boundary faces; the resulting matrices remain positive definite. The discrete constraints are $-G^T{M}_e{e}^n={\rho }^n$ on interior vertices and $D{b}^{n-1\mathrm{/}2}=0$ on all cells.

The minus sign in the electric constraint is essential: an edge oriented from tail to head has incidence entries $-1,+1$, so $-G^T$ is the outgoing-flux divergence. Thus $\rho$ represents physical electric charge with the sign of $\mathrm{div}(ϵE)$ in Transient Maxwell system with piecewise materials and boundary conditions.

For current cochains $j^{n+1\mathrm{/}2}$ and charge rows ${\rho }^n$, the discrete continuity condition means ${\rho }^{n+1}-{\rho }^n-\mathrm{\Delta }t{G}^T(j^{n+1\mathrm{/}2}+S\frac{e^{n+1}+e^n}{2})=0.$ For PEC this is the interior-vertex volume-charge relation; this closure has no separate boundary surface-charge/current cochain.

The cochain identities $CG=0$ and $DC=0$ are not imposed as an extra numerical approximation: they follow from oriented boundary-of-boundary cancellation, including the PEC row and column restrictions, as proved in Incidence curl has zero discrete divergence and gradient curl ↗. No cell geometry picture is used as a proof.

For the incidence matrices $G,C,D$ of Oriented Maxwell cochains, masses, and the Yee grid, one has $CG=0$ and $DC=0$. The conclusion holds for the full periodic complex and for the PEC mixed closure defined there.

Let $M_e$ be the positive definite electric mass matrix, let $S$ be the positive semidefinite conductivity matrix, and let $\mathrm{\Delta }t>0$. Then $M_e+(\mathrm{\Delta }t\mathrm{/}2)S$ is symmetric positive definite, hence invertible. Thus each midpoint-conductivity Yee electric row has a unique solution.

Let $\mathrm{\Omega }$ be a bounded contractible Lipschitz polyhedron with a finite matching affine tetrahedral mesh ${\mathcal{T}}_h$, and use the real lowest-order spaces $W_h^0,N_h^1,N_h^2,N_h^3$ defined in Lowest-order Whitney and Nedelec de Rham spaces. The electric trial field is $E_h(t)\in N_h^1\cap H_0(\mathrm{curl};\mathrm{\Omega })$; the magnetic flux field is $B_h(t)\in N_h^2$, with no imposed zero normal trace on $B_h$. PEC means the closure-defined zero tangential subspace from Zero tangential trace H(curl) space. Periodic analogues use the periodic cell complex. Countable Choice is inherited from the H(div) and affine Piola space definitions in the suppliers; the finite matrix operations below use no additional choice.

Put $V_h:=N_h^1\cap H_0(\mathrm{curl};\mathrm{\Omega })$. Choose an edge basis $\{{\varphi }_i{\}}_{i=1}^{N_e}$ of the electric trial space $V_h$ and a face basis $\{{\psi }_a{\}}_{a=1}^{N_f}$ of $N_h^2$, so $N_e=\dim V_h$ and $N_f=\dim N_h^2$. Define the symmetric mass matrices $(M_e{)}_{ij}={\int }_{\mathrm{\Omega }}{\varphi }_j^Tϵ{\varphi }_{i}dx,(M_b{)}_{ab}={\int }_{\mathrm{\Omega }}{\psi }_b^T{\mu }^{-1}{\psi }_{a}dx\mathrm{.}$ They are positive definite by the coefficient bounds in Transient Maxwell system with piecewise materials and boundary conditions and linear independence of the finite element bases. Let $C$ and $D$ be the coefficient matrices, in these chosen bases, of the restricted map $\mathrm{curl}:V_h\to N_h^2$ and of $\mathrm{div}:N_h^2\to N_h^3$. The moment-commuting supplier gives $DC=0$ on the full spaces, hence also after restriction to $V_h$. Set $j_i(t)={\int }_{\mathrm{\Omega }}{J}_{\mathrm{s}\mathrm{r}\mathrm{c}}(t)\cdot {\varphi }_{i}dx$, with $J_{\mathrm{s}\mathrm{r}\mathrm{c}}\in C([0,T];L^2(\mathrm{\Omega };{\mathbb{R}}^3))$. The lossless, conductivity-free semidiscretization with this source is $M_e\dot{e}-C^T{M}_{b}b=-j,\dot{b}+Ce=0.$ The coefficient ODE is defined for arbitrary initial coefficients. To represent the magnetic constraint of the Maxwell model, require $Db(0)=0$: the Faraday row and $DC=0$ give $D\dot{b}=-DCe=0$, so $Db(t)=Db(0)$ and compatible data remain divergence-free. An incompatible initial divergence is preserved, not corrected by this ODE. This is a spatial Galerkin system only; no time integrator is included. The cochain convention and weak material-interface formulation come from the declared spaces and Transient Maxwell system with piecewise materials and boundary conditions. Smooth interpolation and a uniformly bounded commuting projection are separate prerequisites for the later error estimate and are not asserted by this definition.

Choose a finite complex and boundary closure as in Oriented Maxwell cochains, masses, and the Yee grid, step size $\mathrm{\Delta }t>0$, initial $e^0\in {\mathbb{R}}^{N_e}$ and $b^{-1\mathrm{/}2}\in {\mathbb{R}}^{N_f}$, and source rows $j^{n+1\mathrm{/}2}\in {\mathbb{R}}^{N_e}$. The output is a sequence $e^n$ at integer times and $b^{n+1\mathrm{/}2}$ at half-integer times. For each $n\ge 0$, perform the ordered instructions

1. Update the face flux by $b^{n+1\mathrm{/}2}=b^{n-1\mathrm{/}2}-\mathrm{\Delta }tC{e}^n\mathrm{.}$

2. Solve the electric row $(M_e+\frac{\mathrm{\Delta }t}{2}S)e^{n+1}=(M_e-\frac{\mathrm{\Delta }t}{2}S)e^n+\mathrm{\Delta }t{C}^T{M}_b{b}^{n+1\mathrm{/}2}-\mathrm{\Delta }t{j}^{n+1\mathrm{/}2}\mathrm{.}$

For $S=0$ this is the lossless Yee leapfrog scheme. For $S⪰0$, the midpoint term is the specified trapezoidal treatment of conductivity. The row is well defined for every input because its matrix is invertible by The midpoint-conductivity Yee electric solve is well defined. This procedure makes no stability claim; the lossless stability condition is proved separately.

Assume the Countable Choice branch inherited from Edge-face Maxwell finite element semidiscretization. Let $T\ge 0$. For any finite matching mesh, initial edge and face coefficients $(e_0,b_0)$, and continuous load vector $j:[0,T]\to {\mathbb{R}}^{N_e}$, the source-driven semidiscrete Maxwell system $M_e\dot{e}-C^T{M}_{b}b=-j(t),\dot{b}+Ce=0,(e(0),b(0))=(e_0,b_0)$ has one and only one global $C^1$ solution on $[0,T]$ when $T>0$, with domain-relative one-sided endpoint derivatives. When $T=0$, there is exactly one map with the prescribed initial value; it extends locally to the same coefficient ODE with the load held constant at $j(0)$, and any two such continuations agree on their common interval containing zero. No derivative or IVP solution on the singleton itself is asserted. The construction and uniqueness argument use finite-dimensional linear algebra and the matrix-ODE supplier; no choice axiom stronger than the inherited Countable Choice branch is used.

Assume Countable Choice in the inherited edge/face finite-element space branch. Let $V=H_0(\mathrm{curl};\mathrm{\Omega })$ with the PEC convention in Edge-face Maxwell finite element semidiscretization, and let $V_h=N_h^1\cap H_0(\mathrm{curl};\mathrm{\Omega })$. Define $m(u,v)={\int }_{\mathrm{\Omega }}{v}^Tϵudx,a(u,v)={\int }_{\mathrm{\Omega }}(\mathrm{curl}v{)}^T{\mu }^{-1}\mathrm{curl}udx,a_{+}(u,v)=m(u,v)+a(u,v)\mathrm{.}$ For $u\in V$, the Maxwell curl-energy Ritz projection $R_{h}u\in V_h$ is the solution of $a_{+}(R_{h}u,v_h)=a_{+}(u,v_h)(v_h\in V_h)\mathrm{.}$ Its finite-dimensional well-definedness, orthogonality, and best-approximation property are supplied by Well-definedness and best approximation of the Maxwell Ritz map ↗.

Assume Countable Choice only in the H(div)/Piola space branch inherited from the suppliers below. Let ${\mathcal{T}}_h$ be a finite matching affine tetrahedral mesh of a bounded polyhedron, with shape-regularity constant $\sigma$ and $h={\max }_K{h}_K$. Assume material interfaces are mesh-aligned. Let $u,v$ be piecewise $C^1$ on cell closures, with $\mathrm{curl}u$, $\mathrm{div}v$ also piecewise $C^1$, and with matching tangential traces of $u$ and matching normal traces of $v$. The canonical edge and face moment interpolants $I_h^1,I_h^2$ are those in Whitney edge, face, and cell moments commute with grad, curl, and div. Then $\begin{array}{rl}{\displaystyle }& {\displaystyle \mathrm{\parallel }u-I_h^{1}u{\mathrm{\parallel }}_{L^2(\mathrm{\Omega })}+\mathrm{\parallel }\mathrm{curl}(u-I_h^{1}u){\mathrm{\parallel }}_{L^2(\mathrm{\Omega })}+\mathrm{\parallel }v-I_h^{2}v{\mathrm{\parallel }}_{L^2(\mathrm{\Omega })}+\mathrm{\parallel }\mathrm{div}(v-I_h^{2}v){\mathrm{\parallel }}_{L^2(\mathrm{\Omega })}\le C_{\sigma }h\mathrm{\mid }\mathrm{\Omega }{\mathrm{\mid }}^{1\mathrm{/}2}\underset{K\in {\mathcal{T}}_h}{\max }(\mathrm{\parallel }\mathrm{\nabla }u{\mathrm{\parallel }}_{L^{\mathrm{\infty }}(K)}+\mathrm{\parallel }\mathrm{\nabla }\mathrm{curl}u{\mathrm{\parallel }}_{L^{\mathrm{\infty }}(K)}+\mathrm{\parallel }\mathrm{\nabla }v{\mathrm{\parallel }}_{L^{\mathrm{\infty }}(K)}+\mathrm{\parallel }\mathrm{\nabla }\mathrm{div}v{\mathrm{\parallel }}_{L^{\mathrm{\infty }}(K)})\mathrm{.}}\end{array}$ The constant depends only on reference-element basis norms and $\sigma$. This estimate is only for the canonical smooth-input interpolants; it makes no boundedness claim on all of $H(\mathrm{curl})$ or $H(\mathrm{div})$.

For the Yee update in Lossless or midpoint-conductive Yee leapfrog update, suppose the initial rows satisfy $-G^T{M}_e{e}^0={\rho }^0$ and $D{b}^{-1\mathrm{/}2}=0$. If the source and charge rows satisfy ${\rho }^{n+1}-{\rho }^n-\mathrm{\Delta }t{G}^T(j^{n+1\mathrm{/}2}+S\frac{e^{n+1}+e^n}{2})=0$ at every step, then $-G^T{M}_e{e}^n={\rho }^n$ and $D{b}^{n-1\mathrm{/}2}=0$ for every later step. This is exact constraint propagation; it does not assert stability for a time step violating the CFL condition.

Set $K=C^T{M}_{b}C$. If the retained electric cochain space is zero-dimensional, define ${\lambda }_{\max }:=0$; otherwise set ${\lambda }_{\max }:={\max }_{e\mathrm{\ne }0}\frac{e^{T}Ke}{e^T{M}_{e}e}\mathrm{.}$ For a lossless, source-free update from Lossless or midpoint-conductive Yee leapfrog update, define ${\mathcal{E}}^n:=\frac{1}{2}(e^n{)}^T{M}_e{e}^n+\frac{1}{2}(b^{n-1\mathrm{/}2}{)}^T{M}_b{b}^{n+1\mathrm{/}2},b_c^n:=\frac{b^{n-1\mathrm{/}2}+b^{n+1\mathrm{/}2}}{2}\mathrm{.}$ Then ${\mathcal{E}}^n$ is independent of $n$, and ${\mathcal{E}}^n\ge \frac{1}{2}(1-\frac{\mathrm{\Delta }{t}^2{\lambda }_{\max }}{4})\mathrm{\parallel }{e}^n{\mathrm{\parallel }}_{M_e}^2+\frac{1}{2}\mathrm{\parallel }{b}_c^n{\mathrm{\parallel }}_{M_b}^2\mathrm{.}$ Consequently, the actual strict spectral condition $\mathrm{\Delta }{t}^2{\lambda }_{\max }<4$ makes this modified energy positive on admissible staggered states and controls the electric field and the centered magnetic field. When ${\lambda }_{\max }>0$, at equality the energy loses coercivity and a maximal mode has a repeated root $-1$; above equality a mode grows exponentially. If ${\lambda }_{\max }=0$, the strict condition is automatic for every $\mathrm{\Delta }t>0$ and neither endpoint case occurs.

Let a periodic rectangular cell be discretized by a uniform Cartesian Yee grid with spacings $h_x,h_y,h_z$, $h={\max }_j{h}_j$, and constant scalar coefficients $ϵ,\mu >0$. Let $(E,B)$ be a source-free smooth periodic Maxwell solution on $[0,T]$, with the spatial derivatives of $E$ through total order four, ${\mathrm{\partial }}_t^{4}E$, and ${\mathrm{\partial }}_t^{3}B$ bounded. Sample $E(t_n)$ by exact oriented edge line integrals $e_{\ast }^n$ and $B(t_{n+1\mathrm{/}2})$ by exact oriented face fluxes $b_{\ast }^{n+1\mathrm{/}2}$. Define $r_h^n=M_e\frac{e_{\ast }^{n+1}-2e_{\ast }^n+e_{\ast }^{n-1}}{\mathrm{\Delta }{t}^2}+C^T{M}_{b}C{e}_{\ast }^n,s_h^n=b_{\ast }^{n+1\mathrm{/}2}-b_{\ast }^{n-1\mathrm{/}2}+\mathrm{\Delta }tC{e}_{\ast }^n\mathrm{.}$ With the homogeneous diagonal Hodge weights of Oriented Maxwell cochains, masses, and the Yee grid, there is a constant $C_T$ independent of the cell counts such that, whenever $\mathrm{\Delta }t\le t_n\le T-\mathrm{\Delta }t$, $\mathrm{\parallel }{r}_h^n{\mathrm{\parallel }}_{M_e^{-1}}\le C_T(h^2+\mathrm{\Delta }{t}^2),\mathrm{\parallel }{s}_h^n{\mathrm{\parallel }}_{M_b}\le C_T\mathrm{\Delta }{t}^3\mathrm{.}$ Here $\mathrm{\parallel }r{\mathrm{\parallel }}_{M_e^{-1}}^2=r^T{M}_e^{-1}r$. This is a smooth homogeneous periodic consistency result; it asserts no Taylor estimate across coefficient jumps or a PEC boundary closure.