Society for Industrial and Applied Mathematics: SIAM Journal on Numerical Analysis: Table of Contents

Unfitted Hybrid High-Order Methods Stabilized by Polynomial Extension for Elliptic Interface Problems

Erik Burman — 2026-05-13T07:00:00Z

SIAM Journal on Numerical Analysis, Volume 64, Issue 3, Page 601-630, June 2026.
Abstract. In this work, we study the design and analysis of a novel hybrid high-order (HHO) method on unfitted meshes. HHO methods rely on a pair of unknowns, combining polynomials attached to the mesh faces and the mesh cells. In the unfitted framework, the interface can cut through the mesh cells in a very general fashion, and the polynomial unknowns are doubled in the cut cells and the cut faces. In order to avoid the ill-conditioning issues caused by the presence of small cut cells, the novel approach introduced herein is to use polynomial extensions in the definition of the gradient reconstruction operator. Stability and consistency results are established, leading to optimally decaying error estimates. The theory is illustrated by numerical experiments.

A Novel Augmented Subspace Adaptive Finite Element Method for the Eigenvalue Problem with Discontinuous Coefficients

Manting Xie — 2026-05-18T07:00:00Z

SIAM Journal on Numerical Analysis, Volume 64, Issue 3, Page 631-655, June 2026.
Abstract. This paper proposes an efficient adaptive finite element method (AFEM) for solving the eigenvalue problem with discontinuous coefficients. Different from the existing adaptive algorithms for eigenvalue problems, the method innovatively focuses on solving linearized boundary value problems in each adaptively refined space, complemented by solving small-scale eigenvalue problems on low-dimensional augmented subspaces which are automatically controlled by the algorithm. Notably, it does not require solving the small-scale eigenvalue problem at every iteration, which improves the computational efficiency. Moreover, a novel a posteriori error estimator, which relies on the local oscillations of coefficients near singular points, guides the adaptive refinement process. To further substantiate this method, we give a thorough and rigorous convergence analysis in this paper. Finally, two numerical examples are provided to illustrate the accuracy and efficiency of our new AFEM.

Lévy Score Function and Score-Based Particle Algorithm for Nonlinear Lévy–Fokker–Planck Equations

Yuanfei Huang — 2026-05-20T07:00:00Z

SIAM Journal on Numerical Analysis, Volume 64, Issue 3, Page 656-684, June 2026.
Abstract. The score function for the diffusion process, also known as the gradient of the log-density, is a basic concept to characterize the probability flow with important applications in the score-based diffusion generative modeling and the simulation of Itô stochastic differential equations. However, neither the probability flow nor the corresponding score function for the diffusion-jump process is known. This paper delivers mathematical derivation, a numerical algorithm, and error analysis focusing on the corresponding score function in non-Gaussian systems with jumps and discontinuities represented by the nonlinear Lévy–Fokker–Planck equations. We propose the Lévy score function for such stochastic equations, which features a nonlocal double-integral term, and we develop its training algorithm by minimizing the proposed loss function from samples. Based on the equivalence of the probability flow with deterministic dynamics, we develop a self-consistent score-based transport particle algorithm to sample the interactive Lévy stochastic process at discrete-time grid points. We provide an error bound for the Kullback–Leibler divergence between the numerical and true probability density functions by overcoming the nonlocal challenges in the Lévy score. The full error analysis with the Monte Carlo error and the time discretization error is furthermore established. To show the usefulness and efficiency of our approach, numerical examples from applications in biology and finance are tested.

A Structure-Preserving and Thermodynamically Compatible Cell-Centered Lagrangian Finite Volume Scheme for Continuum Mechanics

Walter Boscheri — 2026-05-22T07:00:00Z

