Papers

A particle beam with constant intensity that propagates into a plasma triggers a bump-on-tail instability. It has been demonstrated that beam modulation can be used to significantly reduce this instability. However, finding the optimal beam profiles that accomplish this is a complex optimization problem that suffers from an extremely rugged optimization landscape and a relatively high-dimensional parameter space. In this paper, we evaluate a number of optimization algorithms and find that population based methods (in particular, CMA-ES and Differential Evolution) perform best. In particular, we investigate beam profiles where the beam is either completely turned ON or OFF. This configuration is easier to accomplish experimentally, but leads to a higher-dimensional and thus more challenging optimization problem. We obtain ON/OFF beam profiles that reduce the instability by more than 80% and thus perform significantly better than what has been reported in the literature.

We show the convergence of a commonly used mesh-free discretization of a dynamic non-linear bond-based peridynamic model with irreversible damage implemented via a bond-breakage criterion for a fixed horizon $\delta >0$. The model chosen here is a mathematically well-posed regularization of a discontinuous peridynamic operator often found in the literature, which assumes that fracture results from strain in tension only. We explicitly include the one-point quadrature used in this scheme in the analysis. The convergence is with respect to the natural norm of the solution space given by $C^1([0,T],L^\infty(\Omega)^d)$, and therefore $L^\infty(\Omega)^d$ in the spatial coordinates. This result heavily constrains the spatial regularity of the continuous solution which we prove afterwards. More specifically, for locally $L$-Lipschitz continuous input data, local Lipschitz bounds are derived for the displacement, regardless of damage. This is consistent with interpreting the peridynamic model used here as an approximation to a truly discontinuous fracture model with a length-scale parameter given by the horizon $\delta>0$. We note that we do not deal with the limit of the solutions as $\delta \searrow 0$.

The escalating environmental footprint of the Information and Communication Technol- ogy (ICT) sectorwith global Data Center electric- ity usage potentially reaching 3% by 2030 and the imperative for real-time responsiveness in smart buildings demand a foundational shift towards Edge Computing. This article proposes a holis- tic architectural and algorithmic approach for cre- ating autonomous, resilient, and privacy-preserving learning systems within smart buildings, advancing the core concept of the building itself as the data center. Our contributions, demonstrating originality compared to conventional centralized systems, are structured around three main pillars: 1) Validation of resilient data acquisition chains using the af- fordable Arduino ecosystem (Uno R4 WiFi/BME280) and secure persistence solutions including centralized archiving (Node-RED/MySQL) and secure remote transfer (Python/rsync/SSH). 2). The proposal of a lightweight, KOptim orchestration framework, inter- operable with the FIWARE standard, for the dynamic management of Service Level Agreements (SLA). 3) The adaptation and experimental valida- tion of the Generalized Hebbian Algorithm (GHA) for Principal Component Analysis (PCA) on low- power microcontrollers,We introduce and compare an iterative learning strategy using data buffers, featuring sorting and explicit preservation of extreme samples to maintain diversity. Experimental results on the high-dimensional Tour Perret sensor data confirm that this novel method, despite introducing a necessary temporal complexity of O(W log W ) due to sorting, achieves convergence and final result quality that is visually equivalent to processing the entire dataset in full memory, reliably identifying two distinct operational states.

We present high order accurate mixed explicit-implicit discontinuous Galerkin (DG) schemes for the solution of hyperbolic systems of conservation laws on Cartesian cut cell meshes. Such meshes may contain arbitrarily small cut cells, which leads to the well-known small cell problem: the CFL condition of explicit schemes would severely restrict the admissible time step. To overcome this difficulty, we use high order explicit fully-discrete one-step ADER-DG schemes only on Cartesian cells away from the embedded boundaries. The zone close to the embedded boundary, which consists of cut cells and their direct face neighbors, is discretized with an implicit space-time DG scheme that is coupled to the explicit scheme directly via the numerical flux. The resulting scheme, called ADER-ZIST (ADER-DG with Zonally Implicit Space-Time DG), is conservative by construction, and its time step is determined only by the size of the Cartesian background cells. Curved embedded boundaries are represented by high order isoparametric sub-triangles inside the cut cells, and the nonlinear system of the implicit cells is solved via an inexact matrix-free Newton-Krylov method. For problems with shock waves, we add an artificial viscosity term whose strength is controlled by an oscillation-eliminating smoothness sensor. ADER-ZIST is assessed on a comprehensive set of test problems for the compressible Euler equations, ranging from smooth flows around straight and curved embedded boundaries to acoustic scattering and shock-dominated flows. The scheme essentially reaches its nominal order of accuracy N+1 for smooth solutions and for polynomial degrees up to N=5, on meshes with minimal volume fractions of the small cut cells down to the order of 10^-7.

