Papers

In this paper, we introduce a deep truncated forward-backward stochastic differential equation (FBSDE) method for high-dimensional partial differential equations (PDEs). Compared with existing deep-learning solvers for fully coupled FBSDEs, where strong coupling may lead to numerical instability, our approach exhibits improved stability. The proposed method combines gradient-truncated iterative decoupling with fictitious-play averaging to separate the forward and backward processes in a coupled framework. This preserves the coupled dynamics while reducing the unstable feedback induced by parameter-dependent forward paths during optimization. Furthermore, we incorporate a pathwise consistency term to create explicit local gradient shortcuts, thereby providing a structural mechanism that may mitigate gradient vanishing. We also derive a residual-based error estimate and establish conditional convergence of the fully discrete numerical approximations under suitable conditions, in which the pathwise consistency loss is not required. Our approach is particularly effective for convection-dominated equations, where the coupled formulation provides a stable representation of nonlinear transport without introducing singular terms into the BSDE. Numerical experiments demonstrate improved accuracy and stability in both low- and high-dimensional problems and robust performance for strongly coupled problems.

The geometrically decomposed bases of the polynomial differential form spaces $\mathcal{P}_r\Lambda^1$ and $\mathcal{P}_r^-\Lambda^1$ on a simplex are products of a barycentric scalar polynomial and a directional or Whitney 1-form. These bases depend on an ordering of the simplex vertices. On unstructured meshes, neighbouring cells need not agree on this ordering, making it non-trivial to achieve conformity in finite element software. Existing strategies compute degree-of-freedom transformation matrices numerically, or impose a global vertex ordering through face-bubble spaces or by preprocessing the mesh. We consider the standard basis of the $\mathcal{P}_r^-\Lambda^1$ (trimmed) space and introduce a new basis for the $\mathcal{P}_r\Lambda^1$ (full) space, for arbitrary polynomial order $r \geq 1$ and dimension $D \geq 2$, and prescribe their basis polynomials as shape functions. We show that the pullback of shape functions under the change of coordinates induced by an arbitrary vertex relabelling $\pi$ in the symmetric group $S_{D+1}$ admits a closed-form, combinatorial expression. For most basis functions this pullback is a single basis function of the relabelled ordering, up to sign. On an explicitly characterised "filter-hit" set, this pullback is a signed sum of at most $D$ (full space) or exactly two (trimmed space) relabelled basis functions. All coefficients are in $\{-1, +1\}$ for all $r$, $D$ and $\pi$. The inverse transformation is obtained by computing the formulas at $\pi^{-1}$, so no numerical inversion of a change-of-basis matrix is needed. Conforming assembly and evaluation on simplicial meshes with arbitrary vertex orderings reduce to index manipulation. Our results are compared with the study of relabelling-invariant bases of Berchenko-Kogan and Licht, and verified in an open-source Julia implementation in the Gridap.jl library.

A mesh-free numerical framework based on Kolmogorov–Arnold Physics-Informed Neural Networks (KAN-PINNs) is developed for the approximation of fourth-order elliptic boundary value problems, with specific application to the biharmonic equation governing thin plate deflection. Unlike conventional Multi-Layer Perceptrons relying on fixed nodal activations, learnable univariate functions parameterized via radial basis functions are deployed on network edges, while high-order differentiability is preserved through hyperbolic tangent activation mappings. The severe numerical stiffness inherent to fourth-order differential operators and fully clamped boundary conditions is addressed through a direct normalized residual formulation coupled with a hybrid, multi-stage AdamW-to-L-BFGS optimization pipeline. An automated 24/7 hill-climbing search protocol is implemented to systematically calibrate boundary penalty weights and optimization schedules. When evaluated on a smooth manufactured benchmark on the unit square, an error reduction factor exceeding $150\times$ is attained over successive iterations, culminating in a final relative $L_2$ error of $1.593 \times 10^{-5}$ ($0.00159\%$) and a training loss of $3.065 \times 10^{-6}$ utilizing only $7,801$ trainable parameters. These results demonstrate that high-order PDE problems can be accurately resolved using compact, mesh-free KAN-PINNs architecture without requiring auxiliary variable transformations or discrete mesh generation.

