Papers

Cryptocurrency markets exhibit extreme volatility and non-stationary dynamics that challenge conventional forecasting methods. Although Large Language Models (LLMs) have shown promise for time series forecasting, the combined effects of adaptation choices remain largely unexplored in financial settings. This study introduces PRICE, a structured approach for adapting LLMs to short-term Bitcoin price forecasting. Built on a 4-bit quantized LLaMA-3 8B model, PRICE investigates how fine-tuning, numerical representation, prompting, inference, and decoding jointly influence forecasting performance. PRICE integrates Parameter-efficient fine-tuning with Low-Rank Adaptation (LoRA), Recursive multi-step inference, Integer-rounded numerical representation, Context-Task-Format (CTF) prompting, and Exact zero-temperature decoding. Ablation studies show that each component contributes to forecasting accuracy and reliability. LoRA enables efficient training on limited hardware, recursive inference improves accuracy, integer-rounded values reduce errors, CTF prompting outperforms Chain-of-Thought, Implicit Chain-of-Thought (iCoT), and few-shot prompting, and zero-temperature decoding improves stability during recursive forecasting. Comparative evaluation against eight transformer-based and time-series foundation models shows that PRICE achieves the lowest forecasting errors on both validation and test sets while maintaining robust performance across evaluation periods. Despite being based on a model primarily pretrained on text rather than time-series data, PRICE achieves competitive or superior performance relative to specialized foundation models. These findings demonstrate that adaptation choices critically determine the accuracy and robustness of LLMs for numerical time-series forecasting.

We consider (feed-forward) Neural Networks (NNs) for the emulation of the solution to singularly perturbed second order boundary value problems, with two small parameters. We describe a shallow NN which exploits available asymptotic expansions for the solution. These additive decompositions into smooth and layer components, allow for derivative estimates which are explicit in the order of differentiation as well as the singular perturbation parameter(s) \cite{melenk, Irene, SX}. Utilizing such decompositions, we propose a simple NN for emulating the solution to such problems using the $\tanh$ activation function together with different training objectives, such as residual or energy minimization. The key idea is to augment the approximation space with suitable exponential functions, similar to enriched spaces in finite element methods, e.g.~\cite{Kellogg}. Numerical examples in one and two dimensions, including a smooth non-tensor-product domain, illustrate the resulting parameter-robust behavior over the tested perturbation ranges.

We revisit the kinetic electrostatic electron nonlinear (KEEN) waves studied by Afeyan et. al in 2014 using a hybrid flow-mapping strategy that combines the characteristic mapping method (CMM) with numerical flow iteration (NuFi). The study extends the classical setup to dynamical ions, comparing their impact on the long-time KEEN dynamics with the static-ion case. To this end, we extend the CMM-NuFI framework with a multi-map strategy, assigning one map to the ion and one to the electron characteristic flow. The resulting problem exhibits a wide separation of spatial and temporal scales, driven by the fine structures generated by the ponderomotive force and by the large ion-to-electron mass ratio, which renders it computationally prohibitive for conventional grid-based methods. Our multi-map CMM-NuFI method efficiently resolves these disparate scales, enabling long-time simulations of KEEN dynamics with fully dynamical ions.

We analyze the numerical approximation of time-harmonic scattering by highly heterogeneous penetrable obstacles. These problems are especially challenging in the high-frequency regime, where the size of the scatterer $L$ is much larger than the wavelength, i.e., the wavenumber $k$ is such that $kL \gg 1$. Here, we further consider the situation where the scatterer contains different materials, with a characteristic size $\varepsilon$ such that $k\varepsilon \ll 1$. We propose a high-order multiscale finite element method, and provide an error analysis that is explicit in both $k$ and $\varepsilon$. Crucially, our error estimates suggest that using a high-order method should reduce the computational cost for large frequencies, which is corroborated by numerical examples.

