Abstract. This paper develops and discusses a residual-based a posteriori error estimator for parabolic surface partial differential equations on closed stationary surfaces. The full discretization uses the surface finite element method in space and the backward Euler method in time. The proposed error indicator bounds the error quantities globally in space from above and below, and globally in time from above and locally from below. Based on the derived error indicator, a space-time adaptive algorithm is proposed. Numerical experiments illustrate and complement the theory.
Abstract. In this paper we analyze a first-order exponential wave integrator (EWI) for the nonlinear Schrödinger equation (NLSE) with a singular potential that is locally in [Formula: see text], which might be locally unbounded. A typical example is the inverse power potential such as the Coulomb potential, which is the most fundamental potential in quantum physics and chemistry. We prove that, under the assumption of [Formula: see text]-potential and [Formula: see text]-initial data, the [Formula: see text]-norm convergence of the EWI is, roughly, first-order in one dimension (1D) and two dimensions (2D), and [Formula: see text]-order in three dimensions (3D). In addition, under a stronger integrability assumption of [Formula: see text]-potential for some [Formula: see text] in 3D, the [Formula: see text]-norm convergence increases to almost [Formula: see text]-order if [Formula: see text] and becomes first-order if [Formula: see text]. In particular, our results show, to the best of our knowledge for the first time, that first-order [Formula: see text]-norm convergence can be achieved when solving the NLSE with the Coulomb potential in 3D. The key advancements are the use of discrete (in time) Strichartz estimates, which allow us to handle the loss of integrability due to the singular potential that does not belong to [Formula: see text], and the more favorable local truncation error of the EWI, which requires no (spatial) smoothness of the potential. Extensive numerical results in 1D, 2D, and 3D are reported to confirm our error estimates and to show the sharpness of our assumptions on the regularity of the singular potentials.
Abstract. In this paper we consider an evolving surface finite element method method for the advection and diffusion of a scalar quantity on a moving closed curve. The diffusion process is controlled by a forcing term that may include a rough term (specifically a stochastic noise) which in particular destroys the classical time differentiability properties of the solution. We provide a suitable variational solution concept and a fully discrete finite element method discretization. Our error analysis appropriately generalizes classical estimates to this weaker setting. We present some numerical simulations that confirm our theoretical findings.
Signature kernels, inner products of path signatures, underpin several machine learning algorithms for multivariate time series analysis. For bounded variation paths, the signature kernel was recently shown to solve a Goursat PDE. However, existing PDE solvers only use increments as input data, leading to first-order approximation errors. These approaches become computationally intractable for highly oscillatory input paths, as they have to be resolved at a fine enough scale to accurately recover their signature kernel, resulting in significant time and memory complexities. In this paper, we extend the analysis to rough paths and show, leveraging the framework of smooth rough paths, that the resulting rough signature kernels can be approximated by a novel system of PDEs whose coefficients involve higher-order iterated integrals of the input rough paths. We show that this system of PDEs admits a unique solution and establish quantitative error bounds yielding a higher-order approximation to rough signature kernels.
ParaOpt is a time parallel method based on Parareal for solving optimality systems arising in optimal control problems. The method was presented in [M. J. Gander, F. Kwok and J. Salomon, SIAM J. Sci. Comput., 42 (2020), A2773--A2802] together with a convergence analysis in the case where implicit Euler is used to discretize the differential equations governing the system dynamics. However, its convergence behavior for higher-order time discretizations has not been considered. In this paper, we use an operator norm analysis to prove that the convergence rate of ParaOpt applied to a linear-quadratic optimal control problem has the same order as the Runge-- Kutta time integration method used, provided that a few auxiliary order conditions are satisfied. We illustrate our theoretical results with numerical examples, before showing an additional test case not covered by our analysis, namely, a nonlinear optimal control problem involving a Schro"\dinger type system.
Deriving sharp a posteriori upper bounds of the Lipschitz constant of deep neural networks is crucial to formally guarantee the robustness of neural network based models. We analyze three existing upper bounds written for the l2 norm. We highlight the interest of working with the l1 and l\infty norms and we propose two novel bounds for feed-forward fully connected neural networks, which can be extended to convolutional neural networks. Several numerical tests confirm the theoretical results, help to quantify the relationship between the presented bounds, and establish the better accuracy of the new bounds. Four numerical tests are studied: one with random matrices; two where the output has an analytical closed form; and one for convolutional neural networks trained on the MNIST dataset. The last numerical test is an application to a reference convolutional neural network. One of our bounds is optimal in the sense that it is exact for the second test uniformly with respect to the depth of the neural network and it is always better than the other ones.
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.
Efficient simulation of the semiclassical Schro"\dinger equation has garnered significant attention in the numerical analysis community. While controlling the error in the unitary evolution or the wavefunction typically requires the time step size to shrink as the semiclassical parameter h decreases, it has been observed--and proved for first-and second-order Trotterization schemes--that the error in certain classes of observables admits a time step size independent of h. In this work, we explicitly characterize this class of observables and present a new, simple algebraic proof of uniform-in-h error bounds for arbitrarily high-order Trotterization schemes. Our proof relies solely on the algebraic structure of the underlying operators in both the continuous and discrete settings. Unlike previous analyses, it avoids Egorov-type theorems and bypasses heavy semiclassical machinery. To the best of our knowledge, this is the first proof of uniform-in-h observable error bounds for Trotterization in the semiclassical regime that relies only on algebraic structure, without invoking the semiclassical limit.
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 sub domains 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.
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.
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 lowerPartial 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 L2-norm. A cut-out argument allows converting this into an L1-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.
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 \scrT -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.
In this paper, we develop and numerically implement a novel approach for solving the inverse source problem of the acoustic wave equation in three dimensions. By injecting a small high-contrast droplet into the medium, we exploit the resulting wave field perturbation measured at a single external point over time. The method enables stable source reconstructions where conventional approaches fail due to ill-posedness, with potential applications in medical imaging and non-destructive testing. Key contributions include: 1. Implementation of a theoretically justified asymptotic expansion, from [33], using the eigensystem of the Newtonian operator, with error analysis for the spectral truncation. 2. Novel numerical schemes for solving the time-domain Lippmann-Schwinger equation and reconstructing the source via Riesz basis expansions and mollification-based numerical differentiations. 3. Reconstruction requiring only single-point measurements, overcoming traditional spatial data limitations. 4. 3D numerical experiments demonstrating accurate source recovery under noise (SNR of the order 1/a), with error analysis for the droplet size (of the order a) and the number of spectral modes N.
We introduce \emph{coarse scrambling}, a novel randomization for digital sequences that permutes blocks of digits in a mixed-radix representation. This construction is designed to preserve the powerful $(0,\boldsymbol{e},d)$-sequence property of the underlying points. For sufficiently smooth integrands, we prove that this method achieves the canonical $O(n^{-3+ε})$ variance decay rate, matching that of standard Owen's scrambling. Crucially, we show that its maximal gain coefficient grows only logarithmically with dimension, $O(\log d)$, thus providing theoretical robustness against the curse of dimensionality affecting scrambled Sobol' sequences. Numerical experiments validate these findings and illustrate a practical trade-off: while Owen's scrambling is superior for integrands sensitive to low-dimensional projections, coarse scrambling is competitive for functions with low effective truncation dimension.
This paper develops divergence-free mixed finite element methods for the Stokes equation. Using H(div)-conforming velocities and discontinuous pressures ensures the inf-sup condition for the velocity–pressure pair and yields pointwise divergence-free velocities. However, this choice makes the vector Laplacian difficult to discretize. Inspired by mass-conserving mixed formulations with stresses, tangential–normal continuous traceless tensor elements are introduced to discretize the vector Laplacian. An inf-sup condition for the weak div operator between the stress and velocity spaces is then proved. Two key properties characterize the scheme. First, the stress–velocity inf-sup stability gives a stable discretization of the vector Laplacian without additional stabilization, unlike discontinuous Galerkin or virtual element methods. Second, the scheme has the property that if a stress field is weakly divergence-free, then it is also strongly divergence-free. This decouples the stress and velocity errors and leads to superconvergence. As a result, optimal-order error estimates are obtained for the stress, while the velocity and pressure converge at rates higher than the approximation orders of the chosen spaces. Numerical experiments confirm the theoretical results.
We consider parabolic evolution equations with Lipschitz continuous and strongly monotone spatial operators. By introducing an additional variable, we construct an equivalent system where the operator is a Lipschitz continuous mapping from a Hilbert space Y imes X to its dual, with a Lipschitz continuous inverse. Resulting Galerkin discretizations can be solved with an inexact Uzawa type algorithm. Quasi-optimality of the Galerkin approximations is guaranteed under an infsup condition on the selected ``test"" and ``trial"" subspaces of Y and X. To circumvent the restriction imposed by this inf-sup condition, an a posteriori condition for quasi-optimality is developed that is shown to be satisfied whenever the test space is sufficiently large.
The Zarantonello fixed-point iteration is an established linearization scheme for quasilinear PDEs with strongly monotone and Lipschitz continuous nonlinearity. This paper presents a weighted least-squares minimization for the computation of the update of this scheme. The resulting formulation allows for a conforming finite element discretization of the primal and dual variable of the PDE with arbitrary polynomial degree. The least-squares functional provides a built-in a posteriori discretization error estimator in each linearization step motivating an adaptive Uzawa-type algorithm with an outer linearization loop and an inner adaptive mesh-refinement loop. We prove R-linear convergence of the linearization iterates for arbitrary initial guesses. Particular focus is on the role of the weights in the least-squares functional of the linearized problem and their influence on the robustness of the Zarantonello damping parameter. Numerical experiments illustrate the performance of the proposed algorithm.
We analyze rates of uniform convergence for a class of high-order semi-Lagrangian schemes for first-order, time-dependent partial differential equations on embedded submanifolds of Rd (including advection equations on surfaces) by extending the error analysis of Falcone and Ferretti [SIAM J. Numer. Anal., 35 (1998), pp. 909--940]. A central requirement in our analysis is a remapping operator that achieves both high approximation orders and strong stability, a combination that is challenging to obtain and of independent interest. For this task, we propose a novel mesh-free remapping operator based on \ell1 minimizing generalized polynomial reproduction, which uses only point values and requires no additional geometric information from the manifold (such as access to tangent spaces or curvature). Our framework also rigorously addresses the numerical solution of ordinary differential equations on manifolds via projection methods. We include numerical experiments that support the theoretical results and also suggest some new directions for future research.
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, H1 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.
In this paper, a theoretical framework is presented for the use of a Kansa-like method to numerically solve elliptic partial differential equations on spheres and other manifolds. The theory addresses both the stability of the method and provides error estimates for two different approximation methods. A Kansa-like matrix is obtained by replacing the test point set X, used in the traditional Kansa method, by a larger set Y, which is a norming set for the underlying trial space. This gives rise to a rectangular matrix. In addition, if a basis of Lagrange (or local Lagrange) functions is used for the trial space, then it is shown that the stability of the matrix is comparable to the stability of the elliptic operator acting on the trial space. Finally, two different types of error estimates are given. Discrete least squares estimates of very high accuracy are obtained for solutions that are sufficiently smooth. The second method, giving similar error estimates, uses a rank revealing factorization to create a “thinning algorithm” that reduces #Y to #X. In practice, this algorithm doesn't need Y to be a norming set.