While most attention logits can be computed in low precision without degrading numerical stability, current attention kernels fail to exploit this phenomenon. We introduce a novel hardware-algorithm co-design in the form of mixed-precision FlashAttention. Our method accumulates key-query products and evaluates their exponentials in 8-bit formats, then adaptively identifies sensitive sub-blocks and recomputes them in 16-bit formats. We propose the specifications for a dedicated accelerator capable of executing this pipeline efficiently. Simulated experiments with Qwen3 and Gemma 3 show that rerouting a selective minority of sub-blocks to high precision is sufficient to recover the baseline model performance.

Collocation methods based on Birkhoff interpolation at Gaussian-type points are well-conditioned for one-dimensional initial and boundary value problems [Wang et al., SIAM J. Sci. Comput. 36 (2014)], but the underlying construction does not extend directly to multiple dimensions. We address the long-standing ill-conditioning of multidimensional collocation methods for second-order elliptic-type problems. The key observation is that the second-order differentiation matrix and its inverse, the pseudospectral integration matrix (PSIM) constructed from Birkhoff interpolation at Legendre-Gauss-Lobatto points, are both similar to symmetric negative definite matrices. This enables stable diagonalisation of the dense, non-symmetric and ill-conditioned differentiation and integration matrices, even for thousands of collocation points, and leads to efficient multidimensional Birkhoff preconditioners. For variable-coefficient problems, the coefficients are incorporated directly into the diagonalisation and preconditioner construction, which is essential for highly anisotropic, high-contrast, oscillatory and degenerate elliptic operators. We provide spectral analysis and extensive two- and three-dimensional numerical experiments, demonstrating substantial reductions in condition numbers and nearly polynomial-degree-independent GMRES convergence while retaining high-order accuracy. The resulting Birkhoff-collocation schemes make multidimensional spectral collocation methods practical for challenging elliptic problems.

We focus on solving the Lyapunov equation $AX + XA^T = F$, where $A$, $X$ and $F$ are square matrices, $A$ is symmetric positive definite (SPD) and sparse, and $F$ is symmetric and of a low-rank. The solution $X$ is then also symmetric and in general dense, but it can be approximated by a low-rank matrix. If $n$ is large, $X$ cannot be computed directly, but it is accessible by using the so-called low-rank arithmetics (LRA). The equation is solved by matrix reformulation of the method of conjugate gradients (MCG). We are interested in the behavior of ranks of the solution approximations $X_\ell$, residuals $R_\ell$, and direction vectors $P_\ell$ during iterations.

We consider harmonic approximation on the torus $\mathbb{T}^d$, the sphere $\mathcal{S}^d$, and the rotation group $\mathrm{SO}(3)$. While several global approaches are available for computing a harmonic expansion from scattered data, their stability and, especially, their runtime depend strongly on the geometry of the nodes, in particular on a sufficiently small fill distance. Even a single large hole in the data may cause instability and long runtimes. We propose harmonic approximation via moving least squares (HAMLS), a global harmonic approximation scheme that avoids the ill-conditioned global solve and works well in both the overdetermined and the underdetermined settings. The main idea is an intermediate transition from the scattered data to a quadrature grid, which is realized via moving least squares (MLS). This replaces one large global problem by many small local ones that can be regularized individually. The global harmonic approximation is then obtained via quadrature using a fast Fourier transform. This makes the global stage fast, stable and non-iterative. For sufficiently smooth functions sampled at well-distributed nodes with small fill distance $h$, we bound the resulting $\mathrm{L}^2$-error by the error of the harmonic approximation obtained from exact values on the quadrature grid, plus a term of order $h^{K+1}$, where $K$ is the polynomial degree employed in the MLS step. Numerical experiments on $\mathcal{S}^2$ with ill-distributed nodes show that HAMLS achieves errors comparable to those of global least-squares approximation, while being significantly faster, especially for larger bandwidths.