We present an explicit quantum circuit that constructs a block encoding of a generator-representable semiseparable (one-pair) matrix directly from its two generators. For an $N\times N$ semiseparable matrix with $N=2^{n}$, the circuit uses $2n+2$ ancilla qubits and has $O({\rm polylog} N)$ depth, and it realizes a $(2\sqrt{N},\,2n+2,\,2 \sqrt{2}\ \varepsilon)$-block encoding, where $\varepsilon$ is controlled by the accuracy of the state-preparation routines that load the generators. The building blocks of the circuit are cyclic shifts, swaps, and Hadamard gates, each admitting an efficient implementation through the binary decomposition of the relevant indices.

Extended-domain spectral collocation expands the solution of a boundary value problem in the Laplacian eigenbasis of an enclosing box and collocates at as many nodes as modes. The square system is exponentially ill-conditioned, and the impossibility theorem for equispaced data rules out stable geometric convergence; yet such solvers converge geometrically to a few orders of magnitude above the unit roundoff. We explain why in one dimension. The collocation solution is the solution operator applied to the interpolant of the data in an operator-dependent test space whose Lebesgue constant must grow exponentially for geometric convergence, so exact solutions can amplify data errors exponentially. In floating-point arithmetic the error of any coefficient vector is bounded by its nodal residual, its distance to a competitor and two frame constants, a leakage and a restricted Marcinkiewicz–Zygmund constant, bounded uniformly in the number of modes at a fixed positive threshold; for a backward-stable solve, under a condition on the perturbed matrix, the bound is uniform in the conditioning for suitably analytic data and, with computed constants, within a factor of about ten of the error floor of a direct solver. For the Poisson operator the rate of the two constants along the tail of the spectrum is given by a constrained equilibrium problem of logarithmic potential theory. For an interval of length $L$ extended by $\delta$ at each end, the problems for $\delta$ and $L^2/(4\delta)$ are dual, the rate is linear for $\delta = L/2$, and it starts with slope minus the fraction of the extended domain occupied by the extension, so for $\delta < L/2$ the leakage constant obeys no power law with the exponent of the end of the spectrum. In experiments, convection–diffusion operators behave like the Poisson operator at moderate Péclet numbers.

In this paper, we propose a linear and decoupled fully discrete scheme that preserves the entropy production property for the dynamic van der Waals model. This model provides a thermodynamically consistent diffuse-interface framework for describing non-isothermal phase transitions in single-component liquid-vapor systems, containing the coupled conservation laws of mass, momentum and energy. We first introduce an equivalent reformulation of the governing equations that avoids the use of chain-rule identities, which generally fail under discrete differentiation. Further, a component-wise convex-splitting method is employed to calculate the discrete chemical potential and temperature. This treatment provides a discrete entropy-stable inequality without a time step restriction. Based on the explicitly treated velocity, the scheme allows the three conservation equations to be solved linearly and sequentially in a decoupled manner. Since the momentum equation does not directly affect the entropy-stability analysis, we use a classical linear semi-implicit scheme to solve it. Numerical experiments verify the expected numerical accuracy. Simulations of droplet evaporation and two- and three-dimensional Leidenfrost effects further demonstrate the proposed method can capture complex non-isothermal phase-transition dynamics while numerically reproducing the discrete entropy production behavior.

This is a lecture note to introduce the innovative concept of numerical smoothness and its application in the error analysis of numerical solutions of evolution equations. We start from the two fundamentally different error splittings. One splitting is popular, which brings us to local error defined along the original solution and error propagation by numerical schemes. The well know Lax Equivalence Theorem is based on this error splitting. The other splitting can only be seen in a few numerical ODE books, which brings us to local error defined along the numerical solution and error propagation by the solutions of the evolution equations. The innovative concept of numerical smoothness arises from the second error splitting. While there are many difficulties when we use numerical stability in error analysis of nonlinear evolution equations, we will show that numerical smoothness works naturally in error analysis/estimation/control with advantages in several ways. We will define temporal and spatial numerical smoothness indicators, prove the necessity of the boundedness of these indicators for having optimal convergence rate, and show computable error estimates as functions of these smoothness indicators. Moreover, these error estimates are optimal in the approximatons' convergence rates and also optimal in the error propagation of the evolution equation. In comparison with error analysis framework of the Lax Equivalence Theorem, from the necessity results and the error estimates proven in the numerical smoothness framework, we have a new equivalence theory: Consistency + Numerical Smoothness <=> optimal error estimates. Unlike the Lax Equivalence Theorem, the new equivalence theory works well for nonlinear problems, nonlinear schemes, and many other complications of problems.