We present a unifying framework to construct local Fortin operators for conforming mixed finite element pairs for the Stokes equations. The operators are constructed to satisfy the divergence-preservation property, local stability and approximation properties, and certain trace-preservation properties. For the latter, some of the finite element pairs require an enrichment of the velocity space. We present the construction for the $P_d-P_0$ element, the Bernardi-Raugel element, the conforming Crouzeix-Raviart element, a modified MINI element and generalised Taylor-Hood elements. Furthermore, we discuss implications of the existence of such a Fortin operator beyond inf-sup stability. These include applications to non-Newtonian fluid flow, problems with inhomogeneous Dirichlet boundary conditions, as well as uniform inf-sup stability relevant for domain approximation and moving domains.

Gaussian kernel sums are the computational core of maximum mean discrepancies (MMDs), kernel gradient flows, Stein variational gradient descent (SVGD), and many other kernel methods. At the same time, softmax attention has received an extraordinary amount of hardware-aware code engineering, culminating in flash attention. We show that Gauss kernel sums with arbitrary, signed weights can be evaluated via flash attention: two small input augmentations turn the normalized softmax reduction into the unnormalized Gauss sum, without writing a single line of custom GPU code. For feature dimension D>8 in fp16, this approach beats compiled PyTorch code as well as PyKeOps kernels (often significantly) in speed, memory-overhead and accuracy. Indeed, its memory scaling remains linear.

Operator splitting provides an attractive approach for the numerical integration of (coupled) port-Hamiltonian systems, as it allows the underlying system structure to be exploited at the level of the individual subproblems. However, the choice of the decomposition is not unique and may strongly affect both the computational efficiency and the preservation of the energy behavior of the original system. In this work, we investigate this interplay systematically and introduce energy consistency as a criterion for assessing splitting methods for port-Hamiltonian ordinary differential equations. We derive sufficient conditions under which a splitting based on a given decomposition inherits the energy behavior of the continuous system and use these conditions to analyze several decomposition strategies for coupled port-Hamiltonian systems. In particular, we compare decompositions that preserve the structure with approaches that exploit lower-dimensional subsystem dynamics or separated time scales. The analysis is complemented by numerical experiments using Strang splitting and its multiple-time-stepping extension. A scalable electro-thermal benchmark with fast electrical and slow thermal dynamics is employed to assess accuracy, energy behavior, and computational efficiency. The results demonstrate that preserving the port-Hamiltonian structure of the subflows is essential for energy-consistent splitting, whereas decompositions that exploit subsystem structure or time-scale separation can provide substantial computational advantages. In particular, the time-scale decomposition yields significant efficiency gains for systems with pronounced multirate characteristics, while structure-destroying decompositions may lead to undesirable energy behavior.

The classical unsolved Gauss circle problem concerns estimating the error between the number of lattice points inside a circle and the area of the circle as its radius tends to infinity. About a century ago, Hardy proposed a conjecture concerning this problem. In this paper, we attempt to provide numerical evidence in support of the Hardy conjecture through large-scale numerical computations.

We propose a novel framework for deriving semi-discrete discontinuous-Galerkin (DG) methods using operator semigroups for scalar conservation laws, then apply it to construct a class of high-order OFDG-type schemes [13] satisfying infinitely many local entropy inequalities with general E-fluxes on non-uniform meshes. Such schemes are further generalized to systems of conservation laws in any number of space dimensions by using entropy stable numerical fluxes in the sense of [1]. Finally, we prove optimal error estimates for smooth solutions to nonlinear scalar conservation laws, and prove strong convergence for discontinuous solutions to strictly convex conservation laws via compensated compactness.

Forsythe's conjecture, published in 1968, asserts that for each restart length $s$, every exact-arithmetic restarted conjugate-gradient iteration on a real symmetric positive definite problem either terminates or has normalised residuals that converge separately along the even and odd restart subsequences. Apart from the classical steepest-descent case, this asymptotic question remained unresolved in full generality for nearly six decades. We give a complete classification by restart length in the original finite-dimensional setting and identify a sharp threshold. For $s=2$ and $s=3$, every problem either terminates or has convergent even and odd residual directions. For every $s\ge4$, there is a diagonal positive definite counterexample of dimension $s+4$ which never terminates and whose even residual directions do not converge. Together with Akaike's theorem for $s=1$, this shows that the conjectured universal conclusion is true precisely for $s\in\{1,2,3\}$ and false for every $s\ge4$. The positive results follow from a degree-independent double-orthogonality identity and an analysis of the low-degree fixed-point sets. At restart length four, rational interval arithmetic and Sturm sequences certify a transverse Hopf point of the leading vector field of the rescaled squared-weight map. Analytic periodic-orbit and shadowing arguments yield the counterexample at restart length four, and degree elevation extends the construction to every larger restart length. The classification for all $s\ge2$ is also formally verified in Lean.

