Papers

In this article, we develop high-order, entropy-dissipative schemes for a hyperbolic system of conservation laws describing the first-order dynamics of two-layer thin-film flows of immiscible fluids with perfectly soluble solute particles. The key aspect is to construct entropy-conservative fluxes in the sense of Tadmor. Here, we consider and compare two different construction processes: an exact formulation and a commonly used simplified approximation in the flux evaluation. While the approximate approach is simpler at the level of construction and is frequently employed in practice (also in other contexts), it is less structured, leading to increased computational cost in the implementation compared to the exact formulation, which has been much more complicated to derive. We employ these fluxes within both an entropy-dissipative finite-difference framework and an entropy-dissipative discontinuous Galerkin spectral element method. Through numerical experiments using two distinct spatial discretizations, we investigate the performance and robustness of the proposed methods and demonstrate that the choice of numerical flux has a noticeable impact on the efficiency of the resulting solvers.

In [H. Lavenant, J. Luckhardt, G. Mordant, B. Schmitzer, L. Tamanini, The Riemannian geometry of Sinkhorn divergences. Ann. Inst. H. Poincaré Anal. Non Linéaire 43 (2026)] we introduced a Riemannian metric $\mathsf{d}_S$ on the space of probability distributions obtained from entropic optimal transport, specifically from the Sinkhorn divergence $S_\varepsilon$. In the present work we discuss how to approximate and compute $\mathsf{d}_S$. Spatially, we prove Gromov–Hausdorff convergence of the metric and convergence of geodesics for increasingly fine Eulerian discretization of the base space. Temporally, we show $\Gamma$-convergence of the chain discretization $N \sum_{k=0}^{N-1} S_\varepsilon(\mu_k, \mu_{k+1})$ to the energy functional defining $\mathsf{d}_S$. We deduce and implement numerical schemes to compute approximations of $\mathsf{d}_S$.

We make highly accurate numerical estimates of basic properties of the invariant measure for Liverani-Saussol-Vaienti intermittent maps, an archetypal slowly mixing dynamical system. We solve the challenge of precise and efficient estimation for a map with non-uniform expansion, where the transfer operator does not have a spectral gap. We do this using an Abel function, which solves the neutral dynamics, and which we can compute accurately via an asymptotic expansion. Our work covers both finite and sigma-finite invariant measure cases. In particular, we obtain the first practical estimates for the parameter far from zero, including the sigma-finite case. This opens the door to much deeper numerical study of intermittent dynamics than was previously possible. A driving motivation for this work is the upcoming proof of optimal rates of memory loss which works equally for finite and infinite invariant measure cases. No prior results on rates of memory loss for dynamical systems with infinite invariant measure exist, with a single exception of a recent paper by I.Chevyrev and A.K. on null recurrent Markov chains.

We develop a variational framework for state-dependent diffusion on adaptive multiresolution representations in which the diffusion operator is generated by the adaptive state itself. The state consists of an admissible multiresolution tree, its active approximation space and basis, and the corresponding coefficient representation. It determines a symmetric nonnegative interaction form and an associated positive semidefinite intrinsic diffusion operator. In contrast to classical adaptive wavelet methods, where a prescribed operator is represented on an evolving approximation space, refinement and coarsening here modify simultaneously the representation, interaction graph, and operator. Because the adaptive hierarchy evolves through discrete topological changes, the coupled dynamics are formulated through a time-discrete variational principle rather than a differential evolution on a fixed space. We establish existence of the discrete updates, a discrete energy inequality, the coefficient-space null mode, and contractivity for frozen adaptive states. For regularized inverse problems, the construction yields Adaptive Multiresolution Diffusion Imaging (AMDI), combining data fidelity, intrinsic diffusion, coefficient sparsity, and tree complexity in a state-dependent energy. Numerical experiments verify the assembled operator identities, examine refinement-commutator decay, and confirm discrete energy dissipation. Adaptive Haar and higher-order multiwavelet calculations demonstrate localization of resolution on heterogeneous data. In denoising, AMDI retains high structural reconstruction quality with less than 10% of the full active representation, with stable behavior across held-out noise realizations.