SIAM Journal on Numerical Analysis, Volume 64, Issue 3, Page 685-707, June 2026.
Abstract. In this work we present a novel structure-preserving scheme for the compatible discretization of the Godunov–Peshkov–Romenski (GPR) model of continuum mechanics written in Lagrangian form. The governing equations fall into the larger class of overdetermined hyperbolic and thermodynamically compatible systems of partial differential equations (PDE). This model admits an extra conservation law for the total energy (first principle of thermodynamics) and satisfies the entropy inequality (second principle of thermodynamics). Furthermore, in the absence of algebraic source terms, the distortion field of the continuum and the specific thermal impulse satisfy a curl-free condition, provided the initial data are curl-free. Last but not least, the determinant of the distortion field is related to the density of the medium, i.e., the system is also endowed with a nonlinear algebraic constraint. In the stiff relaxation limit, the system tends to the compressible Navier–Stokes equations, i.e., the GPR model is able to describe at the same time the dynamics of nonlinear solids as well as the one of fluids. The objective of this work is to construct and analyze a new semidiscrete thermodynamically compatible cell-centered Lagrangian finite volume scheme on moving unstructured meshes that satisfies the following structural properties of the governing PDE exactly at the discrete level: (i) compatibility with the first law of thermodynamics, i.e., discrete total energy conservation; (ii) compatibility with the second law of thermodynamics, i.e., discrete entropy inequality; (iii) exact discrete compatibility between the density and the determinant of the distortion field; (iv) exact preservation of the curl-free property of the distortion field and of the specific thermal impulse in the absence of algebraic source terms. We will show that it is indeed possible to achieve all the above properties simultaneously. Unlike in existing schemes, we choose to directly discretize the entropy inequality, hence obtaining total energy conservation as a consequence of an appropriate and thermodynamically compatible discretization of all the other equations. From this discretization, property (ii) above follows trivially by construction, while (i) leads to provable nonlinear stability, which is an important feature, in particular for the complex PDE system under consideration here. The thermodynamic compatibility and thus nonlinear stability in the sense of total energy conservation is achieved via a very simple and general approach recently introduced by Abgrall, Dumbser, and Maire. by using a scalar correction factor that is defined at the nodes of the grid. This perfectly fits into the formalism of nodal solvers which is typically adopted in cell-centered Lagrangian finite volume methods. The new scheme is run on some academic benchmark problems for computational fluid and solid mechanics to show that the properties also hold in the practical implementation of the scheme.

Uniformly HP-Stable Elements for the Elasticity Complex

Francis R. A. Aznaran — 2026-05-27T07:00:00Z

SIAM Journal on Numerical Analysis, Volume 64, Issue 3, Page 708-736, June 2026.
Abstract. For the discretization of symmetric, divergence-conforming stress tensors in continuum mechanics, we prove inf-sup stability bounds which are uniform in polynomial degree and mesh size for the Hu–Zhang finite element in two dimensions. This is achieved via an explicit construction of a bounded right inverse of the divergence operator, with the crucial component being the construction of bounded Poincaré operators for the stress elasticity complex which are polynomial-preserving, in the Bernstein–Gelfand–Gelfand framework of the finite element exterior calculus. We also construct [math]-bounded projection operators satisfying a commuting diagram property and [math]-stable Hodge decompositions. Numerical examples are provided.

On the Convergence of Split Exponential Integrators for Semilinear Parabolic Problems

Marco Caliari — 2026-05-27T07:00:00Z

SIAM Journal on Numerical Analysis, Volume 64, Issue 3, Page 737-754, June 2026.
Abstract. Splitting the exponential-like [math] functions, which typically appear in exponential integrators, is attractive in many situations since it can dramatically reduce the computational cost of the procedure. However, depending on the employed splitting, this can result in order reduction. The aim of this paper is to analyze different such split approximations. We perform the analysis for semilinear problems in the abstract framework of commuting semigroups and derive error bounds that depend, in particular, on whether the vector (to which the [math] functions are applied) satisfies appropriate boundary conditions. We then present the convergence analysis for two split versions of a second-order exponential Runge–Kutta integrator in the context of analytic semigroups, and show that one suffers from order reduction while the other does not. Numerical results for semidiscretized parabolic PDEs confirm the theoretical findings.

Numerical Analysis for Saddle Dynamics of Some Semilinear Elliptic Problems

Lei Zhang — 2026-05-29T07:00:00Z

SIAM Journal on Numerical Analysis, Volume 64, Issue 3, Page 755-778, June 2026.
Abstract. This work presents a numerical analysis of computing transition states of semilinear elliptic partial differential equations (PDEs) via index-1 saddle dynamics or, equivalently, the gentlest ascent dynamics. To establish clear connections between saddle dynamics and numerical methods of PDEs, as well as improving their compatibility, we first propose the continuous-in-space formulation of saddle dynamics for semilinear elliptic problems. This formulation yields a parabolic system that converges to saddle points. We then analyze the well-posedness, [math] stability, and error estimates of semidiscrete and fully discrete finite element schemes. Significant efforts are devoted to addressing the coupling, gradient nonlinearity, nonlocality of the proposed parabolic system, and the impacts of retraction due to the norm constraint. The error estimate results demonstrate the accuracy and index-preservation of the discrete schemes.

Error Analysis of a Conforming Finite Element Method for the Modified Electromagnetic Transmission Eigenvalue Problem

Jiayu Han — 2026-03-03T08:00:00Z