Sobolev gradient flows for the Gross–Pitaevskii ground state problem are Riemannian gradient descent methods on the sphere under a metric induced by a discrete elliptic operator. We analyze their convergence through the condition number of the Riemannian Hessian at the ground state. Bounds of the Rayleigh quotient of the Riemannian Hessian give a uniform analysis of these methods for different metrics. The uniform analysis applies to the classical metrics, such as the shifted Laplacian and the energy-adaptive metrics. Inspired by this analysis, we also consider three new metrics. One of them is a Thomas–Fermi metric, whose condition number is bounded uniformly in the interaction strength for the one-dimensional harmonic trap and whose iteration count does not grow in the Thomas–Fermi regime in the numerical tests. The performance suggested by our analysis is consistent with three-dimensional numerical tests with up to one billion unknowns, for both Riemannian gradient descent and Riemannian conjugate gradient methods, whose practical performance is also related to the Hessian condition number. The Thomas–Fermi flow needs the fewest iterations, while the metrics inverted exactly by tensor-product solvers are the fastest in computational time on one GPU card.

The precorrected-FFT has been used to rapidly evaluate $N$-body sums arising in electromagnetics and acoustics settings. It has, however, been historically limited to kernels that satisfy a Green's identity. In this work, we demonstrate that, using an improved version of the method of proxy points, the precorrected-FFT can be used to efficiently compute $N$-body sums with a wide variety of kernels. This algorithm requires no analytic information beyond pointwise access to the kernel. Our examples include kernels from non-local PDEs, Matern kernels, anisotropic kernels, and the Green's functions of high order PDEs. We also provide a detailed algorithm for choosing the proxy points and spreading grid, as implemented in the new pcfft package. We show these algorithmic choices result in favorable scaling when compared to other kernel-independent fast summation methods. Finally, we use the algorithm to accelerate a wide variety of integral equations solvers involving non-standard kernels, with applications in surface scattering, geophysics, and microbial community formation.

Approximating the effective anisotropic hyperelastic response of a representative volume element (RVE) requires constitutive expressiveness, mechanical consistency, and computational efficiency. This work develops a polyconvex physics-augmented neural network (PANN) construction based on paired directional features whose orientations, stretch powers, and determinant exponents adapt to the RVE response. The features are composed with convex, componentwise nondecreasing cores, realized as an input-convex neural network (ICNN) or an input-convex Kolmogorov-Arnold network (ICKAN). The complete energies satisfy objectivity, an exactly stress-free reference, two-dimensional polyconvexity, and prescribed growth conditions by construction, and nonnegative energy is certified separately for the selected models. On a multicavity and a single-cavity periodic cell, these energies are compared with an Unconstrained energy, trained identically without the curvature and growth constraints, and with three hyperreduced projection-based reduced-order models (HPROMs), which serve as physics-based references. On the first cell, learning the features lowers the stress error sevenfold, and six features suffice on both. The constrained energies remain below one percent in stress error and in the displacement error of a tensile-coupon FE$^2$ application with a full-order reference, whereas the Unconstrained energy and the HPROMs are more accurate inside the data. Beyond the data, however, the Unconstrained energy loses rank-one convexity, and under confined compression it stops below half of the 500 MPa that the constrained energies reach. These results quantify what the guarantees cost in accuracy and what they secure beyond the data.

We develop a neural transport framework for nonlinear, nonlocal aggregationdiffusion equations arising in chemotaxis, polymer dynamics, and related applications. Based on the Onsager variational principle, we reformulate the governing equation as a constrained optimal control problem, taking the velocity field as the control and evolving the density via the associated continuity equation. The velocity field is parameterized by a neural network trained on the variational objective evaluated along particle trajectories. By coupling the density and generalized score equations into a single system solved via the Euler method, the score function propagates efficiently for loss evaluation. Furthermore, we establish a rigorous error analysis that decomposes the total error into approximation, optimization, Monte Carlo sampling, and time-discretization components, and prove corresponding error bounds. Numerical experiments across several representative models demonstrate the accuracy and effectiveness of the proposed approach.