We develop Toeplitz infinite GMRES for solving large sparse analytic parameterized systems $A(s)x(s)=z$ at many parameter values. The method exploits a block upper triangular Toeplitz structure in the companion Krylov sequence to construct the Arnoldi process and recover solution approximations without storing full Arnoldi vectors. We derive incremental recurrences requiring $\mathcal{O}(np^2+p^3)$ arithmetic and $\mathcal{O}(np+p^2)$ storage for $p$ Arnoldi steps, excluding factorization setup and assuming linear-cost coefficient actions and triangular solves. A dynamic generator-refreshing strategy addresses cancellation in the basic recurrence. For matrix functions admitting a fixed separated representation, we further develop an implicitly rebased method based on a compact, contractive nilpotent matrix. Both variants preserve the Arnoldi process in exact arithmetic, with exact matrix-function actions required for the rebased method. A residual lower bound for square-summable Taylor coefficients explains a scaling obstruction outside the normalized unit disk. Theoretical analysis and numerical experiments illustrate the computational efficiency of the proposed methods.

We present two particle methods for point set registration in bounded domains based on the Fokker-Planck equation. The first method relies on a moving least squares discretization of the S–formulation, in which moving grid points (particles) are advected by the drift with the target distribution, diffusion is resolved on a dynamically evolving particle cloud. This setting naturally leads to strong compression and expansion of the particle cloud. Obstacles are handled by enforcing reflective boundary conditions through a novel ghost point method. The second method is a Monte Carlo solver for the associated Langevin stochastic differential equation. It relies on a local approximation of the logarithmic gradient of the evolving particle density to extract macroscopic osmotic paths from individual stochastic trajectories. Owing to its inherent parallelism, this method exhibits excellent scalability and is well suited for high-dimensional registration problems. We illustrate the main features and performance of both approaches through extensive numerical experiments.

We prove asymptotic bias bounds for stochastic gradient Langevin dynamics in Wasserstein distance of order two. We assume that the negative log-density is strongly convex with a Lipschitz gradient, and that the stochastic gradient estimator is unbiased with an error satisfying a mean-square Lipschitz condition. The bounds are of order $h$ under a fourth moment assumption on the stochastic gradient error and of order $h^{1/2}$ under only a second moment assumption, where $h$ is the stepsize. A spiked-noise example shows that a second moment assumption alone is insufficient for a bound of order $h$ that is uniform over noise distributions with a fixed variance.

This paper presents a cell-centred unstructured finite-volume method for time-harmonic two-phase acoustics. A volume-fraction representation supplies density and acoustic compressibility. Face fractions average available neighbouring planar interface cuts independently of velocity direction, retaining the established one-field face-density closure. Cartesian complex-coordinate stretching provides perfectly matched layers (PMLs). Real-imaginary splitting yields a real-valued block system, discretized with corrected non-orthogonal fluxes and consistent boundary contributions. Verification covers homogeneous waves, layered gas-liquid transmission with PML truncation, baffled-piston radiation with Kirchhoff far-field reconstruction, and resolved rigid-sphere radiation forces. Homogeneous-wave pressure converges at approximately second order on orthogonal meshes and on meshes with non-orthogonal interiors and orthogonal boundary cells. Reconstructed velocity converges at approximately second order on orthogonal meshes and with an observed order of about 1.7 on the latter mesh family. Under aligned refinement, the layered case's relative whole-domain complex-pressure $L_2$ error decreases to $2.657\times10^{-4}$. The resolved-sphere force differs from the Gorkov prediction by at most 1.5% over the tested Rayleigh size range. The results quantify the accuracy and current limitations of the heterogeneous Helmholtz-PML formulation.