SIAM Journal on Numerical Analysis, Volume 64, Issue 2, Page 303-325, April 2026.
Abstract. The modified electromagnetic transmission eigenvalue problem (METEP) arises from the inverse scattering theory and can be used to detect changes of the material properties in nondestructive testing. This paper proposes and analyzes a conforming edge element method for the METEP. We establish a rigorous error analysis of the numerical eigenpairs by proving the uniform convergence of the discrete operator. In particular, as the problem contains two second order equations and is indefinite, we introduce auxiliary problems and show that they satisfy [math]-coercivity, based on which we prove the existence of both the continuous and discrete solution operators to the source problem. We then prove the uniform convergence of the discrete solution operator by reformulating the continuous and discrete solution operators. Optimal error estimates are obtained by investigating the adjoint problems and using the spectral approximation theory for compact operators. The theory is validated by numerical examples with various coefficients for different domains in both two and three dimensions.

Preasymptotic Error Estimates of Linear EEM and CIP-EEM for the Time-Harmonic Maxwell Equations with Large Wave Number

Shuaishuai Lu — 2026-03-03T08:00:00Z

SIAM Journal on Numerical Analysis, Volume 64, Issue 2, Page 326-349, April 2026.
Abstract. Preasymptotic error estimates are derived for the second-type Nédélec linear edge element method and the linear [math]-conforming interior penalty edge element method (CIP-EEM) for the time-harmonic Maxwell equations with large wave number. It is shown that under the mesh condition that [math] is sufficiently small, the errors of the solutions to both methods are bounded by [math] in the energy norm and [math] in the [math]-scaled [math] norm, where [math] is the wave number and [math] is the mesh size. Numerical tests are provided to illustrate our theoretical results and the potential of CIP-EEM in significantly reducing the pollution effect.

A Posteriori Error Control for Nonconvex Problems via Calibration

Benjamin Berkels — 2026-03-04T08:00:00Z

SIAM Journal on Numerical Analysis, Volume 64, Issue 2, Page 350-369, April 2026.
Abstract. In this paper, a posteriori error estimates are derived for the approximation error of minimizers of functionals on the space of functions with bounded variation with a nonconvex lower-order term. To this end, the calibration method by Alberti, Bouchitté, and Dal Maso [Calc. Var. Partial Differential Equations, 16 (2003), pp. 299–333] allows the problem to be reformulated as a uniformly convex variational problem over characteristic functions of subgraphs in one dimension higher. A primal-dual approach is formulated where the duality of divergence and gradient properly incorporates boundary conditions for the primal variable. Based on this, a posteriori error estimates can be derived first for the relaxed problem in the [math]-norm. A cut-out argument allows converting this into an [math]-error estimate for the characteristic subgraph functions apart from the jump interface, whereas the area of the interfacial region is estimated separately. To apply the estimate, we consider as one possible discretization a conforming finite element space for the primal variable and a nonconforming space for the dual variable. Finally, we validate the a posteriori error estimates in numerical experiments for a prototypical nonconvex functional in one and two dimensions as well as depth estimation in stereo imaging, a classical computer vision problem.

Local Time Integration for Friedrichs’ Systems

Marlis Hochbruck — 2026-03-09T07:00:00Z

SIAM Journal on Numerical Analysis, Volume 64, Issue 2, Page 370-390, April 2026.
Abstract. In this paper, we address the full discretization of Friedrichs’ systems with a two-field structure, such as Maxwell’s equations or the acoustic wave equation in div-grad form; cf. [W. Dörfler et al., Wave Phenomena: Mathematical Analysis and Numerical Approximation, Springer, Cham, 2023]. We focus on a discontinuous Galerkin space discretization applied to a locally refined mesh or a small region with high wave speed. This results in a stiff system of ordinary differential equations, where the stiffness is mainly caused by a small region of the spatial mesh. When using explicit time-integration schemes, the time stepsize is severely restricted by a few spatial elements, leading to a loss of efficiency. As a remedy, we propose and analyze a general leapfrog-based scheme which is motivated by [C. Carle and M. Hochbruck, SIAM J. Numer. Anal., 60 (2022), pp. 2897–2924]. The new, fully explicit, local time-integration method filters the stiff part of the system in such a way that its CFL condition is significantly weaker than that of the leapfrog scheme while its computational cost is only slightly larger. For this scheme, the filter function is a suitably scaled and shifted Chebyshev polynomial. While our main interest is in explicit local time-stepping schemes, the filter functions can be much more general, for instance, a certain rational function leads to the locally implicit method, proposed and analyzed in [M. Hochbruch and A. Sturm, SIAM J. Numer. Anal., 54 (2016), pp. 3167–3191]. Our analysis provides sufficient conditions on the filter function to ensure full order of convergence in space and second order in time for the whole class of local time-integration schemes.

