Peer Revenue Fully funded

Papers

Fresh from arXiv, typeset properly. Highlighted words trade on the market.

Mathematics Numerical Analysis

math.NA Watch Numerical Analysis (sign in to use)

Friday 14 August

Numerical Analysis

Mixture of experts surrogate model for the homogenization of open-porous materials

For open-porous materials, incorporating their microstructural properties into mechanical simulations poses a significant challenge for accurately capturing elastic deformation. To deal with this difficulty, multiscale methods are a common tool to couple characteristics of the microstructure of the considered material with the macroscopic material behavior. However, when desiring a high accuracy, these multiscale computations can be computationally very expensive due to the large number of microscopic problems which need to be solved in each compute step. Here, surrogate models that learn the mechanical response of the underlying constitutive model can significantly reduce the computational cost of multiscale approaches. In previous work by some of the authors, beam frame models have been used to model the microstructure of open-porous materials which have been combined with neural network-based surrogate models to approximate the material behavior of a given RVE (repesentative volume element). In this work, we extend our previous study by training a more complex neural network model to predict the mechanical behavior of several RVEs, differing in their maximum pore size and pore-size distribution. Concretely, we focus on mixture of expert (MoE) models and compare different MoE architectures as well as their performance across different RVEs. This novel approach reduces the computational cost of simulating multiple RVEs as the MoE model does not require additional training when new RVEs are considered.
Read the full abstract

All authors Axel Klawonn, Martin Lanser, Lucas Mager, Ameya Rege, Janine Weber-Hamacher

Numerical Analysis

Comment on arXiv:2106.08363v3 [math.NA], E. Abreu, A. Espirito Santo, W. Lambert, and J. Perez, Convergence of a Lagrangian–Eulerian scheme by a weak asymptotic analysis for one-dimensional hyperbolic problems

This Comment concerns arXiv:2106.08363v3 [math.NA] by E. Abreu, A. Espirito Santo, W. Lambert and J. Perez, published in Numer. Methods Partial Differential Equations 39 (2023) 2400-2443. That article builds its scheme on space-time control volumes whose lateral boundaries, called "no-flow curves", solve dsigma/dt = H(u)/u and are presented as new. We show, equation by equation, that this object is the space-time integral curve of the locally conservative Eulerian-Lagrangian method of Douglas, Pereira and Yeh [Comput. Geosci. 4 (2000) 1-40], that the no-flow region is the tube those curves bound, and that its scalar forward-tracked form appears in Mancuso, Pereira and de Souza [TEMA 8 (2007) 269-276, 277-286]. At issue is not moving a control volume, which is generic, but which curve moves it: H'(u) and H(u)/u coincide identically only for a linear flux, and only the latter makes lateral mass flux vanish. The commented article itself calls the 2000 paper the first to introduce space-time local conservation, with an integral tube bounded by integral curves, while its abstract calls the same object introduced by the authors. A full-text corpus documents the changing terminology and attribution. Papers published in 2025 and 2026 use the same construction while citing later work but not DPY. Another 2026 paper gives mixed attribution: it calls the no-flow curve an "extension" of the DPY integral curve, although its own zero-flux definition and ratio ODE show identity, and it states that diffusion and dispersion do not modify the defining vector field. It extends the model, scheme and analysis, not the continuous curve. The later discrete contributions are not challenged.
Read the full abstract

Numerical Analysis

Nodal discontinuous Galerkin methods for non-ideal equations of state: pressure equilibrium preservation and entropy correction

Structure-preserving discontinuous Galerkin (DG) methods typically improve the robustness of high order simulations of real fluids. In addition to conservation, key structures include the preservation of pressure equilibrium and satisfaction of at least one entropy inequality. In this work, we investigate conservative discretizations using exactly pressure equilibrium conserving (EPEC) and approximately pressure equilibrium conserving (APEC) flux differencing DG formulations, as well as entropy stable formulations through the use of minimally dissipative corrections for non-ideal equations of state (EOS). We introduce an analysis of EPEC schemes and a new procedure for designing such fluxes based on a generalization of Tadmor's shuffle condition. We also analyze APEC DG schemes and show that the incorporation of dissipative interface penalization terms does not significantly increase pressure equilibrium errors, especially at higher orders of approximation. Finally, we observe that when combined with APEC flux differencing formulations, entropy correction improves robustness for under-resolved solutions and long-time simulations.
Read the full abstract

All authors Jesse CHan, Hendrik Ranocha, Raymond Park, Joshua Lampert, Eric Ching, Ayaboe Edoh