In this paper, we introduce a new set of $\mathbb{L}_2$-orthogonal polynomials on the hemisphere $\mathbb{S}^d_+:=\{\mathbf{x}\in\mathbb{R}^{d+1}:\|\mathbf{x}\|=1, \mathbf{x}\cdot\mathbf{e}_{d+1}\geq0\}$ and establish $\mathbb{L}_1$ Marcinkiewicz-Zygmund inequalities on the hemisphere for even and odd spherical polynomials. Hemispherical $t$-designs provide equal weight cubature rules on the hemisphere which are exact for polynomials up to degree $t$. We propose variational characterizations of hemispherical $t$-designs by the $\mathbb{L}_2$-orthogonal polynomials. Moreover, based on the Marcinkiewicz-Zygmund inequalities and the $\mathbb{L}_2$-orthogonal polynomials, we prove that for each $n\geq C_d t^d$, there exists a hemispherical $t$-design on the hemisphere with $n$ points, where $C_d$ is a constant depending only on $d$.

Nonlinear partial differential equations (PDEs) underpin computational fluid dynamics, yet resolving nonlinear transport on fine grids remains computationally demanding. Quantum linear-evolution algorithms offer a possible route to large-scale simulation but cannot directly propagate nonlinear coupling. Here we develop a Frozen–Oseen quantum homotopy method (FOQHAM) framework that links transport-aware auxiliary-operator selection to the size of the resulting linear representation. We freeze the full Fréchet derivative at a prescribed flow profile, retaining transport and profile-gradient coupling and close each prescribed finite-order homotopy hierarchy as a fixed affine linear system by product lifting. We formulate its propagation through Linear Combinations of Hamiltonian Simulations (LCHS), without outer homotopy iterations or profile updates. For spatially semidiscrete equations, we establish local approximation bounds and exact closure of the finite-order hierarchy. Classical tests on Burgers, Korteweg–de Vries, Zakharov–Kuznetsov, and one- and two-dimensional isothermal compressible Navier–Stokes equations provide numerical evidence of accurate non-iterative approximations in the tested regimes. The framework thus connects nonlinear PDE approximation to quantum linear evolution and offers a conditional path toward quantum acceleration.

Accurate long-time prediction of Cahn-Hilliard dynamics, particularly in the late-stage coarsening regime, is computationally demanding because it requires simulations over extended physical time. Various machine learning approaches have been developed to speed up the simulations through efficient surrogate evaluations. While the learned map can be accurate over a single learned interval, small prediction error will accumulate when the map is repeatedly applied in a long-time prediction, leading to large numerical errors. In this work, we identify the structure function for the coarsening dynamics of the Cahn-Hilliard equation and propose a two-stage structure-preserving Fourier Neural Operator (FNO) learning framework. Stage~I introduces a bound-conforming FNO (bcFNO) with a soft, magnitude-dependent amplitude penalty. Stage~II retains the bound-conforming objective and supplements it with structure-function regularization, yielding a structure-preserving FNO (spFNO) that corrects late-stage coarsening statistics. Numerical experiments on the various spatial resolutions and reference-solver time steps suggest that the bcFNO stabilizes the long-time prediction by suppressing large amplitude excursions, without changing the low online cost of the FNO. Furthermore, the structure-function term consistently reduces discrepancies in the normalized structure function and provides further correction in associated late-stage coarsening quantities, including the $L^3(t)$ growth trend and the scaling collapse. These results indicate that our multi-stage implementation of physical principles in operator learning could be applied to other multiscale systems with statistical scaling, such as the functionalized Cahn–Hilliard equations and turbulence, provided that the state constraints and statistical observables are adapted to the governing dynamics.

We study how well deep neural networks approximate and learn spectral Barron functions. Recent studies have shown that these function classes can be efficiently approximated by shallow neural networks without suffering from the curse of dimensionality. We complement these results by providing new approximation bounds for deep networks with ReLU activation and establishing the minimax rates for learning these function classes. Specifically, we show that $d$-dimensional spectral Barron functions with smoothness index $s>0$ can be approximated by deep ReLU neural networks with approximation rate $\widetilde{\mathcal{O}} (S^{-\frac{1}{2}-\frac{s}{d}})$, where $S$ denotes the number of nonzero parameters in the network. Using this approximation result, we further show that deep ReLU neural networks can learn spectral Barron functions in a fast rate $n^{-\frac{d+2s}{2d+2s}}$ with $n$ training samples. Finally, we prove that this convergence rate is minimax optimal up to logarithmic factors.