The standard Yee finite-difference time-domain (FDTD) algorithm uses Maxwell's equations to replace first order temporal derivatives with spatial curls. In this work, the relation is extended to arbitrary orders: a Time-Space Theorem establishes the correspondence between higher-order temporal derivatives and successively higher-order spatial curl operators. Based on this relation, an arbitrary-order framework for electromagnetic field approximations is developed. The framework constructs higher temporal orders through higher curl orders and permits different spatial discretizations and curl-estimation techniques, thereby providing a general basis for developing numerical algorithms of arbitrary orders. As one implementation, an arbitrary-order FDTD algorithm is constructed using staggered time and space stepping together with recursive curl estimations. Numerical simulations of an initial-value problem are performed for temporal approximation orders from second to tenth order while the spatial and temporal step sizes are kept fixed. Comparisons with the analytical solution show that the errors decrease rapidly and monotonically with increasing approximation order. These results demonstrate the practical potential of the framework for high-accuracy electromagnetic computation and provide a basis for further development of arbitrary-order numerical methods.

The repeated assembly of the sparse Jacobian matrix constitutes a primary computational bottleneck in the Newton-Raphson power-flow analysis. Through an algebraic equivariance analysis of the Jacobian matrix, this paper proposes a direct sparse Jacobian assembly method based on symbolic permutation. By pre-applying the permutation to the admittance matrix during a one-time symbolic phase, the runtime bus-type projection overhead is entirely eliminated. Consequently, complex Jacobian differentiation and global matrix assembly fuse into a single operation, yielding a single-pass, branchless algorithm that streams values directly into contiguous memory. The framework is realized in the open-source RustPower solver. Benchmarks on IEEE test systems and PEGASE 9241 demonstrate that the proposed assembly is 24.4-30.1x faster than the MATPOWER- style baseline and 6.0-7.8x faster than the Schafer-Braun optimized assembly. End-to-end, the grid-level solve time is 11.7x-154.3x faster than pandapower 3.5.4 with Numba, while leaving numerical results and convergence behavior unchanged. Compared with the native C++ LightSim2Grid 1.0.0, grid-level performance is comparable on IEEE 39 and 1.3x-1.8x faster on IEEE 118 and the PEGASE systems. The method extends directly to general problems where a sparse Jacobian inherits its structure from a static network operator, including transient stability analysis, EMT simulation, and OPF Hessian assembly.

Nested integrals $\int f\left(\int g(\boldsymbol{y},\boldsymbol{x})\mathrm{d}{}\boldsymbol{x}\right)\mathrm{d}{}\boldsymbol{y}$ where $\boldsymbol{y}\in[0,1]^{d_1}$ and $\boldsymbol{x}\in[0,1]^{d_2}$ frequently arise in Bayesian experimental design, medical decision making, and computational finance. Statistical estimators for nested integrals are more computationally demanding than for standard integrals. Recent studies demonstrate that the median randomized quasi-Monte Carlo (mRQMC) method achieves super-polynomial convergence rates for analytic integrands. We extend mRQMC methods to nested integrations and propose a novel double-loop mRQMC estimator with super-polynomial convergence rate if both the inner and outer integrands are analytic. Established estimators for nested integrals, such as double-loop or multilevel Monte Carlo, are by contrast limited to algebraic convergence rates. However, expected information gain (EIG) applications from experimental design typically do not fulfill this condition due to a lack of regularity in the outer integrand. We address this through a theoretically-grounded importance sampling regularization scheme. Although super-polynomial convergence rates remain achievable for the inner integral, only the algebraic rate $\alpha+3/2-\epsilon$ for any $\epsilon>0$ can be theoretically guaranteed for the outer integral in the EIG, where $\alpha\in\mathbb{N}\cup\{0\}$ denotes the order of the mixed derivatives of the regularized outer integrand. Our approach shows substantial improvements over double-loop randomized quasi-Monte Carlo methods based on the sample mean in moderately high-dimensional settings for a simple polynomial example. Moreover, we demonstrate improved convergence behavior of the median-based estimator for the more challenging setting of EIG estimation.