Accuracy of the Ensemble Kalman Filter in the Near-Linear Setting

E. Calvello — 2026-03-17T07:00:00Z

SIAM Journal on Numerical Analysis, Volume 64, Issue 2, Page 391-429, April 2026.
Abstract. The filtering distribution captures the statistics of the state of a possibly stochastic dynamical system from partial and noisy observations. Classical particle filters provably approximate this distribution in quite general settings; however, they behave poorly for high dimensional problems, suffering weight collapse. This issue is circumvented by the ensemble Kalman filter, which is an equal-weights interacting particle system. However, this finite particle system is only proven to approximate the true filter in the linear Gaussian case. In practice, however, it is applied in much broader settings; as a result, establishing its approximation properties more generally is important. There has been recent progress in the theoretical analysis of the algorithm in discrete time, establishing stability and error estimates, in relation to the true filter, in non-Gaussian settings; but the assumptions on the dynamics and observation models rule out the unbounded vector fields that arise in practice, and the analysis applies only to the mean field limit of the discrete time ensemble Kalman filter. The present work establishes error bounds between the filtering distribution and the finite particle discrete time ensemble Kalman filter when the dynamics and observation vector fields may be unbounded, allowing linear growth.

Multipoint Stress Mixed Finite Element Methods for Elasticity on Cuboid Grids

Ibrahim Yazici — 2026-04-07T07:00:00Z

SIAM Journal on Numerical Analysis, Volume 64, Issue 2, Page 430-455, April 2026.
Abstract. We develop multipoint stress mixed finite element methods for linear elasticity with weak stress symmetry on cuboid grids, which can be reduced to a symmetric and positive definite cell-centered system. The methods employ the lowest-order enhanced Raviart–Thomas finite element space for the stress and piecewise constant displacement. The vertex quadrature rule is employed to localize the interaction of stress degrees of freedom, enabling local stress elimination around each vertex. We introduce two methods. The first method uses a piecewise constant rotation, resulting in a cell-centered system for the displacement and rotation. The second method employs a continuous piecewise trilinear rotation and the vertex quadrature rule for the asymmetry bilinear forms, allowing for further elimination of the rotation and resulting in a cell-centered system for the displacement only. Stability and error analysis are performed for both methods. For the stability analysis of the second method, a new auxiliary [math]-conforming matrix-valued space is constructed, which forms an exact sequence with the stress space. A matrix-matrix inf-sup condition is shown for the curl of this auxiliary space and the trilinear rotation space. First-order convergence is established for all variables in their natural norms, as well as second-order superconvergence of the displacement at the cell centers. Numerical results are presented to verify the theory.

Analysis of BDDC Preconditioners for Nonconforming Polytopal Hybrid Discretization Methods

Santiago Badia — 2026-04-09T07:00:00Z

SIAM Journal on Numerical Analysis, Volume 64, Issue 2, Page 456-484, April 2026.
Abstract. In this work, we build on the discrete trace theory developed by Badia, Droniou, and Tushar (Foundations of Computational Mathematics, in press, 2025; doi:10.1007/s10208-025-09734-6) to analyze the convergence rate of the balancing domain decomposition by constraints (BDDC) preconditioner generated from nonconforming polytopal hybrid discretizations. We prove polylogarithmic bounds on the condition number for the preconditioner that are independent of the mesh parameter and the number of subdomains and that hold on polytopal meshes. The analysis relies on the continuity of a face truncation operator, which we establish in the fully discrete polytopal setting. To validate the theory, we present numerical experiments that confirm the truncation estimate and condition number bounds. In particular, we conduct weak scalability tests for second-order elliptic problems discretized using discontinuous skeletal methods, specifically hybridizable discontinuous Galerkin and hybrid high-order methods. We also demonstrate the robustness of the preconditioner for piecewise discontinuous coefficients with large jumps.

Wasserstein Convergence Rates for Stochastic Particle Approximation of Boltzmann Models

Giacomo Borghi — 2026-04-10T07:00:00Z

SIAM Journal on Numerical Analysis, Volume 64, Issue 2, Page 485-509, April 2026.
Abstract. We establish quantitative convergence rates for stochastic particle approximation based on Nanbu-type Monte Carlo schemes applied to a broad class of collisional kinetic models. Using coupling techniques and stability estimates in the Wasserstein-1 (Kantorovich–Rubinstein) metric, we derive sharp error bounds that reflect the nonlinear interaction structure of the models. Our framework includes classical Nanbu Monte Carlo method and more recent developments as Time Relaxed Monte Carlo methods. The results bridge the gap between probabilistic particle approximations and deterministic numerical error analysis, and provide a unified perspective for the convergence theory of Monte Carlo methods for Boltzmann-type equations. As a by-product, we also obtain existence and uniqueness of solutions to a large class of binary collision models.