We develop a positive-conservative Fourier optimization (PCFO) method for computing ground states of Bose–Einstein condensates with higher-order interactions. The ground-state problem admits a convex density formulation, but the singular behavior near zero density poses difficulties for high-accuracy computation. We introduce convex regularizations of the density energy and establish their $\Gamma$-convergence on the computational domain as the regularization parameters vanish. The regularized problem is discretized by a Fourier pseudospectral method with mass conservation and nodal positivity. For the resulting constrained optimization problem, we use an accelerated first-order method for the main energy reduction, followed by a Newton refinement to reduce the remaining optimality residual. Numerical experiments show spectral-type convergence for regularized problems, quantify the effects of the regularization parameters on accuracy and spatial resolution, and demonstrate the effectiveness of the two-stage optimization method.

Nearly incompressible elasticity is often discretized by mixed methods to avoid locking, but the resulting saddle-point systems can be expensive to solve. We consider a weakly symmetric multipoint stress discretization whose stress and rotation unknowns can be eliminated independently over vertex interaction regions, leaving a symmetric positive definite system for cell-centered displacement only. Although smaller, the reduced system contains a parameter-dependent Schur complement that becomes difficult to precondition as the Lamé ratio increases. We show that its energy consists of a uniformly coercive shear part and a dominant semidefinite volumetric part. Standard cell-centered interpolation need not preserve the volumetric kernel and can introduce energy amplified near incompressibility. Motivated by this structure, we develop a two-grid preconditioner combining two complementary coarse subspaces: a conventional displacement subspace and a discrete-curl subspace lying exactly in the fine-grid volumetric kernel. Symmetric vertex-patch smoothing completes the displacement-only solver. The energy splitting and kernel-compatible decompositions yield a condition-number bound independent of the Lamé ratio for each fixed grid pair. Two- and three-dimensional experiments confirm the approximation properties of the locally eliminated discretization and show that the kernel correction prevents the deterioration of the displacement coarse subspace alone, also for discontinuous material coefficients.

The fixed-point fast sweeping methods with weighted essentially non-oscillatory (WENO) local solvers are a class of efficient and high-order accuracy numerical methods for solving steady-state solutions of hyperbolic conservation laws. However, with the classical WENO-JS local solver, the iteration residue of high-order fixed-point fast sweeping scheme often has difficulty to settle down to round-off errors. To achieve the full convergence in a fast sweeping method, the fixed-point fast sweeping methods with non-traditional WENO local solvers based on unequal-sized substencils were designed. However, the WENO schemes based on unequal-sized stencils are more complex and in general more expensive in computational costs than the classical WENO-JS schemes. In this paper, we go back to the classical WENO-JS local solver and develop a new fully convergent fifth-order fixed-point fast sweeping method for solving steady-state problems of hyperbolic conservation laws. Based on recent studies on the nonlinear weighting process around discontinuities of solution, we apply the technique of frozen weights for avoiding unnecessary adjustment of nonlinear weights in the WENO-JS local solver, which freezes the nonlinear weights once the residual sequence in the fast sweeping iterations has stabilized. Different from the existing work on frozen weights, we design a simple and robust approach to judge stabilization of iteration residues and determine the iteration step when the nonlinear weights are frozen in the fast sweeping method. Extensive numerical experiments on a wide range of challenging two-dimensional steady-state problems demonstrate that, unlike the previous fast sweeping method with the fifth-order WENO-JS local solver, the proposed new scheme consistently drives the iteration residues to round-off errors and achieves the full convergence.