Numerical Analysis

Hybrid Dealiasing and Implicit Packing for Real Convolutions

Hybrid dealiasing is an FFT-based method for computing linear convolutions of complex-valued data that reduces the cost of dealiasing by performing zero padding implicitly. We develop two new algorithms that extend hybrid dealiasing to real-valued convolutions. The first algorithm exploits conjugate symmetries in the transformed data and computes each residue contribution directly. The second algorithm employs complex-valued hybrid dealiasing via a new implicit packing technique, which packs real-valued data into complex-valued arrays and performs all unpacking and packing operations implicitly in the transformed domain. Multidimensional versions of both algorithms are obtained by recursive decomposition into lower-dimensional convolutions. Based on these algorithms, efficient routines to compute convolutions of equal-length input arrays are implemented in the open-source FFTW++ library. These routines outperform the standard method of explicit zero padding in one, two, and three dimensions. In many cases, implicit packing without hybrid dealiasing is nearly as fast as the general hybrid dealiasing algorithm, particularly when multithreading.
Read the full abstract

Numerical Analysis

Geometry-Conforming Finite Element Methods for Interface Problems on Fitted and Unfitted Meshes

We develop an arbitrary-degree geometry-conforming finite element (GC-FE) framework for two-dimensional elliptic boundary value and interface problems on curved domains. Using the Frenet–Serret transformation, curved-boundary and interface-fitted segments are represented exactly, while polynomials in Frenet coordinates generate generally nonpolynomial local shape functions in physical coordinates. For interface-unfitted meshes, GC-FE spaces on curved-boundary elements are coupled with geometry-conforming immersed finite element (GC-IFE) spaces on interface-cut elements, with standard polynomial spaces used elsewhere. We establish optimal approximation, inverse, and trace estimates for the GC-FE spaces. For fitted meshes, we prove well-posedness and optimal error estimates in energy and $L^2$ norms for a symmetric interior penalty discontinuous Galerkin discretization. By retaining the prescribed curves exactly, the method avoids the geometric variational crime associated with curved-geometry approximation and requires no corresponding geometric consistency estimates. Numerical experiments confirm the predicted rates, show global accuracy comparable to nodal isoparametric finite elements and smaller true-interface trace errors in the reported tests, and demonstrate the coupled GC-FE-GC-IFE method on interface-unfitted meshes.
Read the full abstract

All authors Yuanhui Lin, Tao Lin, Xu Zhang, Minfu Feng

Numerical Analysis · Cross-listed from physics.optics

Full-field and Bloch-periodic-factor discretizations: Accuracy and phantom modes

One can distinguish four major formulations of the Bloch-mode problem in linear periodic media. First, there are two common complementary ways of slicing the Bloch variety (Bloch wave vectors paired with the corresponding frequencies): frequency vs. wave vector or, conversely, wave vector vs. frequency. This paper deals exclusively with the latter option, directly applicable to lossy media, evanescent waves, and complex band structures. A separate crossroads is the full field (FF) vs. the lattice-periodic Bloch factor (PF) formulations. These are fully equivalent on the continuous level, but common discretization techniques may break this equivalence not just approximately but qualitatively. In the FF problem, the primary unknown eigenvalue is the Bloch phase factor (not the wavenumber), which enters only through opposite-boundary coupling. The PF formulation, on the other hand, injects the Bloch wavenumber into the differential operator and leads to a volume quadratic pencil. Consequently, standard PF discretizations of the type considered here, in contrast with the corresponding FF discretizations, violate the reciprocal-lattice (Brillouin-zone-shift) covariance – which, as the theory and numerical examples in the paper show, may lead to inaccurate or even nonphysical modes. From the mathematical perspective, the paper highlights the role of structure-preserving algorithms. The practical importance of the results and recommendations stems from the wide adoption of PF formulations in optics, photonics, and beyond – such as the computation of propagating and evanescent Bloch modes in photonic and acoustic structures, topological photonics, effective-mass and topological insulator band theories.
Read the full abstract

Numerical Analysis

Reduced Order Modeling and Applications to Inverse State Estimation