We design and analyse an unfitted hybrid high-order (HHO) method for the elastic wave equation in a medium made of two components separated by an imperfect interface of linear slip type, across which the traction is continuous and the displacement jump is proportional to the traction through a compliancy tensor $\bK=\alpha\bI+(\beta-\alpha)\bn\otimes\bn$. The mesh is not fitted to the interface: the discrete unknowns are doubled in the cut cells, the small cuts are cured by a cell agglomeration procedure, and no unknown is attached to the interface. The two specific ingredients of the method are a local symmetric strain reconstruction in each cut subcell, which incorporates the interface condition through the regularised interface stiffness $\bS_h=(h_T\delta^{-1}\bI+\bK)^{-1}$ in the spirit of Hansbo and Hansbo {\em{A finite element method for the simulation of strong and weak discontinuities in solid mechanics.}} {Comput. Methods Appl. Mech. Engrg.}, 193, 2004, and an interface stabilisation built from the same matrix. For the space semi-discrete problem we prove that the discrete bilinear form is coercive and continuous, and we derive an energy-error estimate of order $h^{k+1}$ and an $L^2$-error estimate of order $h^{k+2}$, with constants independent of the compliancy parameters and of how the interface cuts the mesh. The scheme is combined either with the Newmark scheme, which conserves a discrete energy exactly, or with singly diagonally implicit Runge–Kutta schemes of order up to four. Numerical experiments in two dimensions confirm the predicted convergence rates for $k\in\{1,2,3\}$, the robustness with respect to the compliancy over sixteen orders of magnitude, and illustrate the propagation of elastic waves across an unresolved slipping interface.

Numerical investigations of dynamo action have remained active in fluid mechanics over the past decades, presenting numerous challenges and open problems. Meanwhile, the development of structure-preserving methods and finite element exterior calculus (FEEC) inspires a revisit of numerical dynamo studies and an exploration of existing open questions in computation. In this paper, we present a FEEC approach for dynamo problems. In particular, we investigate structure-preserving finite element schemes for computing the spectra and pseudospectra of advection-diffusion operators of differential forms. The schemes and their analysis are based on finite element de Rham complexes.

We consider two coupled linear heat equations on different spatial domains that interact through a lower dimensional interface. This models conjugate heat transfer. The problem is solved using Dirichlet–Neumann waveform relaxation. This allows the subproblems to be solved using separate codes, a so called partitioned approach. Our overall goal is to develop more efficient partitioned methods, and to this end, we want reliable error estimates. Here, we use an exponentially weighted Fourier technique to derive new error estimates in $l^2$ for finite time $T$ in the fully discrete setting. These describe both linear and superlinear behavior. We show that the fully discrete estimate is close to a previously obtained time discrete estimate and independent of $\Delta x_1$ and $\Delta x_2$ when the CFL number in the subsolvers is large. We also show that the convergence behaviour depends on the ratio $\Delta x_1/\Delta x_2$ when the CFL number is small. Our numerical experiments show that the fully discrete estimate is accurate across a wide range of grid sizes $\Delta x_1, \Delta x_2$, time step sizes $\Delta t$ and $T$.

In this work, we study the pricing of American options under stochastic local volatility (SLV) models extended by including stochastic correlation driven by an additional stochastic process. We generalize the class of SLV models by incorporating a flexible stochastic correlation structure. To price options within these extended models, we derive the corresponding reflected forward-backward stochastic differential equations (RBSDEs) and employ data-driven numerical methods to solve them for both pricing and hedging purposes. The RBSDE framework enables the modelling of the future evolution of the option price. Furthermore, we conduct a convergence analysis of the proposed numerical method and present numerical experiments that illustrate the performance of the extended models, as well as the accuracy and efficiency of the RBSDE-based approach.

We propose an antithetic coupling that reduces the computational cost of gradient random-walk methods for one-dimensional convection-diffusion equations. The methods represent the solution's spatial derivative by particles with positive or negative masses, whose cumulative sums reconstruct the solution, and model diffusion by Gaussian displacements. The particles reorder and, for a nonlinear flux, interact. Antithetic sampling negates the displacements of a second simulation, which positively correlates the contributions of particles matched by position when their masses have opposite signs. Our coupling instead gives matched particles of opposite sign identical displacements and retains opposite displacements for equal-sign pairs. We prove that each simulation keeps the distribution of an uncoupled run as the matching changes, so the averaged solution keeps the bias of one simulation. For fixed particle configurations and matching, the coupling minimizes each matched pair's conditional variance contribution after diffusion among all couplings with the same Gaussian marginal distributions. From the same configurations, its variance reduction over antithetic sampling at the final diffusion stage has a closed form for any flux. For the heat equation, a martingale decomposition extends both results to each step's contribution to the terminal variance. In a Burgers' equation test, the coupling meets a tolerance for the spatially integrated mean squared error with 4 simulations, against 34 for antithetic sampling and 65 for independent sampling. In a second test, it meets a common error tolerance in about 40% less time than pairing particles only within each sign.