We develop a Halpern-accelerated majorized alternating direction method of multipliers (ADMM) with possibly indefinite proximal terms for linearly constrained convex composite optimization. Majorization simplifies the subproblems but introduces a forward perturbation to the standard degenerate proximal point representation. We reformulate majorized ADMM as a perturbed degenerate proximal point method and characterize the fixed points of the induced mapping as Karush–Kuhn–Tucker (KKT) solutions. Under range compatibility and metric cocoercivity conditions, the mapping is nonexpansive in the seminorm induced by a possibly singular preconditioner. In finite dimensions, a verifiable block-matrix condition permits indefinite proximal regularization while retaining the required metric properties. A relaxed Halpern iteration then yields an $\mathcal O(1/k)$ fixed-point residual bound, sharp in the general nonexpansive setting, and an $\mathcal O(1/k)$ nonergodic KKT residual bound at the intermediate iterates. Numerical experiments illustrate the predicted residual decay and the benefits of majorization and indefinite proximal terms.

In this article, the problem of completing a tensor from some incomplete knowledge of its entries is treated by adopting a worst-case perspective, given the realistic assumption that the tensor's low-order ANOVA terms are dominant. We survey and leverage some recent all-purpose results from the field of Optimal Recovery to provide solutions on a theoretical level. But the accompanying constructions of optimal completion procedures, which often feature semidefinite programs, are not directly applicable in the tensor case due to the huge dimensions involved. To resolve the issue, we put forward a storage-friendly way to produce low-order ANOVA projections based on the fast Fourier transform (FFT), while exploiting the specificities of the completion problem to efficiently compute regularizers and extremal eigenvalues. Numerical experiments on synthetic tensors and real-world datasets demonstrate the accuracy and scalability of our FFT-based method.

We propose ABrA-GD (Adaptive Bregman Accelerated Gradient Descent), an adaptive algorithm for relatively smooth convex optimization. The algorithm is derived from a computable primal–dual certificate that guarantees progress by comparing the objective value with a lower bound on the minimum of a regularized objective. This certificate enables adaptation to both smoothness and geometry, using only the relative strong convexity constant as a problem-dependent input. We introduce the dual Bregman Length Distortion Factor (BLDF), which measures how dual Bregman lengths change under an anchor shift or rescaling. Under bounded local dual BLDF, ABrA-GD achieves an accelerated $O(1/k^2)$ rate for convex objectives and an accelerated linear rate for relatively strongly convex objectives.

We study the problem of sampling from a continuous density $\pi\propto \exp(-V)$ on $\mathbb R^d$, where $V\in C^2(\mathbb R^d)$ has a $\beta$-Lipschitz gradient and $\pi$ satisfies a logarithmic Sobolev inequality with constant $\alpha^{-1}$, and write $\kappa = \beta/\alpha$. We introduce the Gaussian cloud sampler, which achieves total variation accuracy $\varepsilon$ using $\widetilde O(\kappa d^{1/5}\,\text{polylog}(1/\varepsilon))$ gradient queries in expectation. The algorithm uses first-order rejection sampling (FORS) to correct the law of smoothed Picard HMC trajectories. To do so, we represent the iterates of the ideal Picard iteration via Gaussian clouds, whose centers are never evaluated, and we develop a suite of likelihood correction gadgets which only use samples from this indirect cloud representation.

We present a novel algorithm for calculating the preimage of an affine space through a product ${J}={J}_{T-1}\cdots{J}_0$ of matrices ${J}_t$ of special form: finding the largest input space $\mathbf{X}$ such that $\mathbf{x}\in\mathbf{X}$ implies ${J}\mathbf{x}\in \mathbf{Y}$, where $\mathbf{Y}$ is a given output affine space. These special matrices arise in AD, where the Jacobians $J$ describing the linearized computation have precisely this structure: the product of a series of linearized primitive numeric operations. This allows us to use the new algorithm to formulate Null-A mode preimage AD, which finds the affine preimage through the Jacobian or Jacobian transpose of a numeric computation. This is a generalization of the inverse AD problem of solving ${J}\acute{\mathbf{x}}^{\ast}=\acute{\mathbf{y}}^{\ast}$ or ${{J}}^{T}\grave{\mathbf{y}}^{\ast}=\grave{\mathbf{x}}^{\ast}$. The key is to represent affine spaces in a fashion which lends itself to efficient preimage calculation, in a compositional and quasi-local fashion, through a succession of matrices ${J}_t$. Unlike previous methods, Null-A preimage mode AD allows the ${J}_t$ matrices to be non-square, corresponding to a computer program whose number of active variables swells and shrinks during the computation. When ${J}$ is square and the initial affine space is a single point, this finds the conventional inverse. But in the more general case, having the entire affine space provides freedom which can be leveraged in a problem-specific manner. We apply the method to small problems on-CPU where the $J_t$ are linearized scalar unary or binary numeric functions; and to larger problems on-GPU where the $J_t$ are linearized aggregate array operations like convolution and attention.