Solving parametric partial differential equations (PDEs) repeatedly for many parameter values,as needed in optimal control, uncertainty quantification, and inverse problems, is often prohibitive with classical discretizations. This chapter introduces Reduced Order Modeling (ROM), which builds compressed yet accurate representations of parametric solution sets for fast online evaluation. We formulate the approximation of parametric PDEs as a supervised learning task and review linear approximation, very effective for elliptic and parabolic problems, before turning to nonlinear methods needed when solutions exhibit discontinuities or steep gradients. We then show how to use reduced-order models to efficiently recover the state of a physical system from limited measurements, which is an inverse problem known as state estimation. Throughout, we favor results with theoretical guarantees while also discussing practical methods, and provide a companion Jupyter notebook implementing the main algorithms.
Read the full abstract

Numerical Analysis · Cross-listed from math.AP

Analytic inverse problems with finitely many random measurements

While infinite-dimensional inverse problems are traditionally analyzed assuming continuous data, practical applications rely on finitely many discrete measurements. Recent deterministic approaches establish that unknowns belonging to a $d$-dimensional model class can be stably recovered from finitely many measurements. However, for severely ill-posed problems, such as the Calderón problem and inverse scattering, the known constructions may require a number of measurements that is exponential in $d$. We show that random sampling reduces this count dramatically if one asks only for exact identifiability. By exploiting the analytic geometry of the forward maps, we prove that, whenever the infinite-data problem is injective on the model class, $2d+1$ random scalar measurements determine the unknown uniquely, almost surely. Applications are given to the Calderón problem, with both infinite- and finite-dimensional boundary sampling, and to inverse medium scattering from randomly sampled far-field values.
Read the full abstract

All authors Giovanni S. Alberti, Damiano Poletti, Simone Sanna, Matteo Santacesaria

Numerical Analysis

Coupling-Aware Vanka Smoothing for Multigrid Preconditioning of the Implicit Immersed Boundary Equations

The immersed boundary (IB) method models fluid–structure interaction using the natural Lagrangian and Eulerian formulations of structural mechanics and fluid dynamics, respectively, but explicit time discretization of the IB force imposes a stiffness-dependent upper bound on the time-step size. Treating these forces implicitly removes this restriction, but the resulting coupled linear systems become increasingly difficult to solve as the structural stiffness increases. Algebraically eliminating the Lagrangian degrees of freedom yields a reduced Eulerian velocity–pressure IB system. Here we introduce a coupling-aware Vanka (CAV) smoothing strategy to enable effective multigrid preconditioning of this system. CAV patches are built as unions of standard pressure-centered Vanka patches, with the graph of the Eulerian elasticity matrix determining which patches are combined. Under grid refinement, CAV patch sizes remain bounded, and the computational cost of each multigrid cycle thereby grows linearly with the number of Eulerian degrees of freedom. Tests using target-point, membrane, and beam force laws show that CAV-preconditioned FGMRES reduces the relative residual by ten orders of magnitude in $9--15$ iterations, with little growth under grid refinement. In a nonlinear benchmark modeling flow past a flexible fiber, the average number of FGMRES iterations per Newton solve increases only from $8.6$ to $9.5$ as the Eulerian grid is refined from $32\times32$ to $256\times256$ cells. To our knowledge, CAV provides the first robust multigrid strategy for time-dependent implicit IB formulations.
Read the full abstract

Numerical Analysis

Analysis of a Conforming Finite Element Method for Second-Harmonic Generation Scattering

In the frequency domain, nonlinear acoustic and electromagnetic wave propagation can, in certain regimes, be modeled by second-harmonic generation systems consisting of two coupled nonlinear Helmholtz equations. We analyze a conforming finite element approximation of a scalar second-harmonic generation scattering problem truncated using an exact Dirichlet-to-Neumann boundary condition. The discretization uses continuous piecewise polynomial finite elements of degree $p$. Under small-data conditions involving the incident-field and nonlinear susceptibilities, we prove existence and uniqueness of the continuous and discrete nonlinear solutions in a prescribed small ball. In the same regime, we derive quasi-optimal a priori $H^1$-error estimates by a nonlinear Céa-type argument combining linear Galerkin quasi-optimality with small-data stability estimates for the coupled nonlinearities. We also analyze the convergence of the fixed-point iteration used to compute the discrete nonlinear solution and quantify the combined effects of finite element approximation and nonlinear iteration error. A manufactured-solution experiment illustrates the predicted finite element convergence rates, while additional PML-based scattering computations demonstrate the behavior of the nonlinear fixed-point solver in two and three dimensions.
Read the full abstract

Numerical Analysis

TOGEARI: Interaction-Space Preconditioning for Condensed Finite-Element Systems with IPC Contact