Physics foundation models are an increasingly popular approach for building versatile emulators for spatiotemporal systems. They require tokenisation layers that lift physical field data into a latent representation suitable for dynamic processing. A typical choice for this layer is a convolutional encoder, which lifts field information in a locality-aware way and can be trained very efficiently. However, this choice cannot be applied to a large class of scientific computing problems where data are represented on unstructured meshes, for example simulations on domains with curved boundaries. Such data typically require graph neural network architectures, which are more expensive to train and, in the form currently used, encode information differently from their CNN counterparts, thus limiting the portability of physics emulators between these data regimes. In this work we introduce Mesh-Native Convolution (MNC), a locality-aware tokenisation layer that encodes unstructured mesh-based data in a way that is analogous to CNNs. This approach leads to models that can be trained once on mesh data and then deployed to arbitrary unstructured meshes, even if these are irregular, with minimal retraining. We show that this tokenisation approach can be used to retrofit grid-pretrained foundation models efficiently on unstructured meshes and provide theoretical results to confirm the consistency of this approach with classical CNN encoders. We demonstrate the favourable performance of this retrofitting approach, which in many cases reaches acceptable accuracy even in zero-shot deployment, and outperforms existing baselines after only a few epochs of fine-tuning, thus providing an efficient bridge between the two data regimes.

We present a Fourier-based method for evaluating two-electron four-center Coulomb repulsion integrals over Slater-type orbitals. Spherical Fourier-Bessel coefficients of the convolution factors are precomputed and reused. Since the spherical part of the integral can be evaluated directly, there remains only a one-dimensional quadrature problem in the radial direction. For the ten reference integrals reported here, the largest absolute discrepancy is $1.7\times10^{-9}$. A preliminary resident-cache benchmark computes approximately $4.2\times10^5$ unique integrals per second on an NVIDIA V100 at the contraction stage.

We propose a novel numerical method for transport-dominated problems using Physics-Informed Neural Networks (PINNs). In such problems, the solution can develop large gradients and fine structures over time. This is particularly challenging for classical PINNs, which represent the solution itself. In our approach, rather than the solution, we choose to approximate the flow map of the equation. The solution is then recovered by pulling back the exact initial condition, so that its fine structures are produced by the map rather than represented by the network. Our method is constructed by approximating the flow map of the equation with multiple PINNs, one per time subinterval, and composing them using the semigroup property of the flow, so that each network only has to represent a map close to the identity. To highlight the advantages of our method compared to classical PINNs or classical numerical schemes, we present numerical results on the linear advection equation, the incompressible Euler equation in vorticity formulation, the Vlasov-Poisson equation, and finally a drift-kinetic equation.

We study a family of randomized alternating projection methods that alternate among subspaces spanned by subsets of \(n\) vectors from a prescribed set. We derive an explicit characterization of performance bounds when these subspaces are sampled with probabilities proportional to the volumes subtended by the vectors in their corresponding subsets. For each \(n\), the bound is obtained through an explicit nonlinear transformation of the spectrum of an associated matrix. Our analysis reveals an explicit spectral equalization mechanism that drives the spectrum toward improved conditioning as $n$ increases, establishing an unexpected connection between volume sampling and the Cayley-Hamilton theorem. Furthermore, we introduce a computationally efficient uniform sampling scheme that achieves comparable theoretical guarantees through a specific relaxation strategy. Empirical results show that this relaxation method achieves convergence rates comparable to, and often better than, volume sampling, while avoiding the infeasibility of volume sampling at large scale. Besides the randomized Kaczmarz, these results also directly apply to randomized Gauss-Seidel and coordinate descent methods.

Many tasks in scientific computing and machine learning require the Jacobian or Hessian matrix of a function. Automatic differentiation (AD) computes these derivatives to machine precision, but materializing a dense $m \times n$ Jacobian requires $n$ forward-mode or $m$ reverse-mode AD passes, one per column or row. For a large class of functions, each output depends on only a few inputs, making the derivative matrix sparse. Automatic sparse differentiation (ASD) exploits this structure in four steps: detection of the input-agnostic sparsity pattern, coloring of a graph to group columns or rows that can share an AD pass, compressed differentiation to compute a compressed derivative matrix with one AD pass per color, and finally decompression into the original sparsity pattern. The number of colors, and hence of AD passes, is often independent of the problem dimension: a banded Jacobian with $b$ contiguous bands, for instance, only ever requires $b$ colors, regardless of its size. asdex offers the first standalone ASD toolkit in the popular JAX ecosystem. With $\texttt{asdex.jacobian}$ and $\texttt{asdex.hessian}$, it provides sparse drop-in replacements for $\texttt{jax.jacobian}$ and $\texttt{jax.hessian}$.

