Study Guide: Adaptive FEM for the Wave Equation Grote, Lakkis & Santos (2026)

Study Guide

Adaptive FEM with Explicit Time Integration

for the Wave Equation
M. J. Grote, O. Lakkis & C. S. Santos
Journal of Computational and Applied Mathematics 481 (2026) 117272

Abstract of this guide. This document provides a self-contained background on every

concept required to understand and present the Grote–Lakkis–Santos paper. Topics covered:
the second-order wave equation, finite element discretisation, the Galerkin formulation,
the leapfrog time integrator, the CFL stability condition, a posteriori error estimation,
compatible adaptive mesh strategies, and leapfrog-based local time-stepping (LF-LTS).

1. The Wave Equation 1.2 Why is the wave equation hard to

solve adaptively?
A partial differential equation (PDE) re- • No dissipation: unlike parabolic prob-
lates an unknown function of several variables lems, errors accumulate over time; there is
to its partial derivatives. The wave equation no natural damping.
is a second-order hyperbolic PDE modelling • Explicit methods dominate: leapfrog
how disturbances (acoustic, seismic, electro- (central differences) is cheap and conserves
magnetic) propagate through a medium with energy, but requires a CFL time-step re-
spatially varying wave speed c(x). striction.
• Adaptive meshes + explicit methods
1.1 Governing equation conflict: local refinement forces a tiny
Let Ω ⊂ Rd (d = 1, 2, 3) be a bounded Lip- global time-step, destroying efficiency (see
schitz domain and T > 0 a final time. The §5).
unknown displacement u : Ω × (0, T ] → R
satisfies 2. Finite Element Method
∂ 2u
−∇· c2 (x) ∇u = f (x, t),
 
2
x ∈ Ω, t ∈ (0, T ],
∂t 2.1 Mesh and finite element space
(1)
with 0 < cmin ≤ c(x) ≤ cmax , forcing f ∈ At each discrete time tn = n∆t, the domain
L2 (0, T ; L2 (Ω)), and initial conditions Ω is partitioned into a mesh Mn of elements
K (triangles in 2-D, tetrahedra in 3-D). The
∂u mesh-size function is
u(x, 0) = u0 (x), (x, 0) = v0 (x). (2)  \ 
∂t hn (x) := diam K̄ . (4)
Homogeneous boundary conditions are im- x∈K∈Mn
posed on ∂Ω: The associated finite element (FE) space
∂u is
u = 0 on Γ0 , = 0 on ∂Ω \ Γ0 , (3) Vn = Pk (Mn ) ∩ C 0 (Ω), (5)
∂n
the space of continuous piecewise polynomi-
where Γ0 ⊂ ∂Ω has positive measure.
als of degree k ≥ 1. With dimension Mn =
Well-posedness. For u0 ∈ HΓ10 (Ω),
dim Vn , any V ∈ Vn expands as
2
v0 ∈ L (Ω), a unique solution u ∈ Mn
X
C 0 (0, T ; HΓ10 (Ω)) ∩ C 1 (0, T ; L2 (Ω)) exists. V (x) = vm Φnm (x), (6)
m=1

1
Study Guide: Adaptive FEM for the Wave Equation Grote, Lakkis & Santos (2026)

where {Φnm } are the nodal (“hat”) basis func- nique), enabling purely explicit updates.
tions satisfying Φnm (znℓ ) = δm,ℓ .
2.2 Transfer operator 3.3 Energy norms
When the mesh changes between steps, the Two norms measure the accuracy of the FE
FE transfer operator Πn : C 0 (Ω) → Vn (La- solution. The potential energy norm is
grangian interpolation or L2 -projection) maps p
solutions to the new space: ∥ϕ∥A := a(ϕ, ϕ), (11)
Mn
X and the wave energy norm for a pair ϕ =
v znm Φnm (x).