Aircraft outer mould lines are often heavily constrained early by payload requirements, manufacturability limits, and low-observability considerations. The remaining aerodynamic design space is typically small, highly constrained, and difficult to represent using a concise and physically meaningful set of design variables. A structural modal parameterisation method (MPM) is introduced for such problems. The method constructs a tunable pseudo-structure and solves its eigenvalue problem to obtain a modal basis, from which a selected subset of modes parameterises the geometry. Boundary conditions, stiffness distribution, density, thickness, and added masses are treated as intentional basis-design variables that shape the admissible deformation space prior to optimisation. The resulting design coordinates are independent of any specific flow solver, making the approach broadly applicable. For tightly constrained aerodynamic shape optimisation problems, it provides a practical balance of parameter conciseness, robust constraint awareness, and geometric flexibility.

I present a matrix-free symmetric interior penalty discontinuous Galerkin method for adaptively refined Cartesian meshes on GPUs. At non-matching interfaces, auxiliary shadow cells represent the adjacent coarse polynomial on the fine level. This approach converts non-matching interfaces into matching faces, permitting uniform face evaluation throughout the mesh. I prove that this construction is equivalent to the standard non-matching formulation. Furthermore, this representation yields a local geometric multigrid method in which frozen shadows provide inter-level boundary data and carry residual contributions to coarser levels. A primal–dual pairing eliminates shadow assembly from the Krylov iteration. Numerical experiments with cubic elements ($p=3$) on an NVIDIA A100 show stable multigrid convergence under increasing refinement depth and efficient GPU execution.

We develop an energy-based reduced-order model for micromagnetic magnetization dynamics that couples a convolutional autoencoder to a structured latent neural ordinary differential equation. Motivated by the precessional-dissipative structure of the Landau-Lifshitz-Gilbert equation, the latent vector field is generated from the gradient of a learned scalar potential through an antisymmetric operator and a symmetric positive-semidefinite dissipative operator. This potential is learned in nonunique latent coordinates and is not identified with the Gibbs free energy, but decreases monotonically along autonomous continuous-time solutions, while the antisymmetric component permits motion along its level sets. The encoder, decoder, latent energy, and operators are trained jointly on short trajectory windows using latent and decoded-rollout losses alone, without time-derivative supervision, physical-energy labels, or dissipation penalties. At inference, an initial state is encoded once, evolved in latent space, and decoded only at the requested output times, enabling substantially cheaper trajectory prediction than the micromagnetic solver used to generate the training data. We compare quadratic, deep, and additive deep-quadratic latent energies on two datasets parameterized by field amplitude and generated for the two applied-field directions of the NIST $\mu$MAG Standard Problem 4. Dissipative-only and antisymmetric-dissipative models achieve comparable accuracy on short training-style windows but differ substantially on uninterrupted rollouts, for which the antisymmetric-dissipative models provide markedly more accurate trajectory predictions. The deep-quadratic energy gives the best overall accuracy for both field directions and exhibits slower error growth when rollouts are extended to twice the training horizon.

In this manuscript, scalar norm-type accelerators are used for first time in order to increase the order of convergence of known iterative methods for solving the matrix equation $X^{-1}-A=0$, while decreasing the number of matrix-matrix products. Specifically, our proposed scheme upgrades the order of convergence in a $50\%$, achieving the sixth-order of convergence with the computational cost of fourth-order Hyperpower method, in terms of matrix-matrix multiplications. In the main results, not only the convergence and its order is proven, but also its efficiency and stability. This new technique provides good numerical results, even compared with existing fourth-, sixth- and eighth-order schemes in large matrices of up to $10^6$ entries, but also promises the opening of a new kind of procedures with good performance and scalability. Finally, an application to digital image restoration is made to show the behavior of the method, compared with existing ones.