Symmetric and unitary matrices are among the most important classes of structured matrices admitting efficient and stable eigenvalue algorithms. Notable examples include structure-preserving variants of the QR algorithm and their square-root-free counterparts. For unitary matrices, these algorithms were historically derived from recurrence relations for orthogonal polynomials on the unit circle. A matrix-based derivation using core chasing was later given in [3]; this family of algorithms is backward stable, but the underlying symmetries needed to cast out the square roots were not identified, and the question of structured backward stability was left open. In this paper we present a matrix-based derivation of the symmetric unitary core-chasing algorithm and use it to obtain a square-root-free variant. We give a matrix-based proof of backward stability and show that strict structured backward stability is impossible in general. Open-source Fortran code and numerical experiments demonstrate that the combination of exploiting the symmetry and removing the square roots reduces computation time without degrading accuracy.

Geodynamic simulations involve solving the Stokes equations, where boundary conditions play a crucial role in governing internal dynamics and surface topography evolution. A common approach is to model the interface in contact with air as a free surface. This interface is stress-free and deforms in response to internal force balances, which is essential for understanding processes ranging from rifted margins to subduction. However, a free surface can introduce instabilities and oscillations in the numerical solution, commonly known as the "drunken sailor" instability, forcing the use of small time steps. Geodynamic models therefore often employ stabilization terms based on the gravitational effect of free-surface motion. To date, these terms have been derived primarily on a heuristic basis. Here, we apply variational calculus to the energy-minimization formulation of incompressible Stokes flow. Using the Gateaux derivative, we derive the complete weak form for a deforming free-surface domain and show that the stabilization terms arise naturally from the free-surface assumption. We consider the linear-viscosity case and identify two boundary-integral contributions: one associated with internal viscous deformation and the other with gravitational body force. Two benchmark cases demonstrate the consistency of the derived formulation with previously proposed stabilization approaches. Our results provide a mathematical foundation for the commonly used gravity stabilization and identify an additional viscous boundary contribution.

Predicting the structure and response of inhomogeneous polymer electrolytes requires a description of ion correlations that retains molecular-scale accuracy while remaining transferable across spatial scales and geometries. We develop a neural density functional for electrolytes that preserves spatial symmetries, thermodynamic integrability and the Noether identities, with perfect screening recovered in stable, noncritical bulk states. Its nonlinear density dependence captures the concentration-dependent correlations missed by a pair closure, including a crossover from enhanced to suppressed long-wavelength number fluctuations at strong coupling. The functional describes density profiles at an untrained salt concentration and predicts bulk structure factors and the long-wavelength number response. Trained solely on planar density and internal-force profiles from molecular dynamics, the functional predicts ionic structure in larger domains and in two-dimensional external fields. On the same ion data, it is more accurate than three other neural density-functional architectures and keeps its accuracy with a quarter of the training runs, where the errors of the best alternative grow by about two thirds. The spatial transferability provides a necessary foundation for connecting molecular correlations to continuum predictions at larger scales.

Finite element simulations carried out on meshes issued from real-world geometries - patient-specific anatomies, scanned industrial parts, fractured media - almost never satisfy the shape-regularity assumptions under which classical a posteriori error analysis is written. Practitioners are therefore left without a quantitative answer to a very legitimate question: how wrong is my simulation when my mesh is coarse and contains badly shaped elements? In this work, we study to which extent equilibrated flux estimators can provide such an answer. Indeed, they deliver a guaranteed and fully computable upper bound of the energy-norm discretization error, with a constant equal to one; and, crucially for the present study, its derivation uses no shape-regularity, quasi-uniformity or angle condition. The only geometric ingredient is the Poincaré constant on a single cell, which Payne and Weinberger showed to be at most hK /pi on any convex cell K, however distorted (hK being the cell size). Mesh quality affects the local efficiency - the sharpness of the estimator - but never its reliability, so that the failure mode under distortion is pessimism rather than false confidence.