TOGEARI constructs a compressed contact correction for condensed finite-element equations with incremental potential contact (IPC). A core-conditioned interaction matrix ranks contact combinations by their mechanical responses; a raw-factor proxy acquires only the retained responses. For a symmetric positive-definite core, the reference analysis gives the exact spectrum, scaled inverse error and optimal worst omitted interaction. The selected responses are maintained through signed contact changes using supported operator images and a small coarse matrix. A complementary Schur spectrum and perturbation bounds connect this construction to the full current linear equation, damped local Newton steps and compatible episode reuse. On a 325,260-coordinate finger/cup problem, eight selected directions reduce Krylov work from eleven columns to five, matching full contact. Two reversed-order load studies give matching anchored displacements and accepted normal reactions, with 28.8-29.3 percent less warm time and 5.1-5.3 percent less complete time than per-step direct refactorization. A completed twelve-episode run saves 18.8 percent of total time against direct and 2.8 percent against a held factor. Full response caching has comparable time with 7.5 times the response storage. A compressed repeat terminates at an energy line search; separate diagnostics identify constitutive cancellation and terminal-backtracking sensitivity. Larger-contact, held-out-load and frictional comparisons show where compression or preparation loses its advantage. The complete Newton equation remains the residual operator throughout; the convergence guarantees retain their stated reference and local hypotheses.
Read the full abstract

All authors Yanlin Liu, Chao Huang, Kaixiang Yao, Yao Shen

Numerical Analysis

Adaptive multigrid for high-order discontinuous Galerkin methods based on the full approximation scheme

We propose an adaptive multigrid (MG) method for discontinuous Galerkin formulations of elliptic problems using Brandt's full approximation scheme (FAS). Unlike common approaches, this method achieves local $hp$-refinement of hexahedral meshes without the need for hanging nodes. The core component of the FAS-MG method is an overlapping Schwarz smoother, which is optionally accelerated by a Krylov method. This smoother is designed for unstructured curvilinear meshes but maintains a tensor-product structure for fast diagonalization. Numerical experiments demonstrate the exceptional efficiency of the FAS-MG method. Dedicated studies confirm its robustness against high aspect ratios, element deformation, and irregular mesh topology. We also verify its capability for dynamic parallel mesh adaptation using the wave-front benchmark of Červený, Dobrev, and Kolev (SIAM J. Sci. Comp. 41, 2019). Finally, we present preliminary results of extending the method to incompressible Navier-Stokes problems.
Read the full abstract

Numerical Analysis

Unconditionally stable and energy conserving discretization of the dynamic von Kármán equations

A fully discrete approximation of the dynamic von Kármán equations combines nonconforming Morley finite element methods for spatial discretization with an energy conserving modified unconditionally stable Newmark second- order time-stepping scheme. Brouwer's fixed-point theorem establishes existence of a solution to the fully discrete scheme and further uniqueness and stability estimates follow for small loads. Optimal order a priori error estimates in the piecewise energy norm with quadratic convergence in time are derived for the fully discrete scheme. The results of the numerical experiments validate the theoretical error bounds.
Read the full abstract

All authors Carsten Carstensen, Neela Nataraj, Ricardo Ruiz-Baier, Aamir Yousuf

Numerical Analysis · Cross-listed from math.MG

A counterexample to Kusner's conjecture on equilateral sets

We disprove Kusner's 1983 conjecture that every equilateral set in $\ell_p^n$ with $2<p<\infty$ has at most $n+1$ points: there exist $58$ points in $\mathbb{R}^{56}$ whose pairwise $\ell_5$ distances are all equal, so the maximum equilateral-set size satisfies $e(\ell_5^{56})\ge58>57$. This is the first equilateral set of more than $n+1$ points in $\ell_p^n$ for any finite $p\ge2$. The construction persists on an open interval of exponents around $5$; since Ge, Xu and Zhou recently proved the conjecture for $2\le p\le4$, the infimum of exponents at which it fails lies in $[4,5)$. The configuration is the unique solution of an explicit polynomial system with rational coefficients in a rational box, established in exact arithmetic.
Read the full abstract

Numerical Analysis

An $h$-adaptive Tetrahedral Spectral Element Method with Applications to Kohn-Sham Density Functional Theory