The $\ell_p$-Lewis weights of a matrix are defined by a fixed-point equation. For $p<4$, Cohen and Peng [CP15] showed that iterating an equivalent rearrangement of this equation computes Lewis weights to high precision; for $p\geq4$, prior high-precision methods instead use optimization-based approaches. We show that the direct Lewis fixed-point iteration, appearing in the thesis of Lee [Lee16], computes Lewis weights to high precision for every $p>2$. For a matrix $\mathbf{A}\in\mathbb{R}^{m\times n}$ partitioned into row blocks $\mathbf{A}_{[1]},\ldots,\mathbf{A}_{[k]}$, we compute, for $p>2$, coordinatewise $\varepsilon$-approximate $\ell_p$ block Lewis weights in $O\left(p\log\frac{p\sqrt{\sum_{i=1}^k\operatorname{rank}(\mathbf{A}_{[i]})}}{\varepsilon}\right)$ rounds of exact leverage-score-vector computations. For ordinary Lewis weights, this becomes $O\left(p\log\frac{p\sqrt{m}}{\varepsilon}\right)$, improving the $O\left(p^2\log(m/\varepsilon)\right)$ bound of Gribling, Sidford, and Zhang [GSZ26] for $p\geq4$. Our main observation is that each direct Lewis update contracts the KL divergence to the true weights by a factor of $1-\frac{2}{p}$. We also give an alternate explanation of this contraction through volume sampling and entropic independence. Synthetic experiments closely match our predicted local contraction rates and iteration counts grow approximately linearly with $p$; real-data experiments illustrate the information-concentration tradeoff of finite-$p$ block Lewis designs.

We study probabilistic numerical methods for solving nonlinear PDEs from a Bayesian nonparametric perspective. Given noisy evaluations at random collocation points, we place a truncated Gaussian series prior on the unknown solution and establish contraction at the minimax nonparametric rate, up to a logarithmic factor. Our main results give Gaussian approximations of the posterior in positive-order Sobolev spaces and, under suitable conditions, in the uniform topology. This contrasts with classical ill-posed inverse problems, where Bernstein–von Mises theorems typically require substantially weaker topologies. Here, the observation operator is differential rather than smoothing, and inversion of its linearisation gains regularity, making these strong-topology results possible. The posterior may be centred at either the posterior mean or the posterior mode. We further prove that the Gaussian Laplace approximation is asymptotically equivalent to the true posterior at a $\sqrt{N}$-scale.

What are the transversality conditions for an optimal control problem when the boundary conditions are defined by differential equations? This seemingly bizarre question is motivated by trajectory optimization problems in the $N$-body system. The question, however, is more fundamental and goes beyond problems in astrodynamics to nonintegrable dynamical systems in general. The main contribution of this paper is the development of generic initial- and final-time transversality conditions for optimal control problems whose boundary conditions are defined in terms of differential equations with side conditions. The mathematical definition of differential boundary conditions are part of the foundations developed in this paper. To support the new fundamentals, the concept of coordinated/uncoordinated clock times and weak adjoint covectors are introduced. In the case of uncoordinated clock times, the new transversality conditions reveal that there exists a special situation where a weak adjoint covector is orthogonal to the vector field of the boundary differential equation. This condition is sharply different from the classical statement of orthogonality with respect to the endpoint manifold. The theorems developed in this paper are generic. An application of the theorems to several cases in the three-body problem are described in separate papers.