We present a novel numerical integration method for non-periodic functions over the high-dimensional unit cube $[0,1]^s$ by a specially crafted (unequally-)weighted quadrature rule using transformed lattice points, with the additional option of subsampling. The method is designed for integrands whose mixed derivatives up to order $\alpha$ are square-integrable with respect to a product Chebyshev density. The ingredients are (i) periodizing the integrand, (ii) approximating the smooth periodic function by a kernel interpolant at lattice points, (iii) integrating exactly the product of the kernel interpolant and the nonsmooth density, together with subsampling alternatives for (ii) and~(iii). With $|J|$ corresponding to the number of function evaluations, we achieve the convergence rate close to the order $|J|^{-(\alpha-1/2)}$ and $|J|^{-\alpha}\,($1.110721$)^s$, where the implied constants are independent of $s$ under favorable conditions. We can interpolate these two results to trade between the convergence rate and the growth in $s$. We provide numerical experiments to demonstrate our theory.

The scalar auxiliary variable (SAV) method replaces the nonlinear part of the free energy by a positive scalar $r(t)=\sqrt{\mathcal{E}_{_\mathcal{N}}[\phi]+C}>0$, thereby yielding linear, unconditionally energy-stable schemes for gradient flows. At the discrete level, however, the standard backward Euler and Crank–Nicolson discretizations provide no guarantee that the computed $r^{n+1}$ remains positive, an inconsistency with the continuous definition that contradicts the square-root ansatz and may compromise long-time robustness. Although the SAV method has been widely applied, this subtle but consequential issue has received little attention. We first characterize this failure quantitatively by deriving a sharp criterion and a sufficient condition on the time step size, and construct an explicit counterexample showing that sign loss occurs for parameters of practical relevance. Rather than modifying the definition of $r$ as in existing positivity-preserving variants, we retain the square-root form and reformulate the discrete evolution from $r_t$ to $(r^{2})_t$, which converts the scalar equation into a convex quadratic with a strictly negative constant term, always yielding a unique positive root. For the Crank–Nicolson scheme, the product-form discretization $r^{n+1}r^{n}$ preserves this quadratic structure, while conventional alternatives do not. The resulting schemes incur the same computational cost as the original SAV method and are proved unconditionally energy-stable. Numerical experiments for the Cahn–Hilliard equation confirm the predicted positivity, energy stability, and convergence rates.

We present a continuous Galerkin, Deferred Correction (cG-DeC) method for the equations of Lagrangian hydrodynamics. The proposed scheme combines a high-order continuous finite element discretization with an explicit Deferred Correction time integrator, yielding a formulation that avoids the inversion of a global sparse mass matrix at each update. To stabilize the method in the presence of strong shocks we introduce a modification of the hyperviscosity model that is compatible with the cG-DeC framework and likewise does not require a global solve to construct the finite element approximation to the polyharmonic operator. We further present an alternative reformulation of the DeC iteration which admits a simple recursive algorithmic structure and provides insight into previously observed convergence behavior of explicit DeC applied to hyperbolic partial differential equations. Conservation of momentum and total energy of the resulting fully discrete scheme is analyzed. A set of numerical experiments illustrates the accuracy and robustness of the cG-DeC scheme. Comparisons with a continuous Galerkin Runge-Kutta (cG-RK) method are provided.

Iterative methods are extremely popular for solving linear ill-posed inverse problems. While deterministic convergence, or more precisely, semi-convergence of such methods is widely studied, the problem of quantifying uncertainty in the resulting reconstructions due to random noise in the observed data has received far less attention. In this work, we study uncertainty quantification for the Landweber iteration in linear inverse problems using the Radon transform as a motivating example. We interpret the Landweber reconstruction statistically by analyzing how uncertainty in the measured data is propagated through the reconstruction map. This perspective is closely related to generalized fiducial inference, where uncertainty about the parameter is induced by inverting the relation between the observed data and the unknown (fixed) quantity. We combine the resulting stochastic uncertainty with a theoretical bound on the regularization bias to construct confidence intervals for the reconstructed solution. The evaluation of the resulting uncertainty estimates require repeated operations with the forward operator that can become expensive in large-scale problems. To reduce this cost, we use randomized singular value decomposition to obtain a low-rank approximation of the normal operator associated with the Radon transform and incorporate this approximation in quantifying uncertainty. {{In particular, this requires that the additional randomization error introduced due to the use of such randomized techniques be included in our approach for uncertainty quantification.}} Numerical results show that the proposed method provides accurate reconstructions and reliable uncertainty estimates at substantially reduced computational cost.