High-order $h$-adaptive spectral element methods on tetrahedral meshes provide an effective framework for resolving localized singularities and multiscale structures in complex three-dimensional geometries. However, their development is often hindered by difficulties in maintaining $C^0$ continuity across refinement interfaces and efficiently transferring solutions between adaptive meshes. Such limitations are particularly relevant in demanding applications such as all-electron Kohn-Sham density functional theory, which place stringent requirements on the accurate resolution of both nuclear singularities and multiple physical scales. In this paper, we present an efficient $h$-adaptive tetrahedral spectral element framework. To address the continuity challenge, we develop an adaptive strategy that combines element orientation alignment with geometric red-green refinement, thereby eliminating the need for algebraic hanging-node constraints while preserving inter-element continuity. Furthermore, an efficient topology-based point-location algorithm is introduced to accelerate interpolation between adaptive meshes. Numerical experiments on Poisson and Laplacian eigenvalue problems confirm the spectral convergence of the proposed method. Applications to all-electron Kohn-Sham equations further demonstrate its capability to accurately resolve nuclear singularities. Moreover, parallel performance studies exhibit excellent scalability, with matrix assembly and adaptivity modules generally achieving speedups above 15 times and the proposed interpolation algorithm attaining speedups ranging from 25 to 35 on 64-core configurations compared to the single-core performance. These results indicate that the proposed framework provides an accurate, robust, and efficient solution for large-scale, high-resolution simulations.
Read the full abstract

Numerical Analysis

Nonparametric Schrödinger Bridge Time Series Generator: Algorithm, Convergence Analysis and Applications

We conduct a convergence analysis for the Schrödinger Bridge Time Series (SBTS) data generator. Starting from a regularized formulation in which the data ensemble is mixed with a standard multivariate Gaussian distribution with a prescribed probability, we prove that the Euler-Maruyama discretization converges to the mixed target distribution with half-order convergence rate, provided that the ensemble size and kernel bandwidth are chosen appropriately. We further show that the regularized distribution converges to the original target distribution as the mixing probability tends to zero. The analysis simultaneously accounts for the ensemble approximation error, kernel approximation error, and time-discretization error, and therefore provides a full distributional convergence result for the SBTS generator. Empirically, we further examine the flexibility of the method by replacing the Wiener reference measure with the path measure induced by a more general SDE. The numerical experiments show that the schemes based on both the original Wiener reference measure and the SDE-induced reference measure achieve comparable performance, demonstrating the robustness and stability of the SBTS framework.
Read the full abstract

Numerical Analysis · Cross-listed from stat.ME

Semi-supervised Concordance Learning for Optimal Individual Treatment Regimes

Finding the optimal individualized treatment rule that maps individual characteristics or contextual information to treatment assignments has been extensively investigated in existing literature, with widespread practical applications. This paper considers the estimation of optimal treatment regimes within a semi-supervised data framework (exemplified by electronic medical record data). In such settings, only a tiny proportion of observations have observed outcome labels, owing to high labeling costs, time limitations, data privacy concerns, and other constraints, while covariates and treatment assignments are available for all study subjects. We develop a semi-parametric inference method for optimal treatment regimes, which leverages outcome- unlabeled samples with complete covariate and treatment information to enhance estimation efficiency. The proposed estimation framework consists of two key steps: first, flexible nonparametric imputation via single-index kernel smoothing; second, subsequent estimation of the optimal treatment regime based on concordance-assisted learning. We establish the consistency and asymptotic normality of our proposed estimators. Numerical simulation studies demonstrate that our method achieves higher efficiency and stronger robustness relative to fully supervised estimators under finite-sample settings. We further validate the practical value of our proposed framework using the MIMIC-III and ACTG175 datasets.
Read the full abstract

Numerical Analysis

PDE Realization and Structure-Aware Solvers for a Compact Two-Stage Fourth-Order IMEX Method

This paper develops a PDE realization and solver framework for a compact two-stage fourth-order two-derivative implicit–explicit time discretization for stiff split problems. The key difficulty is a consistency–cost coupling: full-field temporal differentiation is required to preserve the mixed explicit–implicit interactions of the time integrator, but the same interactions widen the implicit stage operators and can make each solve substantially more expensive. We show that a second-order ADER/Cauchy–Kowalevski local evolution is sufficient for the fourth-order outer composition, establish smooth reconstruction consistency and a fully discrete error balance, and derive the leading Lie-bracket defect produced by self differentiation of the split fields. An inexact-stage analysis gives asymmetric midpoint and endpoint residual tolerances that preserve fourth-order accuracy. For the widened stages, the complete mixed action is retained matrix-free while only dominant stiff physics is approximated in the inverse. This yields quadratic and shifted preconditioners, Fourier and multilevel realizations, a semilinear reaction–diffusion reduction, and source-local elimination for relaxation systems. Exact quadratic cancellation, a diffusion-dominated Fourier estimate, and an $\varepsilon/\Delta x$ bound for Jin–Xin relaxation explain the main solver mechanisms. Numerical ablations verify the consistency and tolerance results, while a two-dimensional Brusselator study on two grids shows favorable error-versus-wall-time behavior over a useful accuracy range. A classical stiff-front benchmark also identifies the separate spatial shock–source limitation. The results support a regime-dependent efficiency claim rather than a universal speedup.
Read the full abstract