We consider reduced basis approximations of non-equilibrium steady states of parameter-dependent Lindbladians, whose ambient matrix space grows exponentially with system size. We give a unified structural view of the problem through its characterizations as an eigenvalue problem or as a trace-constrained source problem. Building on these characterizations, we propose a residual-driven reduced basis framework for both formulations. Under the assumption of a unique steady state, we derive residual-based state-error estimates, and establish (sub-)exponential convergence of the greedy algorithms under real-analytic parameter dependence of the problem data. Numerical experiments on boundary-driven Fermi–Hubbard and dissipative XYZ chains support the analysis and suggest that the trace-constrained greedy method offers the best overall balance of accuracy and computational cost.

In this work, we study the fast Fourier transform (FFT) of discrete data with crystallographic symmetry, for which standard FFT algorithms incur substantial redundant computation. We establish a rigorous mathematical theory of the symmetric fast Fourier transform (SFFT) for $n$-dimensional crystallographic groups and propose the SFFT and symmetric inverse fast Fourier transform (SIFFT) algorithms. We derive symmetry-reduction formulas for both linear and translational symmetries and unify them into a recursive framework. We analyze the arithmetic complexity to obtain theoretical speedup factors relative to the standard FFT, and we confirm the accuracy and efficiency of the proposed methods by numerical experiments on representative three-dimensional space groups.

Structure-preserving numerical methods for gradient flows have been extensively developed when the governing equations and free energies are known. However, accurately predicting dynamics from observational data while preserving intrinsic physical properties remains challenging. Although neural operators, such as Fourier neural operators (FNOs), provide powerful data-driven approximations, their predictions do not inherently guarantee physical structure preservation. To address this challenge, we integrate deep learning with structure-preserving exponential time differencing (ETD) methods. Instead of directly learning the evolution operator, we introduce a *trajectory-inferred free energy* (TIFE), which reconstructs the unknown energy density from observed trajectories through a Duhamel-based learning framework. The learned energy determines the variational force and stabilization parameter, enabling structure-preserving numerical evolution. We further extend the framework to jointly identify the free energy and diffusion coefficient. Under suitable assumptions, we establish the maximum-bound principle, unconditional dissipation of the original learned energy, and error estimates incorporating discretization and inference errors. Numerical experiments demonstrate accurate energy recovery and reliable long-time predictions.

In this paper, we propose a new approach called preflattening to speed up the resolution of nonlinear algebraic systems by Newton-type methods. Preflattening is meant to be for nonlinear systems what preconditioning is for linear ones. The idea is to transform the system into an equivalent one that is more favorable to Newton's method. To this end, we introduce the notion of flatness, which is intended to play for nonlinear systems a role similar to that of the condition number for linear systems. For a physical model of interest in the geosciences, namely the chemical equilibrium problem, we combine this strategy with the nonlinear elimination technique and demonstrate the value of the resulting approach through numerical simulations.

We present a boundary integral equation (BIE) framework for two-dimensional time-harmonic elastic transmission problems. Using the Helmholtz decomposition, the elastodynamic displacement fields are expressed in terms of scalar pressure and shear potentials, so that only the four classical Helmholtz boundary integral operators are needed. The key step is a rewriting of the Navier traction operator in the local normal–tangential frame of the interface, in which second-order normal derivatives are eliminated in favour of tangential derivatives and curvature terms. This leads to a one-parameter family of indirect combined field integral equations, depending on a coupling constant $\alpha\ne 0$, which we prove to be uniquely solvable. The principal part of the resulting operator is, however, defective, and the discretized systems are severely ill-conditioned. Using the pseudodifferential calculus of the Helmholtz BIOs, we construct explicit left and right analytical preconditioners and prove that the preconditioned operator is a compact perturbation of the identity. Numerical experiments based on high-order Nyström discretizations confirm that the preconditioner clusters the spectrum: the number of GMRES iterations becomes independent of the discretization and is reduced by one order of magnitude or more, also in the high-frequency regime and for non-convex scatterers.

We study the approximation by a finite volume scheme of the stationary Gross-Pitaevskii equation with angular momentum rotation. We prove that, under standard mesh regularity assumptions, the discrete global minimizers converge to the continuous one as the mesh diameter tends to zero. Moreover, under a spectral gap assumption on the Hessian of the discrete energy, we prove that solutions to a normalized continuous gradient flow converges exponentially fast to such discrete minimizers. Additionally, we show that this previous discrete spectral gap assumption can be derived from similar hypotheses in the continuous setting. We support our theoretical results with numerical simulations, which broaden the scope of the analysis.