We study transient solute dispersion in pressure-driven micropolar flow through a concentric annulus representing a catheterized artery. The catheter is impermeable to solute, whereas the arterial wall combines irreversible removal with reversible surface retention. The exact steady micropolar velocity field is coupled to the Gill-Sankarasubramanian generalized-dispersion framework; the resulting bulk-surface hierarchy is advanced with Rannacher-damped Crank-Nicolson time stepping to determine $K_0(t)$, $K_1(t)$ and $K_2(t)$. The exchange coefficient is independent of the micropolar parameters and satisfies $-K_0(0^+)=2(\beta+\theta Da)/(1-\lambda^2)$, whereas $K_2-Pe^{-2}\sim\sigma_v^2t$ initially, so the leading short-time dispersion is hydrodynamic and independent of wall kinetics. In the weak-reaction regime, micropolarity modifies convection approximately linearly and shear dispersion quadratically through the velocity-amplitude factor. A narrow-gap analysis with $\varepsilon=1-\lambda$ gives $\bar v\sim(2-N_c)\varepsilon^2/6$ and $K_2-Pe^{-2}\sim(2-N_c)^2\varepsilon^6/7560$ for a fixed pressure gradient, revealing a sixth-power suppression of shear dispersion as the catheter approaches the arterial wall. Increasing $\lambda$ from $0.01$ to $0.30$ reduces effective convection by about a factor of $2.45$, while chemically passive shear-induced dispersion falls by a factor of $24$. The coupled formulation also yields the exact partition $\Phi_m+\Phi_s+\Phi_a=1$, separating mobile, reversibly retained and irreversibly absorbed solute. The reconstructed field quantifies reaction-modulated transverse non-uniformity, while the results distinguish hydrodynamic effects of confinement and microrotation from kinetic wall effects.

This paper develops a high-order cut energy-based discontinuous Galerkin (CutEDG) method for second-order acoustic wave equations on complex domains and across stationary material interfaces. The method employs unfitted Carte- sian meshes and addresses the small-cut-cell problem by introducing ghost-penalty stabilization directly into the local bilinear forms defining the EDG energy, while leaving the original spatial coupling and numerical-flux structure unchanged. We establish semidiscrete energy stability and a priori error estimates showing that high-order accuracy is retained. The formulation is further extended to interface problems with discontinuous wave speeds. Numerical experiments confirm the pre- dicted convergence rates and demonstrate robustness with respect to arbitrarily small cut cells. In particular, the explicit time-step restriction is governed by the background mesh size rather than the smallest physical cut cell, and the stabilized energy Gram matrices exhibit cut-independent conditioning.

This paper presents approximate computational methods for modeling surface plasmon propagation on a graphene sheet with time- and space-dependent material properties. Starting from Maxwell's equations, the evolution of current density under a spatially and temporally modulated Drude weight is formulated as a partial integro-differential equation (PIDE). Assuming small-amplitude perturbations of the Drude weight, a regular perturbation expansion is developed. The solvability and Fourier invertibility of the resulting order-by-order integral equations are established in appropriate Sobolev spaces. A numerical solver combining a Volterra integral equation formulation with the fast Fourier transform (FFT) is introduced, along with pointwise error estimates and computational complexity analysis. Additionally, an analytical transform-based solution using contour integration and Laplace inversion is derived for traveling-wave Drude weight modulations. Numerical simulations demonstrate diverse wave phenomena, including traveling, standing, and exponentially growing plasmonic current densities.