The universal approximation property of dropout neural networks does not by itself describe the network size required for an accurate random realization. In this work, we study approximation of the unit ball of $W^{n,\infty}([0,1]^d)$ by ReLU networks whose edges are retained independently with probability $p$. The approximation error is measured uniformly over the input domain, and the guarantee holds with probability at least $1-\delta$ for a single sampled network. We construct networks of constant depth and size $\widetilde O_{n,d}(p^{-9}\varepsilon^{-\max\{d/n,2\}} \log(1/\delta))$. The construction combines bounded local subnetworks, localization on a successful approximation event, and a multiscale Taylor decomposition. Conversely, Sobolev capacity imposes a lower bound on the number of surviving edges, while approximation of a fixed affine function requires an output-layer cost of order $((1-p)/p)\varepsilon^{-2}\log(1/\delta)$ at sufficiently high confidence. For fixed $p\in(0,1)$ and $\delta<\min\{1/2,1-p\}$, the upper and lower bounds match in the accuracy exponent under a fixed or logarithmic depth budget. When $d\leq2n$, they also match in confidence up to logarithms of accuracy. We extend the lower bounds to $W^{n,r}$ targets with $L^s$ error, and distinguish this extension from the upper bound for $W^{n,\infty}$. The optimal retention dependence and logarithmic factors remain open.

This paper develops a new moment-matching order reduction method for multi-input quadratic-bilinear descriptor systems. The proposed approach accommodates multiple inputs acting on both the differential and algebraic equations and constructs separate projection spaces associated with different input channels and input combinations. The method uses a matrix-free algorithm, allowing the reduced-order matrices and tensors to be computed without explicitly assembling or storing the full-order system matrices and tensors. This feature makes the proposed method particularly suitable for large-scale computational fluid dynamics problems. Numerical experiments carried out for two-dimensional flow problems demonstrate that the resulting reduced-order models accurately reproduce the transient responses of the full-order models under multiple time-varying inputs.

This work introduces a hyperbolic relaxation system designed for the simulation of kinetic equations in the low-Mach number limit. The methodology is built upon a micro-macro decomposition of the scaled BGK model, which reformulates the kinetic distribution function into a coupled system consisting of a macroscopic equilibrium part and a microscopic non-equilibrium remainder. By projecting the microscopic deviations onto a set of orthogonal polynomials, we derive a closed moment relaxation system. The resulting system is a version of Grad's 13 moment system with a linear hyperbolic part and a relaxation adapted to the incompressible limit. We prove the model's structural stability in the incompressible Navier-Stokes limit. Moreover, we develop a high order Asymptotic-Preserving (AP) numerical framework using Implicit-Explicit (IMEX) Runge Kutta schemes for temporal accuracy and finite difference WENO reconstructions as well as central difference approximations for high order spatial resolution. The proposed scheme ensures uniform stability and consistency across different physical regimes, automatically degenerating into a consistent high order discretization of the incompressible thermal Navier-Stokes limit as the scaling parameter vanishes. Numerical experiments in one and two dimensions corroborate the theoretical findings, demonstrating significant reductions in computational cost and robust performance across a wide range of Knudsen and Mach numbers.

We develop multifidelity methods for constructing triangular transport maps from samples, when high-fidelity data are scarce but lower-fidelity data are more abundant. Using this set of multifidelity data, we approximate a triangular transport map that bijectively maps between a tractable reference density and the high-fidelity target distribution. We introduce two strategies to leverage low-fidelity data: a hierarchical approach that composes maps between adjacent fidelity levels, and a non-hierarchical method that incorporates low-fidelity information through monotonicity-preserving corrections to the map parameterization. Numerical experiments compare these strategies with single-fidelity transport and demonstrate how the proposed multifidelity approaches can improve map estimation from limited high-fidelity data. To illustrate the broader utility of the learned maps, we also deploy them in a downstream amortized simulation-based inference task. This example shows that multifidelity improvements in map estimation can translate to improved conditional sampling and uncertainty quantification when high-fidelity data are scarce.