Πn v(x) = (7)
m=1
[ϕ0 , ϕ1 ]⊤ is
q
2.3 Coarse and fine subspaces ∥ϕ∥erg := ∥ϕ0 ∥2A + ∥ϕ1 ∥2L2 (Ω) . (12)
Each mesh decomposes into Mn = Mcn ∪
Mfn (coarse and fine parts), giving a splitting These norms are physically meaningful: the
Vn = Vnf ⊕ Vnc . The fine-mesh interpolator first measures potential (strain) energy, the
Πfn : C 0 (Ω) → Vnf is used inside the local second total mechanical energy.
time-stepping operator (see §8).
4. Leapfrog Time Integration
3. Galerkin Formulation

3.1 Bilinear form and elliptic operator 4.1 The leapfrog (Störmer–Verlet)
scheme
The elliptic bilinear form associated with (1)
is The leapfrog method is an explicit, second-
Z order time integrator in which displacement
a(u, v) = c2 (x) ∇u·∇v dx, u, v ∈ H 1 (Ω). U and velocity V are staggered in time:
Ω
(8) U evaluated at integer times tn = n∆t,
For each FE space Vn we define the discrete 1


elliptic operator An : Vn → Vn via the Riesz V evaluated at half-integer times tn+ 1 = n + 2
∆t.
2

representation
Including mesh change and the LTS operator
⟨An ϕ, v⟩ = a(ϕ, v) ∀ v ∈ Vn . (9) Aen (defined in §8), the full scheme reads:

In matrix form, An corresponds to the stiffness U 0 := Π0 u0 , (13a)

matrix K; the L2 inner product corresponds h i
U 1 := Π1 U 0 + Π0 v0 ∆t + F 0 − A
e0 U 0 ∆t2 ,


to the mass matrix M .
(13b)
3.2 Semi-discrete system h i
U n+1 := Πn+1 2U n − Πn U n−1 + F n − A
en U n ∆t2 ,