We develop a convergence theory for geometric multigrid with vertex-patch smoothers applied to cut finite element discretizations of the Poisson problem. The framework addresses non-inherited level forms and the mismatch between the physical and active domains. Using the discrete extension property, we prove two-level convergence bounds uniform in the mesh size and the cut geometry, and W-cycle bounds under an additional smallness assumption on the two-level rate. The numerical experiments intentionally use the stronger V-cycle, for which no convergence bound is claimed here. The convergence constants degrade with the degree $p$. Lowering the ghost penalty improves iteration counts. An aligned two-cell model exhibits a semidefiniteness threshold of order $p^{-2}$, whereas the visibility scale of a degree-$p$ cut mode decreases exponentially. Experiments at the model threshold reduce the iteration counts, but do not establish an assembled-operator threshold.

We partially answer an open problem, posed by Nick Higham, concerning the conditioning of solutions to Sylvester and Lyapunov equations. The question arises in the backward stability analysis of numerical algorithms for these equations. We first show that the solution to the Sylvester equation $AX-XB = C$ can be arbitrarily ill-conditioned even if $A, B, C$ and the Kronecker sum $I \otimes A - B^T \otimes I$ are all perfectly conditioned. We then derive general a priori bounds on the condition number of the solution, as well as bounds for the Sylvester equation when $A$ and $B$ are diagonalizable. We also provide lower bounds involving matrix exponentials and Zolotarev numbers. For the Lyapunov equation $AX+XA^T = -C$, we obtain upper bounds in two settings: (i) when $A$ is symmetric positive definite while $C$ is symmetric negative definite, and (ii) when $A$ is strictly dissipative and $C$ is symmetric positive definite.

The geometric morphology of deposited filaments can significantly influence the structural performance and stability of 3D concrete-printed (3DCP) structures. However, most finite element (FEM)-based approaches for buildability assessment represent printed layers as simplified rectangles, potentially limiting predictive accuracy. This study proposes a geometry-informed modelling framework that integrates the deep-learning-based filament shape prediction tool ShapeGen3DCP with a layer-activation FEM approach to investigate the effect of realistic filament geometries on buildability. The framework generates geometry-aware numerical models directly from material and process parameters, eliminating the need for experimental filament characterization or computationally intensive fluid-flow simulations. Validation against experimental data and a parametric study of rectilinear walls demonstrate that extrusion parameters and the resulting filament geometry can significantly influence buildability predictions. Realistic filament representations are particularly important for free-flow deposition, whereas layer-pressing strategies are less sensitive to geometric simplifications. Among the investigated representations, an elliptical approximation provides an effective balance between geometric fidelity and modelling simplicity. When rectangular representations are preferred to enable regular computational meshes for faster simulations, defining their dimensions based on volume conservation improves prediction reliability compared with calibrating them using either the maximum filament width or the interlayer contact width. Overall, the proposed methodology demonstrates the importance of incorporating filament geometry into 3DCP simulations and provides practical guidance for selecting efficient and accurate geometric representations for buildability assessment.

Transport models underpin the design and optimisation of chemical processes, yet their constitutive laws are typically hand-picked, their experimental initial conditions are assumed exact, and their computational cost restricts optimisation to a handful of process variables. Here we introduce a differentiable hybrid modelling framework that embeds neural networks within a GPU-accelerated population balance equation solver to discover constitutive laws and assimilate initial conditions directly from experimental crystallisation data. Using thirteen real experiments with multi-dimensional size and shape measurements, we show that treating measured initial conditions as exact yields a physically inconsistent rate law despite improved predictive accuracy, whereas assimilating them recovers a physically consistent, non-linear law and reduces trajectory error sevenfold relative to conventional modelling. End-to-end differentiability further enables gradient-based optimisation of high-dimensional operating profiles, identifying a 300-variable temperature profile that reduces crystal aspect ratio by 25% at constant mass, facilitating superior product attributes. Our framework augments, rather than replaces, mechanistic modelling, offering a general route to learning transport models from data and optimising processes with them.