We study fast numerical methods for the Poisson equation on a disk within the FFTRR (Fast Fourier Transform Radial Recurrence) framework, which is built on Green's function representations. Classical FFTRR schemes apply FFTs in the azimuthal variable and evaluate mode-by-mode radial recurrences, achieving a complexity of \(O(MN\log N)\) on an \(N\times M\) uniform grid. However, they require a uniformly spaced azimuthal mesh of \(N\) points. In this work we develop a Nonuniform (NUFFTRR) solver that admits nonuniform grids in both the radial and azimuthal directions while retaining the favorable structure of the original FFTRR formulation. The azimuthal analysis-synthesis step is implemented using either a dense NUDFT least-squares solver or one of two NUFFT-based iterative schemes: a Toeplitz solver using circulant-preconditioned conjugate gradients (PCG), and a preconditioned conjugate gradient for least squares (PCGLS) solver with Pipe–Menon density compensation. These yield azimuthal complexities \(O(N^3 + MN^2)\) for the NUDFT variant and \(O(K_{\mathrm{iter}} MN\log N)\) for the NUFFT-based variants on an \(N\times M\) grid and Krylov iteration count \(K_{\mathrm{iter}}\). Numerical experiments on nonuniform meshes demonstrate that the proposed method is robust, spectrally accurate in the azimuthal variable, and competitive in runtime with existing fast Poisson solvers. We make use of vectorization and batched BLAS/GPU-accelerated linear algebra operations to eliminate explicit loops while also formulating azimuthal and radial steps entirely in terms of dense array operations, FFTs, and NUFFTs, allowing for straightforward GPU acceleration. The implementation is released as an open-source Python package at https://github.com/CharliePyle4/NUFFTRR_Poisson, and its methodology can be extended directly to related elliptic problems such as the Helmholtz equation.

We describe adapol, a Python package for fitting Matsubara Green's functions, self-energies, and hybridization functions by a compact sum of simple poles with real pole locations and matrix-valued residues. This pole-fitting step arises in many contexts, including hybridization fitting for quantum impurity solvers based on exact diagonalization, perturbation theory, and tensor networks, the compression of existing pole expansions such as the discrete Lehmann representation, diagrammatic methods based on discrete pole representations, and the analytic continuation of Matsubara Green's functions. adapol uses a modified version of the AAA rational approximation algorithm to estimate pole locations, and non-convex optimization to refine them. This procedure yields accurate and compact fits in a black-box and noise-robust manner. adapol provides a simple, self-contained interface with extensive documentation and examples, as well as an interface to the TRIQS package.

We study the inversion of the geodesic ray transform from finitely many local averages of its data. Under suitable geometric assumptions, we prove Lipschitz stability on any fixed finite-dimensional reconstruction space when the measurement partition is sufficiently fine. We develop convergent reconstruction algorithms based on Steepest Gradient Descent and Conjugate Gradients and implement them using piecewise constant finite element spaces. Numerical experiments in 2D and 3D illustrate the influence of geometry and measurement discretization on reconstruction quality.

We study a reduced multigrid Schur preconditioner for the fictitious-domain method with distributed Lagrange multipliers applied to elliptic interface problems. The method replaces the costly background solve by a multigrid approximation and treats the resulting correction only in the multiplier coupling space through a Woodbury formula. Numerical experiments show a marked improvement over the baseline preconditioner and consistently rapid convergence under mesh refinement and coefficient contrast.

We develop a pressure-free virtual element method for the surface Stokes problem on polygonal approximations of closed surfaces of arbitrary genus. The method extends the planar VEM Stokes complex to discrete surfaces, yielding a nonconforming complex with an exactly tangential, $H(\operatorname{div})$-conforming velocity space, commuting interpolation, and the correct discrete cohomology. To obtain the weak consistency needed for optimal convergence, we introduce an edge correction that preserves the complex without adding degrees of freedom; the resulting construction reduces to the original complex in the planar setting. The virtual formulation requires no penalties, and a local divergence-preserving reconstruction provides a computable pressure-robust load. We also construct explicit macroelement representatives whose energy and load induce an alternative formulation on the VEM space, algebraically equivalent to a direct macroelement Galerkin method. The original virtual formulation and the macroelement method are stable and attain optimal first-order convergence in the broken $H^1$ norm and second-order convergence in the $L^2$ norm for the velocity. Numerical experiments illustrate these convergence rates, the role of the edge correction, and pressure robustness.