Numerical Analysis

Vanishing Clutter: Fast and Accurate Shape Imaging via Active Cloaking

Imaging a target sample in near-field scanning optical microscopy (NSOM) is fundamentally limited by measurement artifacts and data contamination from multiple probe-sample scattering. The probe, essential for subwavelength resolution, inherently perturbs the local field, degrading the signal-to-noise ratio and rendering the inverse problem for quantitative shape reconstruction highly ill-posed. We first establish the well-posedness of the corresponding forward model, providing a rigorous foundation for subsequent imaging. We then repurpose active cloaking–conventionally the antagonist of imaging–as an enabling mechanism to eliminate probe-induced interference. Rather than directly reconstructing the sample from corrupted data, we actively cloak the probe by formulating an optimal control problem and prove the existence and stability of its minimizers. Leveraging the theory of localized anomalous resonance in layered plasmonic structures, we derive an exact closed-form minimizer, thereby circumventing the computationally prohibitive iterative solution of the optimal control problem. The resulting cloaking-driven interference removal yields a virtually probe-free measurement environment, enabling fast, artifact-free shape reconstruction. Extensive numerical experiments demonstrate accurate shape reconstruction and dramatic acceleration–often by orders of magnitude–over conventional iterative methods applied directly to probe–contaminated data without cloaking–based preprocessing, validating the robustness and transformative potential of the proposed approach for high-fidelity subwavelength imaging.
Read the full abstract

Thursday 13 August

Numerical Analysis

An Efficient Solver for Finite Element-based Constrained Transport in 3D Magnetohydrodynamics Applied to Magnetic Confinement Fusion

We present an efficient solver framework for the stiff magnetic wave coupling arising in resistive magnetohydrodynamics (MHD) on realistic tokamak geometries. The approach builds on an implicit-implicit (IMIM) time-splitting that separates fast magnetic waves and anisotropic heat transport from slower acoustic dynamics while retaining full coupling (Krzysik et al. 2026). Within this formulation, the magnetic wave subsystem appears as an anisotropic curl-curl operator, enabling the use of scalable auxiliary-space Maxwell (AMS) multigrid solvers. To exploit this structure at the discrete level, we employ curl-conforming finite element spaces for the magnetic field and design the velocity space to preserve the curl-curl structure induced by the Lorentz-force coupling. The resulting compatible discretization preserves the discrete magnetic divergence constraint while producing linear systems directly amenable to efficient AMS-based solvers. We demonstrate solver efficiency as well as the accuracy and stability of the resulting structure-preserving discretization on fully nonlinear three-dimensional tokamak test cases.
Read the full abstract

All authors Golo A. Wimmer, Konstantin Lipnikov, Ben S. Southworth, Xian-Zhu Tang

Numerical Analysis

Approximating matrix functions by block Krylov methods with randomized vectors

The need to evaluate expressions of the form $f(A)\mathbf{b}$, where $A$ is a square matrix, $f$ is a function, and $\mathbf{b}$ is a vector, arises in several areas of applied mathematics. When the matrix $A$ is very large, it is usually not attractive to evaluate $f(A)$. Instead, $f(A)\mathbf{b}$ often is approximated by computing an estimate in a Krylov subspace that depends on $A$ and $\mathbf{b}$, and only requires that $f$ be evaluated at a small matrix. This paper explores the application of several variants of randomized block Krylov methods to the approximation of $f(A)\mathbf{b}$. Computed examples suggest that block Krylov methods with an initial block vector that contains $\mathbf{b}$ as well as a few randomly generated vectors may require less computing time and reduce the number of Krylov steps than standard Krylov methods.
Read the full abstract

All authors Josh Kane, Lucas Onisk, Lothar Reichel, Giuseppe Rodriguez

Numerical Analysis

An adaptive superconvergent mixed finite element method for exactly symmetric linear elasticity based on local residual minimization