We give a deterministic algorithm that solves a nonsingular linear system $Ax=b$, where $A\in\mathbb{R}^{n\times n}$ has $m$ nonzero entries and condition number $\kappa$, to any relative residual tolerance $0<\epsilon\le1/2$ using $\widetilde{O(}mn\log(\kappa/\epsilon))$ bit operations for inputs with logarithmically many bits per entry. For sparse, polynomially conditioned systems with $m=\widetilde{O}(n)$, this gives an $\widetilde{O}(n^2)$ algorithm for any inverse-polynomial accuracy, improving on the algorithm of Peng and Vempala, as sharpened by Nie, whose running time in this regime is approximately $O(n^{2.2707})$ with the best current matrix multiplication exponent, and largely closing a gap between the idealized performance of the conjugate gradient method in exact arithmetic and the running time achievable with finite-precision computation that has persisted for over 70 years. The algorithm is surprisingly simple. For integer inputs, we apply Dixon's lifting algorithm to the perturbed system $(A+I/R)x=b$ for a suitable integer $R$. After scaling, the matrix of this system is $RA+I\equiv I\pmod R$, so its modular inverse is trivial, and each lifting step needs only a sparse matrix-vector product with $A$ on $O(\log R)$-bit numbers. Fast rational reconstruction then recovers the exact solution of the perturbed system, which is an $\epsilon$-accurate solution of the original one. Normalization and rounding extend the result to fixed-point and floating-point inputs, with floating-point outputs represented using short integer significands and a common encoded exponent. A computable certificate removes the need for prior knowledge of $\kappa$.

Ship structural vibrations contribute to noise, fatigue, and equipment damage, while dynamic-compliance topology optimization can produce pathological designs near resonance. This study extends neural-reparameterized topology optimization using a convolutional Kolmogorov-Arnold network (KATO) to forced-vibration design with active input power (AIP) as the objective. Applications include a 100 Hz engine-supporting deck panel and an 18 Hz thruster foundation frame. Helmholtz PDE filtering and Heaviside projection control feature sizes and manufacturing tolerance. Across both deck families, all eight optimized layouts reduce AIP relative to size-optimized references and, after finite-depth extrusion, also achieve lower static compliance. For unrestricted, manufacturing-aware, and stress-aware frame variants, KATO matches GCMMA in AIP within 0.5 dB while yielding 22-36x lower static compliance after matched-volume binary re-analysis. In a near-resonant 300 Hz case, both methods reduce initial AIP by more than 32 dB; KATO maintains a connected design, achieves 59x lower binary static compliance, and reduces maximum AIP over 1-500 Hz by 2.7 dB. KATO runs 6.4-10.4x faster than GCMMA for the implemented stress-aware formulations. The results demonstrate neural AIP-driven topology optimization as an efficient approach for designing connected, feature-size-controlled ship structures with improved forced-vibration performance.

Rounding errors in finite element computations can lead to a complete loss of accuracy, stalled convergence, and incorrect results. Moreover, the effects of rounding errors accumulated within automatically generated and compiled kernels are difficult to analyze a priori. We present the first software framework for automated rounding error estimation within finite element kernels. The proposed methodology is based on an a posteriori technique called Running Error Analysis (REA), where a forward error estimate is automatically computed concurrently with the value. An open-source implementation is provided for the FEniCS Form Compiler (FFCx), based on a C++ backend for generating type-generic templated kernels over a custom arithmetic type that tracks both the value and its error estimate. We demonstrate REA on two examples. First, we use it to detect catastrophic cancellation in the assembly of a Neo-Hooke hyperelastic model in the small deformation regime. A series expansion circumvents the cancellation problem and the computed error estimates show this. Second, we study the assembly of the Laplace operator on a near-degenerate mesh. We demonstrate that our REA implementation typically incurs only 2-4x performance overhead. Applications of this work include robust reduced-precision computations in embedded systems, numerical debugging of new, possibly ill-conditioned or unstable PDE formulations, and guiding the design of mixed-precision kernels.