On triangular meshes, we analyze local mesh refinement based on longest-edge bisection equipped with the serial longest-edge propagation-path closure. Ties are resolved by terminal priority: if the incoming shared edge is a longest edge of the neighboring triangle, the pair is declared terminal and that edge is bisected. For every adaptive mesh sequence $\mathcal{T}_0, \mathcal{T}_1, \ldots, \mathcal{T}_L$ with a sequence of marked sets $\mathcal{M}_0, \mathcal{M}_1, \ldots, \mathcal{M}_{L-1}$, we prove the cumulative closure estimate $$\#\mathcal{T}_L-\#\mathcal{T}_0 \lesssim\sum_{\ell=0}^{L-1}\#\mathcal{M}_\ell.$$ The proof derives single-mark locality from the uniform multiplicative gap between possible descendant diameters implied by finite similarity classes, and converts this locality into the cumulative estimate through a Binev–Dahmen–DeVore-type charging argument.

We study a time fractional convection-diffusion-reaction equation in a bounded domain $\Omega\subset\mathbb{R}^2$. A stabilized numerical scheme satisfying a discrete maximum principle is constructed by combining the conforming linear finite element method with the algebraic flux correction method. The resulting semi-discrete scheme is nonlinear, and its well-posedness is established. Assuming nonsmooth initial data, we derive error estimates for the semi-discrete scheme using energy arguments. For the temporal discretization, we employ the L1 method, obtaining a fully discrete scheme for which we prove well-posedness and the discrete maximum principle. We also present numerical experiments that validate the order of convergence as well as we test our schemes to solutions that possess layers.

This paper addresses the numerical solution of nearly incompressible linear elasticity boundary value problems and their associated eigenvalue problems. A mixed finite element formulation based on Mini element is proposed to circumvent the locking phenomenon that plagues standard low-order elements in the nearly incompressible limit. For the boundary value problem, we establish the well-posedness of mixed variational formulation and derive a priori error estimates that are uniform with respect to the Lamé constant $\underline{\lambda}$, thereby proving the method's locking-free property. For the eigenvalue problem, we develop an efficient non-nested augmented subspace algorithm designed within the mixed finite element framework. A comprehensive convergence analysis is provided for the discrete eigenvalue approximation and the proposed iterative solver, demonstrating that the convergence rates remain independent of $\underline{\lambda}$. Numerical experiments on both problems confirm the theoretical conclusions, showing optimal convergence rates and robustness as $\underline{\lambda} \to \infty$. The results validate the effectiveness of Mini element and proposed augmented subspace algorithm for reliable and efficient computation in the nearly incompressible regime.

In our study, Darbo's fixed point theorem and the coupled fixed point theorem have been extended. The two theorems were generalized using $C-class$ mapping. Using the Darbo type theorem, we gave an existence result for a functional integral equation along with an appropriate illustration. In addition, we explored the numerical solution of a given functional integral problem using modified homotopy perturbation.

We consider finite-dimensional nonlinear inverse problems arising from finite element discretizations of elliptic inverse coefficient problems such as the Calderón problem with finitely many measurements and unknowns. Such inverse coefficient problems are notorious for their nonlinearity and ill-posedness, and numerical solvers tend to depend strongly on good initial values. In this work, we develop a new locally convergent algorithm with an explicit residual criterion that ensures convergence to the inverse problem solution, and a globalized variant that is guaranteed to automatically switch to the faster locally convergent algorithm after finitely many global search steps.

Rigorous results show that feedforward neural networks can overcome the curse of dimensionality in the numerical approximation of high-dimensional partial differential equations (PDEs), but comparatively little is known about residual neural networks (ResNets) in the nonlinear PDE setting. We prove that ResNets overcome the curse of dimensionality in the numerical approximation of solutions of semilinear heat equations with globally Lipschitz continuous, gradient-independent nonlinearities: under polynomial growth and network approximability hypotheses on the PDE data, there exist $\eta\in(0,\infty)$ and ResNets $\Psi_{d,\varepsilon}$, $d\in\mathbb{N}$, $\varepsilon\in(0,1]$, with at most $\eta d^{\eta}\varepsilon^{-\eta}$ parameters whose realizations approximate the solution in dimension $d$ with an $L^2$-error of at most $\varepsilon$. The proof represents one deterministic realization of a multilevel Picard estimator by a ResNet whose shortcut connections transmit the spatial variable and a scalar accumulator, while the residual branches successively add the summands of the estimator. For ridge-sum initial conditions, admissible sigmoidal activations, and globally Lipschitz truncations of the nonlinearity, we obtain, for every $\xi>0$, the explicit bound $C_\xi d^{4+\xi}\varepsilon^{-(3+\xi)}$ on the number of parameters.