We introduce an a posteriori error estimator for exactly symmetric mixed finite element discretizations of linear elasticity, based on recasting a Stenberg-type postprocessing scheme as a local residual minimization problem in a discrete dual norm. This construction yields, as a dual variable and at no additional computational cost, a Riesz representative of the associated local residual, from which we build an a posteriori error indicator. We establish a reliability estimate, with the dependence on the Lamé parameters tracked explicitly, and a local efficiency estimate, without auxiliary bubble functions, in the standard compressible regime, as well as an alternative reliability estimate, based on a robust stability property and an Oswald averaging operator, with a constant that remains bounded as $\lambda\to\infty$; the local efficiency estimate holds, with the same bounded behavior, uniformly in both regimes. Notably, a single indicator and a single comparison norm serve both regimes, in contrast with existing hypercircle-based estimators, which require a distinct construction for the incompressible limit. Numerical examples, including a benchmark with a known singular solution and one without an analytical solution, validate the theoretical findings.
Read the full abstract

Numerical Analysis

A new approach to the reaction-diffusion systems modelling

We consider a new methodology for modelling the reaction-diffusion systems based on systems of ordinary differential equations. In contrary to the specialised numerical methods like straight line method, this new methodology is positioned as a pure alternative at the model level for partial differential equations. In its description, the new method is largely similar to the finite volume method, but unlike the latter, it uses statistical simplifications and principles of geometric probability to describe diffusion. The main objectives of this approach are to simplify the qualitative analysis of reaction-diffusion systems and to improve the efficiency of numerical model implementation. The first objective is successfully addressed, as it becomes possible to use the apparatus of classical dynamical systems theory for a qualitative analysis of model dynamics based on the systems of ordinary differential equations. The second objective is only partially addressed, as the gain in efficiency while maintaining acceptable accuracy for numerical implementation will be significant only for certain simple initial distribution of molecules and for specific diffusion coefficients. Furthermore, to formulate criteria for practical applicability, we separately evaluate the modelling error using this new methodology.
Read the full abstract

Numerical Analysis

Stochastic resistive Hall–MHD with current fluctuations: martingale weak solutions via a convergent structure-preserving finite element method

We study the stochastic resistive Hall–magnetohydrodynamic (Hall–MHD) system on bounded convex polyhedral domains, subject to nonlinear perfectly conducting boundary conditions. The system is driven by multiplicative Gaussian forcing in the momentum equation together with structured curl-type noise in the induction equation arising from stochastic perturbations of the resistive and Hall parts in the generalised Ohm law. We develop a fully discrete, linearly implicit, structure-preserving mixed finite element method based on a compatible discrete de Rham complex, which preserves the magnetic Gauss law exactly. We prove subsequential convergence of the discrete solutions to a finite-energy weak martingale solution, thereby obtaining a constructive existence result for the stochastic Hall–MHD system with curl-type noise. Numerical experiments illustrate the stochastic dynamics, exact divergence preservation, and convergence of the method.
Read the full abstract

Numerical Analysis

Performance Evaluation of an Adaptive Quadrature and a Double Exponential Formula Using Arbitrary-Precision Floating-Point Arithmetic

Using arbitrary-precision arithmetic provided by the GNU Multiple Precision Floating-Point Reliable Library, we implement AQE11D—that is, Ninomiya's adaptive 9-point Newton–Cotes rule extended with a sequence of higher-order rules—and Takahasi and Moris' double exponential (DE) formula. We evaluate them for Kahaner's 21 test problems. For both absolute tolerances $10^{-50}$ and $10^{-100}$, AQE11D attains target accuracy on all 21 problems; however, for strong endpoint singularity such as $1/\sqrt{x}$, it requires about $5.4\times10^{7}$ function evaluations at $10^{-100}$, roughly $7\times10^{4}$ times as many as the DE formula. The formula converges on 18 problems at both tolerances, demonstrating its strength against endpoint singularities but also its failure, as it stands, on problems with a singularity inside the integration interval.
Read the full abstract

Numerical Analysis

Analysis of Error Propagation in Autoencoder-Based Reduced-Order Neural Ordinary Differential Equations

Neural ODE reduced-order models often achieve comparable local prediction accuracy, yet their long-horizon extrapolation behavior can differ substantially. To analyze this discrepancy, we develop a path-integral identity that separates local discrepancy injection from amplification in the learned latent dynamics. The associated multi-step Jacobian norms quantify transport sensitivity and distinguish different propagation regimes. Experiments on the Burgers and Gray–Scott systems exhibit two distinct patterns of error evolution. In Burgers systems, prediction errors remain bounded and are primarily associated with persistent local discrepancies. In contrast, Gray–Scott systems exhibit pronounced amplification during extrapolation, where Jacobian norms serve as sensitivity diagnostics rather than direct indicators of physical prediction accuracy.
Read the full abstract