Computing spectra of differential operators on surfaces is crucial in many applications, ranging from medical image analysis to nanoscale quantum structures. However, discretization of the infinite-dimensional operator and the underlying surface can lead to inaccurate and misleading results. This paper introduces the SurfSpec suite of algorithms to compute eigenvalues and spectral measures of operators on surfaces with error control and convergence guarantees. We combine methods based on complex contour integration and residuals to compute discrete spectra and high-frequency eigenmodes, including of nonlinear eigenvalue problems. We then compute convolutions of spectral measures with rational kernels to deal with continuous spectra. In both cases, the key tool is the operator resolvent applied to functions, which is computed using high-order methods that solve surface PDEs. We illustrate our algorithms on a variety of differential operators and surfaces, including a triceratops.

Molecular mixing models for the conditional scalar diffusion (CSD) term in the filtered mass density function (FMDF) transport equation used in Large Eddy Simulation (LES-FMDF) approach are evaluated against direct numerical simulation (DNS) of spatially evolving jets. A database of DNS of turbulent mixing slot jets of $C_{7}H_{16} / O_2$ and $C_{7}H_{16} / Air$ is developed. The formulation includes the compressible form of the governing equations, a generalized multicomponent diffusion model, a cubic real gas equation of state and realistic property models. Simulations are conducted over a wide range of initial pressures ($1 atm < P_0 < 100 atm$) and jet width based Reynolds numbers of 850 and 1300. FMDF of mixture fraction at various filter widths are obtained from the simulation at several spatial locations within the flow. The CSD term in the exact transport equation of FMDF is calculated from the DNS. An a priori analysis of Interaction by Exchange of Mean (IEM), Modified Curl (MC) and Mapping Closure (MAPPING) mixing models for CSD is conducted. At high pressures, significant deviations are observed between the CSD measured from the DNS and those predicted by the models. The source of the deviations is determined to be the mixing frequency used in the models. The significance of the mixing frequency and its variation with change in pressure and diffusion models is studied. The primary objective of this work is to assess the applicability of aforementioned mixing models to flows at large pressures and determine the causes for deviation in prediction.

We study the finite element approximation of the Dirichlet problem for the Bessel $(p,s)$-Laplacian, $\operatorname{div}^{s} \left(|\nabla^{s} u|^{p-2}\nabla^{s} u\right)$, built from the Riesz fractional gradient $\nabla^s$ and posed on the Bessel potential space $X^{s,p}$, obtained by complex interpolation between $L^p$ and $W^{1,p}$. For $p = 2$ this operator reduces to the fractional Laplacian, while for $p \neq 2$ it differs from the fractional $p$-Laplacian arising from real interpolation. Since a direct Galerkin discretization requires reassembling a dense stiffness matrix at every nonlinear iteration, we propose an augmented Lagrangian formulation based on the projected fractional gradient $B_h = \mathbf{\Pi}_h \nabla^s$, in which the only dense matrix is assembled once and the nonlinearity decouples elementwise. We analyze the resulting consistency error, establish a priori convergence rates in $X^{s,p}$, conditional on a discrete inf-sup condition for the augmented Lagrangian scheme, and present numerical experiments that illustrate the convergence of the method and its ability to deal with degenerate and strongly nonlinear problems.

The Runge phenomenon in polynomial interpolation is often considered a classical analogue of the double descent phenomenon in machine learning. In this note, we explore overparameterized polynomial interpolation in three popular polynomial bases: Monomial, Chebyshev and Legendre basis with coefficients that are minimal in the $\ell^2$-norm (and, for the monomial basis, also those minimal in the $\ell^1$-norm). We present our results primarily for equidistant and Chebyshev data points, but many results are independent of the exact form of sampling.

In two and three dimensions, we analyze a mixed formulation of the nonlinear Biot equations in terms of the poroelastic stress tensor, the displacement, the pressure, and the rotation tensor. To establish the well-posedness of the problem, we introduce suitable stabilization terms, leading to an augmented mixed formulation in which the associated parameters are chosen appropriately. Owing to the nonlinearity of the problem, we employ a fixed-point strategy to prove the existence and uniqueness of solutions. We propose a finite element scheme for which we establish existence and uniqueness of a discrete solution, together with a priori error estimates. Additionally, we design and analyze an a posteriori error estimator of the residual type which results to be reliable and efficient. Finally, we present a series of numerical experiments in two and three dimensions in order to assess the performance of the proposed method.