Complex Scaling for the Helmholtz Equation with Dirichlet Boundary Conditions in a Perturbed Half-Space

Charles Epstein — 2026-04-30T07:00:00Z

SIAM Journal on Numerical Analysis, Volume 64, Issue 2, Page 510-536, April 2026.
Abstract. We present a new analysis of complex scaling applied to the classical double layer potential for the solution of the Helmholtz equation with Dirichlet boundary conditions in compactly perturbed half-spaces in two and three dimensions. The kernel for the double layer potential is the normal derivative of the free-space Green’s function, which has a well-known analytic continuation into the complex plane as a function of both target and source locations. Here, we prove that—when the incident data are analytic and satisfy a precise asymptotic estimate—the solution to the boundary integral equation itself admits an analytic continuation into specific regions of the complex plane and satisfies a related asymptotic estimate (this class of data includes both plane waves and fields induced by point sources). We then show that, with a carefully chosen contour deformation, the oscillatory integrals are converted to exponentially decaying integrals, effectively reducing the infinite domain to a domain of finite size. No other modifications of the domain or the governing equations are introduced. We illustrate the performance of the scheme with two and three dimensional examples and discuss its extension to other boundary conditions and open waveguides.

A Deformation-Based Framework for Learning Solution Mappings of PDEs Defined on Varying Domains

Shanshan Xiao — 2026-04-30T07:00:00Z

SIAM Journal on Numerical Analysis, Volume 64, Issue 2, Page 537-564, April 2026.
Abstract. In this work, we establish a deformation-based framework for learning solution mappings of PDEs defined on varying domains. The union of functions defined on varying domains can be identified as a metric space according to the deformation; then the solution mapping is regarded as a continuous metric-to-metric mapping, and subsequently can be represented by another continuous metric-to-Banach mapping using two different strategies, referred to as the D2D subframework and the D2E subframework, respectively. We point out that such a metric-to-Banach mapping can be learned by neural networks, and hence the solution mapping is accordingly learned. With this framework, a rigorous convergence analysis is built for the problem of learning solution mappings of PDEs on varying domains. As the theoretical framework holds based on several pivotal assumptions which need to be verified for a given specific problem, we study the star domains as a typical example, and other situations could be similarly verified. There are four important features of this framework: (1) The domains under consideration are not required to be diffeomorphic, and therefore a wide range of regions can be covered by one model, provided they are homeomorphic. (2) The deformation mapping is unnecessary to be continuous, and thus it can be flexibly established via combining a primary identity mapping and a local deformation mapping. This capability facilitates the resolution of large systems where only local parts of the geometry undergo change. (3) In fact, the recent methods (Geo-FNO, DIMON, etc.) belong to the D2D subframework. We point out that the D2D subframework introduces regularity issues, whereas the proposed D2E subframework remains free from such problems. From a comprehensive perspective, the D2E subframework is better. (4) If a linearity-preserving neural operator such as MIONet is adopted, this framework still preserves the linearity of the surrogate solution mapping on its source term for linear PDEs, and thus it can be applied to the hybrid iterative method. We finally present several numerical experiments to validate our theoretical results.

A Taylor–Hood Finite Element Method for the Surface Stokes Problem Without Penalization

Alan Demlow — 2026-04-30T07:00:00Z

SIAM Journal on Numerical Analysis, Volume 64, Issue 2, Page 565-600, April 2026.
Abstract. Finite element approximation of the velocity-pressure formulation of the surface Stokes equations is challenging because it is typically not possible to enforce both tangentiality and [math] conformity of the velocity field. Most previous works concerning finite element methods (FEMs) for these equations thus have weakly enforced one of these two constraints by penalization or a Lagrange multiplier formulation. Recently in [] the authors constructed a surface Stokes FEM based on the MINI element which is tangentially conforming and [math] nonconforming, but possesses sufficient weak continuity properties to circumvent the need for penalization. The key to this method is construction of velocity degrees of freedom lying on element edges and vertices using an auxiliary Piola transform. In this work we extend this methodology to construct Taylor–Hood surface FEMs. The resulting method is shown to achieve optimal-order convergence when the edge degrees of freedom for the velocity space are placed at Gauss–Lobatto nodes. Numerical experiments confirm that this nonstandard placement of nodes is necessary to achieve optimal convergence orders.