A common approach to representing disk-like and sphere-like geometries in isogeometric analysis is to use collapsed-edge singularities in the underlying parameterization. The resulting tensor-product B-spline spaces on such singular parameterizations lack the required smoothness to be used in isogeometric analysis. Polar splines (Toshniwal et al., CMAME 2017) rectify this issue by extracting a smoother subspace, and have been successfully used to discretize high-order partial differential equations, as well as to perform structure-preserving discretizations of the de Rham complex.The latter is particularly relevant for mixed formulations that appear in applications such as electromagnetism and fluid dynamics. The main contributions of this paper are twofold. We combine the polar spline construction with hierarchical splines to obtain adaptively-refinable polar splines. Moreover, we do so in the context of structure-preserving methods. That is, we show how to construct hierarchical polar spline $k$-form spaces, prove that the corresponding basis functions are linearly independent, and prove that they form a cohomologically-correct discrete de Rham complex. These contributions provide a mathematically sound foundation for adaptive structure-preserving simulations on polar geometries. We also provide numerical experiments that validate the theory and illustrate the practical behavior of the construction with problems where preserving the structure is mandatory. These include a study of spurious harmonics that highlights how cohomology-breaking refinements differ between the polar and tensor-product settings, a discretization of a Maxwell eigenvalue problem on a locally-refined polar mesh, and an adaptive convergence study for a Hodge–Laplace problem.

We study the semi-discrete quadratic Wasserstein energy. The energy is nonsmooth at collisions of sites. We prove local Lipschitz continuity on the full configuration space, together with global semiconcavity, coercivity, and dissipativity; show that every global minimizer is interior and collision free; and establish $C^2$ regularity on the collision-free configuration space. The gradient is expressed through the barycenters of the balanced Laguerre cells, while the Hessian is given by an explicit facet formula and satisfies a global one-sided bound. We also solve the one-dimensional problem explicitly in each ordering chamber and give a two-site example on the unit square with non-minimizing Lloyd fixed points. For $d\ge 2$, we then formulate an entropy-regularized relaxed control of the Langevin temperature. The controlled dynamics is strongly well posed, nonexplosive, and collision free. Its value function is a classical interior solution of the exploratory Hamilton-Jacobi-Bellman equation; the Laplacian of the value function is locally $C^1$, which yields a locally Lipschitz optimal temperature feedback. Independently of this optimal-control result, for every fixed Borel temperature rule bounded away from zero, and every sufficiently small step size, the associated Gaussian Euler chain is geometrically ergodic with a full-support invariant law. The raw iterates do not converge, whereas the best-so-far energy converges almost surely to the global minimum and the running record approaches the set of global minimizers.

We first prove spectral convergence of the random feature method (RFM) for multidimensional targets in Sobolev, Gevrey, ultra-analytic, and bandlimited classes. The analysis establishes general high-probability approximation estimates in the interpolation scale generated by a kernel integral operator. On a single event determined only by the sampled features, one random space approximates every target in a prescribed source ball; moreover, for each target, a single coefficient vector defines an approximant that attains spectral accuracy simultaneously in all admissible error norms. For both regularity-adapted frequency distributions and uniform distributions on growing frequency windows, the resulting rates range from super-exponential to algebraic, depending on the regularity of the target. Second, we establish abstract error estimates for strong- and weak-form RFM discretizations, thereby converting the preceding approximation bounds into convergence estimates for multidimensional second-order elliptic boundary value and eigenvalue problems. Finally, for random feature matrices (RFMtxs), we prove super-exponential singular-value decay with Fourier features and exponential decay with $\tanh$ features, together with corresponding condition-number lower bounds. The analysis identifies a common mechanism: the same spectral approximation that yields high accuracy also drives severe ill-conditioning.