Galerkin projection onto Vn converts (1) into
(13c)
M Ü(t) + K U(t) = F(t), (10) n+1 n
1 U − Πn+1 U
V n+ 2 := . (13d)
a second-order ODE system in the nodal coef- ∆t
ficient vector U(t).
The centered difference approximations are
Mass lumping. For explicit time integration,
M must be diagonal. Mass lumping replaces U n+1 − U n−1 U n+1 − 2U n + U n−1
∂U n := , ∂ 2 U n := .
M by a diagonal matrix M̂ (row-sum tech- 2∆t ∆t2
(14)

2
Study Guide: Adaptive FEM for the Wave Equation Grote, Lakkis & Santos (2026)

4.2 Advantages of explicit methods 6. A Posteriori Error Estimation

• No linear system solve: O(Mn ) work
per step.
• Natural parallelism: only local matrix–
6.1 A priori vs. a posteriori
vector products.
• Energy conservation: for f ≡ 0 and a An a posteriori error estimate bounds the
fixed mesh, leapfrog conserves a discrete error ∥U n − u(tn )∥ using only the computed
energy exactly. solution—no knowledge of the exact solution is
required. This contrasts with a priori bounds,
5. CFL Stability Condition which require regularity information about the
exact solution.
5.1 The condition
6.2 Main estimate
For explicit time integration, the Courant–
Friedrichs–Lewy (CFL) condition is a nec- The paper derives the fully computable bound
essary stability condition:

hmin (
∆t ≤ CCFL , (15) 2N
cmax
X
n
max ∥U − u(tn )∥erg ≤ ∥e(0)∥erg +C 2 ζm + max
0≤n≤N 1≤n≤N
m=1
where hmin is the smallest element diameter (17)
and CCFL depends on the element shape and where all terms on the right are computable
polynomial degree. Violation causes exponen- from the numerical solution.
tial blow-up of the numerical solution.
5.2 The bottleneck with local refine- 6.3 Error indicator components
ment
Elliptic error indicators
If the mesh is locally refined so that hfine =
hcoarse /p in a small sub-region, then (15) forces The residual-based elliptic indicator on mesh
Mn is
hfine hcoarse
∆t ≤ C = C , (16)
cmax p cmax
X n 2
imposing p× more time steps globally, even E[w, Vn ]2 := h2K An w−∇·(c2 ∇w) +1
L2 (K) 2
hK c
though the fine region may occupy only a tiny K∈Mn

fraction of Ω. (18)
The two terms measure the interior residual
Classic remedy. The standard answer has and the inter-element flux jump, respectively.
been to use unconditionally stable implicit
methods. The Grote–Lakkis–Santos paper Time accumulation indicators
challenges this by introducing explicit LF-LTS:
small steps only where needed.
The indicator ζm accumulates contributions
from time discretisation, mesh change, LTS,
5.3 Implicit vs. explicit comparison and data approximation:

Property Implicit Explicit+LTS

Stability Unconditional Local CFL only
Z t(m+1)/2 h
2
2 i1
Cost/step O(Mn1.5 ) O(Mn ) ζm := µn0 +ϑn0 (t) + αn +µn1 +δ n (t)+ϑn1 (t)
tm/2
Local refine. No restriction Efficient
(19)
Parallelism Limited Excellent

3
Study Guide: Adaptive FEM for the Wave Equation Grote, Lakkis & Santos (2026)

Summary of all indicators In 1-D with piecewise linear elements,

h i
Symbol Meaning y − ŷ = ym − 2 (ym−1 + ym+1 ) Φnm (x), (22)
1

εn0 , εn1 Elliptic residual indicators

n n
which vanishes when ym equals the average of
µ0 , µ1 , Mesh-change indicators
n its neighbours. Coarsening is permitted only
µ2 (nonzero only when mesh 0 1
where max(βK , βK ) < tolC .
changes)
α n LTS error (using A en instead
8. Local Time-Stepping (LTS)
of An )
n
ϑ0 , ϑ 1 n Time-discretisation error
(O(∆t2 )) 8.1 Motivation and idea
δ n (t) Data approximation LTS resolves the CFL bottleneck: instead of
(∥F n − f (t)∥L2 ) using the small time-step ∆τ = ∆t/p every-
where, it applies p sub-iterations only inside
6.4 Dörfler (bulk) marking the fine sub-region Ωfn while the coarse region
Elements are marked for refinement using the advances with the global step ∆t.
Dörfler criterion with parameter θ ∈ (0, 1): 8.2 The modified operator
select the smallest set M such that The p sub-steps are encoded in the modified
X
2
X
2 elliptic operator
ηK ≥ θ ηK . (20)
en := An Pp ∆t2 Πfn An ,


K∈M K∈Mn A (23)
The paper uses θ = 0.8. where the polynomial Pp is defined via the
Chebyshev polynomial Tp of degree p:
7. Adaptive Mesh Strategy   
2 x
Pp (x) = 1 − Tp 1 − 2 , Tp (x) = 2x Tp−1 (x)−
x 2p
7.1 Compatible meshes (24)
2
Two meshes M1 and M2 are compatible if For p = 2: T2 (x) = 2x − 1, giving P2 (x) =
every element of M1 is either an element of 1 − x/16 and Aen = An − ∆t2 An Πf An .
16 n
M2 or a union thereof, and vice versa. Com- 8.3 Stabilized version LF-LTS(ν)
patibility is maintained by using bisection re-
To guarantee stability for all values of p and
finement from a common macro triangulation
∆t, the Chebyshev polynomial is damped by a
M0 .
small parameter ν > 0 (typically ν = 0.01):
• Refinement (K → K− ∪ K+ ): interpola-  
2 Tp (δp,ν − x/ωp,ν )
tion to new nodes is exact for polynomials Pp,ν (x) = 1− , (25)
⇒ no additional error. x Tp (δp,ν )
• Coarsening (K− , K+ → K̂): removing with δp,ν = 1 + ν/p2 and ωp,ν =
the midpoint node introduces a small error, 2Tp′ (δp,ν )/Tp (δp,ν ). Setting ν = 0 recovers the
controlled by the coarsening pre-indicators original LF-LTS(0).
βK .
8.4 Why Chebyshev polynomials?
7.2 Coarsening pre-indicators • Chebyshev polynomials minimise the maxi-
Before coarsening elements K− , K+ ⊂ K̂, the mum error on an interval (optimal polyno-
potential information loss is estimated: mial approximation).
• The p sub-iterations emerge naturally from
0 1
βK := (Πn+1 −Id) y A,K , βK := (Πn+1 −Id) y L2 (K)
the. three-term recurrence—no explicit sub-
(21) loop is needed.

4
Study Guide: Adaptive FEM for the Wave Equation Grote, Lakkis & Santos (2026)

• Only unknowns in Vnf are touched during Parameter Value Role

sub-stepping (the Πfn projector acts as a
tolH 20 Refinement threshold
mask).
tolC 10−4 (1-D); 0.01 (2-D) Coarsening threshold
• Multiple refinement levels are handled by θ 0.8 Dörfler fraction
nesting, each with its own local p. ν 0.01 LTS stabilisation
8.5 Computational savings

Experiment Memory Compute

1-D traveling wave ≤ 5% of dofs ≈ 0.95%
1-D splitting wave ≤ 6% of dofs ≈ 1.2%
2-D L-shaped domain significantly reduced significantly reduced

Key result. The adaptive LF-LTS algorithm

achieves over 100× reduction in computa-
tional effort compared to a uniform mesh us-
ing the finest mesh size hfine .

9. The Full Adaptive Algorithm

9.1 Initialisation
1. Start with macro mesh M0 .
0
2. Estimate ηK = hK ∥∂xx u0 ∥L2 (K) and analo-
gously for v0 .
3. Refine via Dörfler marking (20) ⇒ initial
mesh M1 .
4. Set U 0 = Π1 u0 ; compute U 1 by Taylor ex-
pansion.
9.2 Time-stepping loop
The experiments confirm optimal conver-
At each step n → n + 1, the algorithm pro-
gence: O(h) in the energy norm and O(h2 ) in
ceeds:
L2 (Ω). The a posteriori estimator converges
at the same rate O(h) as the true error, con-
1. Tentative step. Set Mn+1 = Mn ; compute
firming reliability and efficiency.
U n+1 via (13c) using LF-LTS.
2. Error check. Evaluate εn+1 0 , εn+1
1 . If
n+1
max(εi ) > tolH /N : refine Mn+1 , trans-
fer solutions, repeat.
3. Accept. Accept U n+1 on current mesh.
4. Coarsen. Compute βK ; coarsen elements 10. Notation Reference
0 1
with max(βK , βK ) ≤ tolC .
5. Transfer. Move U n , U n+1 to updated Vn+1 .
9.3 Parameters and convergence

5
Study Guide: Adaptive FEM for the Wave Equation Grote, Lakkis & Santos (2026)

Symbol Meaning Explicit Time Integration for the Wave Equation

– J. Comput. Appl. Math. 481 (2026) 117272.
Ω ⊂ Rd Spatial domain
u(x, t) Exact solution
Un FE approx. at tn = n∆t
V n+1/2 FE velocity approx.
c(x) Wave speed
f (x, t) Source term
Mn FE mesh at step n
hK Diameter of element K
Vn FE space at step n
Mn = Degrees of freedom
dim Vn
Φnm m-th nodal basis function
Πn FE transfer operator
An Discrete elliptic operator
Aen LTS-modified operator
a(u, v) Bilinear form (8)
∥ · ∥A Potential energy norm
∥ · ∥erg Wave energy norm
∆t Global time step
∆τ = ∆t/p Local sub-step
p Coarse/fine step ratio
Tp (x) Chebyshev poly. of order p
εn0 , εn1 Elliptic indicators
ζm Time accumulation indicator
µni Mesh-change indicators
αn LTS error indicator
ϑni Time-error indicators
δ n (t) Data approximation indicator
βK Coarsening pre-indicator
tolH Refinement tolerance
tolC Coarsening tolerance
θ Dörfler fraction
ν LTS stabilisation parameter

Functional spaces

L2 (Ω) : square-integrable functions,

H 1 (Ω) : functions with L2 gradient (Sobolev),
HΓ10 (Ω) : H 1 functions vanishing on Γ0 ,
L2 (0, T ; X) : Bochner space (X-valued L2 functions).

Grote, Lakkis & Santos – Adaptive FEM with