Numerical Analysis · Cross-listed from math.OC

Adaptive Schauder Stochastic Mirror Descent in Banach Spaces

In this paper, we introduce an adaptive regularization strategy for stochastic mirror descent (SMD) to solve a class of risk functional minimization problems in infinite-dimensional Banach spaces. This regularization strategy centers on using a Schauder basis to construct a nested family of finite-dimensional subspaces, with the dimension chosen adaptively according to the sample size $n$. We then restrict each SMD subproblem to the corresponding subspace and project the stochastic gradient onto its dual space. This yields closed-form solutions to the SMD subproblems and coordinate-wise updates of the basis coefficients, enabling an implementation with low computational and storage complexity. The subspace dimension also serves as a regularization parameter that balances approximation and optimization errors. For risk functional minimization in $\mathcal{L}^p$ spaces with $1<p<\infty$, we construct Bregman distances adapted to the geometry of the underlying Banach spaces using $\max\{2,p\}$-convex functionals induced by their uniform convexity. At the non-uniformly convex $\mathcal{L}^1$ endpoint, we instead construct a locally strongly convex functional based on the entropy function. By developing a new analytical framework, we establish a convergence rate of $\mathcal O\left(n^{-\min\{\frac12,\frac1p\}}\right)$, up to logarithmic factors. In the misspecified setting, where the minimizer satisfies only weaker regularity conditions, we prove that the risk functional still converges to its minimum value. Finally, we apply the method to statistical inverse problems and illustrate its empirical performance through numerical experiments in both settings.
Read the full abstract

Numerical Analysis · Cross-listed from math-ph

Transparent Boundary Conditions for the Heat Equation on Metric Graphs

We study transparent boundary conditions (TBCs) for the time-dependent heat equation in a branching, quasi-one-dimensional domain modeled as a metric star graph. By combining the classical concept of TBCs with the theory of partial differential equations (PDEs) on networks, we derive an exact vertex condition that ensures unobstructed thermal flow across junctions. Specifically,we derive a sum rule for the diffusion coefficients that eliminates thermal backflow at the vertex.We demonstrate the validity of this analytical model numerically using a Crank-Nicolson finite-difference method, which confirms smooth, unobstructed heat propagation through the network.These results provide a practical mathematical framework for tunable control and optimization of thermal diffusion in low-dimensional structures. In particular, combining the well-known concept of the TBCs and theory for PDE on metric graphs, we propose a mathematical model providing a control tool for thermal diffusion in networks.
Read the full abstract

Numerical Analysis

A Counterexample to Robust Second-Order Convergence of the Strang Projector-Splitting Integrator

The classical Strang projector-splitting integrator is widely observed to converge with order two, whereas the error analysis that remains uniform as the smallest singular value retained in the low-rank approximation tends to zero proves only order one. We show that this gap is intrinsic under the standard assumptions. We construct $3\times3$ matrix differential equations that are $C^2$ in time and smooth in the matrix variable, with rank-two initial data that satisfy uniform boundedness, Lipschitz, tangency-defect, and regularity bounds. Nevertheless, the exact-subflow Strang method has a nonzero $h^2$ term in the local error over one periodic forcing cycle consisting of four Strang steps. Repetition of that cycle rules out a global second-order bound whose constant and stepsize threshold are independent of the retained singular values. The mechanism is a rapid rotation of the factor directions associated with the small singular value: first-order consistency is preserved, but the changing projection spaces prevent the cancellation normally expected from a symmetric Strang composition. Hence the robust first-order result cannot, under these assumptions alone, be upgraded to robust second order for the classical projector-splitting method.
Read the full abstract

Numerical Analysis

A robust and optimal method for a sixth-order elliptic singular perturbation problem

In this paper we propose a robust nonconforming finite method for a gradient-elastic Kirchhoff plate (GEKP) model with a small parameter $0<\iota<1$. While the method employs an $H^3$ cubic finite element for the discrete space, it employs an interpolation mapping into an $H^2$ nonconforming finite element space in both the lower-order terms of the bilinear form and the right-hand side. Notably, the numerical method is robust with respect to the parameter $\iota$ and uniformly convergent to order $h$ without any stabilized terms. Finally, some numerical experiments are performed to verify the theoretical findings.
Read the full abstract