We introduce a simple reaction-network model of a quorum-sensing population that couples the cell density $x$ to the signal density $w$. In the four-channel cell-signal network studied here, scaling signal production and removal by the same factor $r$ leaves all deterministic equilibria, and their stability types, unchanged. Nevertheless, we show that $r$ shifts the quasipotential barrier $\Delta V(r)$ for rare transitions towards extinction, and thus, under metastable exit assumptions, the mean time to reach a fixed neighbourhood of the extinction state on the exponential scale $e^{N\Delta V(r)}$. We compute the barrier by minimization of the path action with the signal retained as a fluctuating coordinate, and we compare it with exact stochastic simulation of population-threshold crossing times regressed in $N$. Over a range of $r$, the minimum-action barriers satisfy $\Delta V(r)=\Delta V_\infty+O(1/r)$, where $\Delta V_\infty$ is obtained by eliminating the signal first. At $r=1$, the barrier is 55% larger than $\Delta V_\infty$. The saddle barrier $\Delta V(r)$ is also the least action needed to enter the basin of extinction, but the density threshold can be crossed more cheaply: at $r=0.5$ the cheapest crossing costs 7.7% less action, keeps the signal high, and is usually followed by recovery. Simulated arrival times in this neighbourhood, which include failed attempts, grow with slopes within two fitted standard errors of the saddle barrier at every tested rate, and within $0.004$ of it if the logarithmic prefactor term is omitted. Under either regression model they exclude $\Delta V_\infty$ at $r\le2$ by at least $4.4$ fitted standard errors, without using the action solver.

We construct exact representations of finite element functions of arbitrary polynomial degree on simplicial meshes in any spatial dimension by two-hidden-layer multilayer perceptrons (MLPs). The first hidden layer uses $\mathrm{ReLU}^0+\mathrm{ReLU}$, and the second uses the degree-dependent activation $\mathrm{ReLU}^k$. The neural networks can realize both zero-skeleton and prescribed-skeleton representatives of discontinuous finite element functions, where the former vanish on the mesh skeleton and the latter admit independently prescribed polynomial data on each relatively open lower-dimensional mesh subsimplex; in particular, continuous finite element functions can be represented pointwise on the closed domain. In the continuous case, a half-open decomposition reduces the second-hidden-layer width by assigning each subsimplex to a single incident element. The number of nonzero parameters grows linearly with the number of elements for fixed spatial dimension, polynomial degree, and output dimension, and it is of the same order as the number of global degrees of freedom of the finite element space in the continuous case, provided that the mesh is shape-regular. All parameters in the neural network can be explicitly computed via the information of the mesh and finite element function without training.

We derive a method of analytic function expansion using barycentric rational interpolation, in which an arbitrary analytic weight function defines the interpolation and nodes. Truncating these expansions yields successive approximations, and the derivation of the expansions gives the approximation error as a byproduct. We demonstrate the method with cosine and the Bessel functions of the first kind for several choices of weight functions, and we examine the accuracy and convergence of the approximate interpolations as a function of the weight and the number of nodes used in the interpolation.

A technical condition arising in the analysis of atomistic models for dislocations is studied. Key estimates for the far-field strain behaviour proved in Ehrlacher, Ortner, Shapeev (2016) rely on a summation-by-parts argument for linearised forces, under conditions of sufficient decay and of vanishing net force in an infinite system. In particular, the latter condition is required for an application of this theory to the standard far-field dislocation predictor, but the vanishing of the associated force sum has not been fully justified. Here, the missing verification of this condition is provided and its role in determining the far-field decay of the corrector is clarified. It is moreover shown that for more general physically compatible predictors, the linearised force sum need not vanish, leading to slower decay and potentially divergent finite-domain approximations.