Papers

We develop structure-preserving methods for the compressible Euler equations with gravity in vector-invariant form and potential temperature as a prognostic variable within the flux-differencing discontinuous Galerkin spectral element method (DGSEM) framework. By discretizing the nonconservative terms as symmetric and antisymmetric products, we derive two-point numerical fluxes that conserve both the thermodynamic entropy and the total energy. Moreover, we design an entropy-stable numerical flux that is well-balanced for both isothermal and isentropic background states. All properties are shown to carry over to the high-order DGSEM on general curvilinear meshes. Several numerical examples confirm the theoretical findings and show the robustness and accuracy of the scheme for use in modern dynamical cores for atmospheric flows.

We construct computable error estimators, which double as loss functions for neural networks, for a class of parametric nonlinear partial differential equations with a monotonicity property, and prove that they are globally reliable and efficient. The value of such a loss function is bounded above and below by the squared error in the natural trial norm, for every trial function, not merely for those near the exact solution; this global property rests on monotonicity. The construction rests on splitting the nonlinear operator into a linear part and a strongly monotone closure, and on measuring the linear residual in a discrete dual norm. Since the closure contributes a dual norm that admits no closed form when the trial norm is stronger than an $L_2$ norm, the estimator is built around a computable surrogate for it, required only to satisfy a pairing bound and a Lipschitz bound. The main theorem then yields two-sided bounds with explicit constants and covers two instances. A first-order system least-squares estimator on conforming trial spaces is the first instance studied: its two-sided bound holds on the whole trial space and therefore applies to arbitrary approximations, including those not in a discrete finite element space. The second instance is a discontinuous Petrov-Galerkin estimator on trial spaces of finite element functions for each parameter value, built with broken test spaces, for which no conformity is required and the dual norm is computed by independent element-local problems. All assumptions are verified for a model class of nonlinear fluxes, with the constants tracked explicitly in terms of the parameter range. Being computable for an arbitrary input, both estimators serve as variationally correct loss functions for neural network approximations of parameter-to-solution maps.

In this paper, we introduce a class of function-preserving neural network operators designed to preserve a prescribed strictly positive function. We examine the approximation behaviour of the proposed operators on the space of continuous functions and show that they are positive, linear, well-defined, and bounded. Specifically, we demonstrate uniform convergence on compact intervals and get quantitative estimates of the approximation error in terms of the modulus of continuity. To support the theoretical analysis, we present several numerical experiments focusing on the preservation of exponential and logarithmic functions. These examples show how well the proposed operators perform and compare them with existing neural network operators. Also, an application to denoising modulated signals is provided, demonstrating the practical relevance of the proposed approach.

Given incomplete measurements of a single physical field in a coupled system with unknown parameters, can we infer its full physical state and identify the underlying parameters? This problem is challenging because multiple coupled fields must be reconstructed simultaneously from limited observations of only one, while the system parameters are unknown. In this work, we propose a machine learning framework for full-field reconstruction and parameter identification of unknown physical systems from sparse observations of a single physical field. Specifically, the cross-attention encoder propagates sparse sensor observations onto a regular grid to construct a sensor-conditioned latent representation, while a Fourier neural operator (FNO) decoder captures global spatial dependencies to reconstruct all coupled physical fields. The network parameters and unknown physical parameters are jointly optimized by minimizing observation losses, governing equation residuals, and boundary/initial condition constraints. The proposed approach is validated on two- and three-dimensional lid-driven cavity flows, a two-dimensional cylinder wake, and a two-dimensional non-ideal magnetohydrodynamics problem, demonstrating the recovery performance of unobserved fields and physical parameters from incomplete observations.

We present a structure-aware construction of quantum BPX preconditioners for Q1 finite element discretizations of constant-coefficient reaction–diffusion equations on the unit cube, including the Poisson equation. Building on existing quantum finite element and finite difference methods, we express the preconditioned operator through explicit one-dimensional maps, tensor products, and norm-preserving interlevel transfers. Symmetric preconditioning gives multilevel coefficients, and its favorable conditioning alone does not ensure efficient recovery of a state in the original finite element nodal basis. We instead use the local $L^2$-orthonormal basis underlying the construction as the output representation of the same finite element function. Its coefficient norm equals the continuous $L^2$ norm. Assuming direct access to the preconditioned right-hand side, we prepare this orthonormal coefficient state with a represented function error at most $\varepsilon$. The expected query complexity is \[ \mathcal{O}\!\Big(dL^3\frac{\|f\|_{L^2}}{\|u\|_{L^2}} \log\!\Big(2+\frac{L\|f\|_{L^2}}{\varepsilon}\Big)\Big), \] where $u$ is the exact solution, $f$ is the right-hand side, $d$ is the spatial dimension, and $L=\log_2(1/h)$ for mesh width $h$. The bound holds for $0<\varepsilon\le\|u\|_{L^2}/2$ and depends only polylogarithmically on mesh refinement for fixed continuous data. We also analyze linear-functional estimation and illustrate the construction and continuous error bounds with a one-dimensional UnitaryLab experiment.

The Guided Local Hamiltonian Problem (GLHP) is the problem of ground state energy estimation and preparation given a guiding state having good overlap with the true ground state. The problem is known to be BQP-complete for small error but is classically tractable for larger error. In this work, we identify an orthogonal regime in which the problem is classically tractable, namely when the true ground state $\psi$ has polynomially bounded $\ell_1$-norm. In particular, we show that the Truncated Power Method (TPM), a classical algorithm introduced by Kirby et al. in the quantum context, has a runtime that is polynomially comparable to quantum algorithms for GLHP in all relevant parameters, except for an extra factor of $\|\psi\|_1^2$. In other words, if $\|\psi\|_1^2$ is polynomially bounded, then TPM dequantizes quantum algorithms for GLHP.

We develop and analyze several families of fully discrete, high-order methods in the finite difference time domain (FDTD) framework for a nonlinear dispersive electromagnetic model for ultrashort pulse propagation in a nonlinear optical medium. The model consists of the time domain Maxwell's equations coupled to two ordinary differential equations; the nonlinear Kerr-Debye equation and the linear dispersive Lorentz equation, in one space dimension. The proposed methods combine staggered spatial discretizations of $2M$th-order, for integer $M>0$, with second-order explicit leapfrog or implicit trapezoidal time integrators. To handle the stiff relaxation in the nonlinear Kerr-Debye equation, we adopt a modified exponential integrator, introduced in \cite{PengJCP2020}, which ensures asymptotic-preserving (AP) behavior in the limit of vanishing relaxation time. We prove that the resulting schemes satisfy a discrete analogue of the continuous energy identity and preserve the positivity (PP) of the nonlinear susceptibility. For the schemes based on the leapfrog time integrator, we establish energy stability under a CFL condition, while the schemes based on the trapezoidal time integrator are unconditionally energy stable. High-order temporal accuracy in the trapezoidal time integrator based schemes is further achieved via Richardson extrapolation yielding $2M$th-order accurate methods in both space and time. Numerical experiments demonstrate the theoretical convergence rates and the ability of the proposed schemes to capture nonlinear wave phenomena such as soliton propagation, self-steepening, and odd harmonic generation. This work extends previous AP-PP and energy stable approaches in a discontinuous Galerkin framework to the FDTD setting, offering a computationally efficient, flexible and robust tool for simulating electromagnetic wave propagation in nonlinear dispersive optical media.

We develop and analyse a conforming piecewise affine finite element method for the Grad–Mercier equation, a nonlocal elliptic model for magnetic equilibria whose source depends on the global rank of a flux function. The principal difficulty is that the unique continuum equilibrium has a maximum plateau of prescribed positive area. The solution therefore samples the rank law only above the plateau mass, while perturbations that separate equal values may distribute this mass over several numerical ranks and thereby alter the non-strict rank source. We exploit the continuum well-posedness structure to replace the rank law by an equivalent clipped, nonnegative law, and replace Lebesgue measure by a mass-lumped nodal measure. Equal nodal values are assigned a common rank interval, thereby retaining the full lumped mass of each tied level. The resulting weighted rearrangement gives a strictly convex discrete energy and hence a unique finite element equilibrium. On meshes satisfying a discrete comparison principle, the discrete equilibrium is nonnegative and bounded above by a computable barrier. A Poisson-preconditioned rearrangement iteration converges on each fixed mesh and admits computable a posteriori energy and energy-norm bounds. As polygonal finite element domains converge to the physical cross-section, the discrete equilibria converge strongly in $H^1_0(\Omega)$ and uniformly on $\overline\Omega$ to the unique Grad–Mercier equilibrium. Reconstructed equilibrium sources converge weak-* in $L^{\infty}(\Omega)$ to the original nonlocal source.

A geometrically-exact planar beam finite-element framework is developed for the analysis of beams with hyperelastic constitutive models. A conventional three-field beam model is first formulated with full geometrical nonlinearity and a linear elastic material law. The kinematics are then enriched by an additional global field that permits deformation through the in-plane thickness. The out-of-plane direction is incorporated separately through plane-strain and plane-stress constitutive reductions, allowing the cross-sectional deformation assumptions and constitutive behaviour to be formulated independently. Numerical results for the present benchmarks show that the plane-stress condition gives essentially coincident centreline deformations for the Poisson's ratio values considered. In contrast, the plane-strain condition, which is widely employed in the literature, makes the formulation sensitive to Poisson's ratio coupling, with the beam becoming significantly stiffer at higher Poisson's ratios. Comparisons between Saint Venant-Kirchhoff and Neo-Hookean hyperelastic models further demonstrate close agreement in the global deformation and in-plane thickness response, with small constitutive differences becoming more apparent in the local out-of-plane stretch. The proposed approach offers a simple and consistent framework to integrate hyperelastic material models into planar geometrically-exact beam models.

Surrogate modelling is important to enhance the computational efficiency of uncertainty quantification for physical processes involving random parameters. The advantage of Fourier neural operator (FNO) compared to conventional surrogates is its cross-resolution prediction capability. In the current study, a cross-resolution surrogate based on FNO is built for two-phase transient Darcy flows in porous media with random parameters. A multiresolution Sobolev learning method is proposed for FNO to enhance the cross-resolution prediction accuracy by incorporating first and second-order gradients in the loss function using training samples of different resolutions. The gradients are approximated by resolution-related finite difference schemes compatible with the grid for numerical simulation. Experimental results demonstrate the effectiveness of the new method in training FNO for building cross-resolution surrogates.

We analyse a divergence form stochastic differential equation featuring a Leray-Lions flux function that satisfies standard Carathéodory, growth, monotonicity, and coercivity conditions. We consider two types of noise: multiplicative Lipschitz continuous noise and a transport noise. Our approach utilises the gradient discretisation method, a generic framework encompassing various numerical schemes such as conforming and nonconforming finite element methods, hybrid mimetic mixed methods, and hybrid high order methods. For both types of noise, we propose a stable generic scheme, prove that scheme solutions converge to a strong solution in the probabilistic sense. For the multiplicative Lipschitz continuous noise, we establish the rate of convergence.

We prove a geometric hierarchy of Fourier layer energies for the weakly nonlinear cubic Schrodinger equation on $\mathbb T_\pi\times\mathbb T_1$, retaining its nontrivial resonant dynamics. For initial Fourier support in a rectangle, form nested rectangles with half-widths $(2j+1)$ times the initial ones; layers are the differences between consecutive rectangles. On $0\le t\le T\delta^{-1}$, the energy in layer $j$, measured by its squared $H^1$ norm, is at most $C_N\delta^{2j}$ for $1\le j\le N$; the energy outside the $N$th rectangle is at most $C_N\delta^{2N+2}$. These bounds hold at every fixed depth $N$ under a uniform $H^s$ bound, $s>2$, with constants depending on $N$. The Sobolev bound holds for every fixed $T$ in the defocusing case. We also obtain layer estimates for support bounded in one coordinate and for geometrically decaying initial Fourier tails. The proof uses a filtered Poincare-Dulac normal form: input layers select which nonresonant interactions are removed, while the same layer weights control the interactions left in the equation. Finite coordinate changes and stability estimates transfer this hierarchy to the exact solution.

We develop a family of compact, dispersion-corrected marker-and-cell (MAC)-based finite-difference discretizations for the elastic Helmholtz equation. While dispersion correction is well established for scalar Helmholtz problems, fewer approaches exist for the elastic problems, and most use non-compact stencils. To this end, we introduce a generalized MAC (GMAC) framework that enables compact discretization of the elastic Helmholtz equation. We prove that even though GMAC is non-symmetric, this discretization produces no numerical dissipation and no shear-mode separation. For certain weights, GMAC has lower polarization error than standard MAC. We also introduce an elastic dispersion correction (EDC) process, extending real-shifted-wavenumber dispersion correction methods from acoustic to elastic Helmholtz problems. We prove and demonstrate numerically that GMAC is second-order accurate. Furthermore, we show that EDC substantially reduces dispersion for standard MAC and can provide additional correction in some instances of GMAC. Dispersion analysis shows that the resulting method outperforms existing compact finite-difference discretizations in dispersion accuracy and is competitive with some non-compact schemes. We give numerical evidence for improved phase alignment in mildly heterogeneous media.

This paper develops a mesh-free numerical method for the one-phase Mullins-Sekerka flow, in which the interface bounding a set moves as the Wasserstein gradient flow of its perimeter. One difficulty in computing this flow is that two pieces of the interface with opposite normals form a hidden boundary when components merge. Our approach represents the interface by an oriented point cloud, which we treat as an oriented varifold. We introduce a quantity that measures how much nearby normals cancel each other and use it to detect the hidden boundary. Based on this quantity, we construct a perimeter estimator, prove that its expected error decays at an algebraic rate in the number of points, also in the presence of a hidden boundary, and use it in a minimizing movement scheme for the flow in two dimensions. Numerical experiments show that the scheme computes mergers of components through the hidden boundary, while the perimeter term in the scheme decreases at every step and the area is nearly conserved. As hidden boundaries also appear in other perimeter-driven interface problems, we expect our oriented representation to be applicable to them as well.

A second-order exponential low-regularity integrator is proposed for the incompressible Navier–Stokes equations. The method requires only linear solves and satisfies a discrete energy-decay property. Second-order temporal convergence is established under a uniform $H^4$ solution bound, with an error constant that remains bounded as the viscosity tends to zero. For the Fourier–Galerkin discretization, the optimal spatial convergence rate is established for fixed viscosity. The method can also be extended to conforming finite element methods while retaining the discrete energy-decay structure. Numerical results confirm the convergence orders and robustness with respect to the viscosity.

We develop a high-order discontinuous Galerkin formulation of the general synthetic iterative scheme (GSIS) for steady rarefied gas mixture flows. The resulting GSIS-DG discretizes the mesoscopic kinetic equations for monatomic gas mixtures and the macroscopic synthetic equations with the same polynomial degree and solves them alternately. The kinetic solution supplies nonequilibrium corrections to the stress and heat flux, while the synthetic equations accelerate the convergence toward the steady state. Asymptotic analysis shows that GSIS-DG recovers the continuum limit without requiring the spatial mesh to resolve the molecular mean free path. Numerical results over a wide range of Knudsen numbers demonstrate high-order spatial accuracy and rapid steady-state convergence, with more than an order-of-magnitude reduction in computational cost relative to conventional kinetic iteration in the near-continuum regime. Comparisons with second-order finite-volume GSIS further demonstrate that GSIS-DG achieves comparable or better accuracy on coarser physical meshes with fewer total spatial degrees of freedom.

The adjoint method is a standard approach for computing gradients in PDE-constrained optimization and central to parameter estimation, model calibration, and optimal control. This tutorial presents a discrete, implementation-oriented introduction to the method with an emphasis on conceptual understanding. Starting from an inverse problem governed by a partial differential equation, it develops the corresponding optimization problem and derives the adjoint gradient from the perspective of iterative descent methods. The roles of the state, the parameters, the objective, and the constraint are made explicit, and the adjoint equation is obtained directly from the gradient computation rather than postulated. Particular attention is given to explicit derivations, implementation details, and the connection between the mathematical development and the resulting computational algorithm. The tutorial also examines the influence of discretization on the adjoint construction, explains the origin of the backward-in-time adjoint recursion, and presents practical procedures for gradient verification. A vibrating-string inverse problem serves as a running example, illustrating the complete workflow from model formulation to adjoint-based optimization. The presentation is intended for advanced undergraduate students, graduate students, and practitioners seeking a step-by-step understanding of adjoint gradient computation.

Gaussian elimination with partial pivoting is the standard algorithm for solving a dense $n\times n$ system of linear equations $Ax=b$. The “growth factor” $\rho\,$ for this process may be as large as $2^{n-1}$, and when $\rho$ is large, the algorithm is unstable. Nevertheless, long experience has shown that Gaussian elimination is resoundingly stable in practice. Experiments indicate that among random matrices with independent normally distributed entries, growth factors $\rho \gg n^{1/2}$ are in a precise sense exponentially rare. Here we take a step toward confirming this observation by proving it not for the full growth factor, which depends on all the entries of the upper-triangular factor $U$, but for the corner entry $u_{nn}^{}$.

This paper refines randomized block coordinate descent by introducing a coordinate-level batch-sampling distribution characterization, yielding a fine-grained randomized batch-sampling coordinate descent framework, which applies to any block coordinate descent method under an arbitrary static stochastic sampling rule for solving large-scale linear least-squares problems. Based on this framework, we derive new linear convergence rate bounds that are invariant to the column scaling of data matrices. Theoretical analysis and numerical experiments validate that the new bounds are significantly sharper than existing results and better aligned with the practical convergence rate. Furthermore, we proposed a learning-based sampling strategy by parameterizing the batch-sampling distribution and optimizing its parameters through gradient-based training. Preliminary numerical results on computed tomography (CT) reconstruction problems demonstrate that the resulting method converges significantly faster than variants using other fixed sampling rules, and outperforms tested adaptive coordinate descent methods and full gradient descent in computational efficiency, confirming that learned sampling distributions can effectively capture problem-specific structure to accelerate convergence within the framework.

Numerical Linear Algebra (NLA) has consistently played a vital role in advancing science by providing tools to solve fundamental problems encountered in scientific and engineering applications. Over the decades, it has continually evolved to meet the demands driven by successive waves of scientific discovery. For instance, during the 1950s and 1960s, substantial efforts were devoted to developing methods for solving eigenvalue problems that emerged from the rapidly growing field of aerodynamics. This led to the discovery of the LR and QR algorithms. Later the attention turned to the solution of sparse linear systems that were common in applications like computational aerodynamics. Today we are experiencing yet another wave of major scientific advancement and NLA is once more at the heart of its development. This Machine Learning (ML) wave is proving to be utterly disruptive in science and engineering. Many tools in ML particularly Large Language Models (LLMs) are grounded in matrix and tensor methods. As we are approaching Artificial General Intelligence (AGI), it is clear that matrix methods will be called to play an even more significant role. For the numerical linear practitioner the speed of the current change makes it particularly challenging to adapt. This is a survey article that centers on machine learning techniques, with a particular focus on large language models. It has two main objectives. The first is to clarify the core concepts behind Large Language Models in a manner accessible to specialists in numerical methods. The second is to examine the key Numerical Linear Algebra concepts employed by LLM techniques, while also highlighting several significant recent contributions of NLA to the field.

An $\mathcal H^2$ matrix represents well-separated interactions using nested low-rank bases while retaining nearby interactions explicitly. We consider the problem of constructing an explicit, strongly admissible $\mathcal H^2$ approximation using only products with the matrix and its adjoint, and develop two randomized recovery algorithms on an adaptive geometric hierarchy. Carried recovery begins with a single pair of dense Gaussian test matrices, without hierarchy-specific probing, and updates the resulting forward–adjoint sketches across levels. All operator samples can therefore be acquired in one batch. Resketched recovery applies the same recovery procedure using fresh Gaussian test matrices at each level. For this variant, we prove a global quasi-optimal bound on the expected Frobenius error relative to the best approximation with the prescribed hierarchy, rank bounds, and retained near-neighbor pattern. Carried recovery requires only $O(k+p)$ applications of the matrix and its adjoint, independent of tree depth, whereas resketched recovery requires $O((k+p)\log(N/k))$ on trees of logarithmic depth. Here $N$ is the matrix dimension, $k$ the maximum prescribed box rank, and $p$ the oversampling parameter. These bounds assume bounded near-neighbor counts and $O(k)$ local dimensions. Experiments on Laplace, Dirichlet-to-Neumann, and Helmholtz operators in two and three dimensions with up to $1.28$ million unknowns show that, with sufficient oversampling, carried recovery achieves nearly refinement-independent error and accuracy comparable to resketched recovery while using many fewer operator applications. Its compressed application cost scales linearly for a fixed rank.

We establish cellwise pressure monotonicity for a mixed finite element discretization of the rotation-based Biot system. The spatial discretization is based on compatible discrete de Rham complexes, using standard lowest-order elements for rotation, displacement, Darcy flux, and pressure. For the temporal discretization, we employ the backward Euler method. Exactness of the discrete complexes reduces the mechanics contribution to the pressure Schur complement to a scaled pressure mass matrix, yielding a positive diagonal mechanics capacity. For the Darcy part, positive circumcentric facet weights yield a weighted cell graph Laplacian after flux elimination, while the mass-lumped product remains uniformly equivalent to the standard Darcy product. The diagonal mechanics capacity and the weighted cell graph Laplacian make the condensed pressure matrix Stieltjes. Consequently, the pressure propagator is entrywise nonnegative, which yields a cellwise discrete maximum principle. Numerical experiments confirm the predicted convergence rates and show that the proposed discretization suppresses spurious cellwise pressure oscillations.

This paper presents a locally conservative Enriched Linearized Neural Network (ENN) method for the Darcy flow model. A linearized shallow ReLU$^k$ network, whose hidden-layer parameters are fixed on a quasi-uniform set of the sphere, carries the approximation, and piecewise constant functions on an auxiliary subdivision of the domain are added to ensure local mass conservation. The subdivision may consist of general, possibly curved, elements. Since neural network functions do not satisfy the inverse inequality, the stability and error analysis of the standard discontinuous Galerkin interior penalty formulation do not apply. The missing inverse inequality is remedied by an edge identity combined with a Galerkin least-squares (GLS) term, together with a single neural approximant that is optimal in the $L^2$, $H^1$ and $H^2$ norms simultaneously. We use the nonsymmetric interior penalty (NIPG) form, which is coercive for every positive penalty parameter. With $n$, the width of the shallow neural network, used in ENN, we prove an error estimate in the EG norm that is of the optimal $H^1$ order, i.e., $O(n^{-(r-1)/d})$ for solutions in $H^r$, $r\ge2$, and $O(n^{-(s-1)/d})$ for solutions in $H^s$, $\frac32<s\le2$, for any shape-regular subdivision, when the GLS parameter satisfies $\min\{h,n^{-1/d}\}^2\lesssim\tau\lesssim n^{-2/d}$ and the edge terms are scaled with $\min\{h_e,n^{-1/d}\}$, i.e., by the resolution of the network rather than of the mesh. A numerical Darcy flux is constructed, which is locally conservative, regardless of the size or the number of the elements. Sample numerical results demonstrate the correctness of the theory. In particular, a fixed $4\times4$ subdivision gives the same accuracy as a subdivision refined with the network, and the mass loss stays at round-off unlike a non-conservative neural network method.

Despite the growing popularity of physics-informed neural networks (PINNs), their applicability in the long-time integration of partial differential equations (PDEs) remains constrained. We argue that this problem stems from the lack of consideration of temporal causality in the original PINN formulation, resulting in a bias towards satisfying governing equations at later times before learning the initial condition and hence leading to erroneous solutions. To this end, we propose a novel method that seamlessly integrates temporal causality into the training process. Drawing inspiration from classical numerical methods where the temporal causality is reflected, we divide the time domain into nonoverlapping subintervals, assign a unique neural network to each subinterval, and construct a loss function founded on the integral form of PDEs within these subintervals. The proposed networks undergo sequential training, beginning with the initial time step. Our method demonstrates significant improvement in accuracy for long-time simulations of various PDE problems where the original PINN method fails while it requires less computational cost and memory compared to the PINN method. A parallelization algorithm is provided to further enhance the computational efficiency, showing a significant speedup for solving time-dependent PDEs.

We develop arbitrarily high-order structure-preserving parametric finite element methods for geometric PDEs including surface diffusion and volume-preserving mean curvature flow based on local pullback flow maps. The central idea is to reformulate the geometric PDEs on an arbitrary intermediate hypersurface and represent the subsequent evolution through a local pullback flow map. Using admissible local pullbacks, we derive two formulations from the corresponding pullback identities for the Laplace-Beltrami operator, referred to as the conformal and equidistribution formulations. We discretize both formulations using arbitrary-degree isoparametric finite elements and Runge-Kutta collocation methods with positive weights. For the conformal formulation, weighted averages of the normal-Jacobian vector ensure exact preservation of the enclosed area or volume, while algebraic stability of the Runge-Kutta method additionally guarantees dissipation of the perimeter or surface area. For the equidistribution formulation, a path-averaged normal-Jacobian vector and a geometric discrete gradient guarantee both exact preservation of the enclosed area or volume and dissipation of the perimeter or surface area, without requiring algebraic stability. Numerical experiments confirm the expected high-order accuracy and structure-preserving properties.

In this paper, we first study the singular value distribution of the unscaled matrices arising from Jacobi weighted histopolation. For the Jacobi weight with parameters $(\alpha,\beta)$, and using the Jacobi basis with parameters $(\alpha+1,\beta+1)$, we previously proved that, for $\alpha,\beta>1/2$ and quasi-uniform meshes, the growth of these matrices in the sense of singular value distribution is slower than any positive power of the logarithm of the matrix size. However, the numerical behavior of the unscaled matrices suggested a nontrivial singular value symbol, which we identify explicitly in this paper. Starting from the standard asymptotic formula for Jacobi polynomials, we derive an asymptotic representation of the weighted cell averages defining the entries of the histopolation matrices. This representation allows us to analyze the associated Gram sequence by means of GLT theory. We identify its GLT symbol and, consequently, we obtain the singular value distribution of the histopolation matrices for any $\alpha,\beta>-1$. We also characterize the set $\mathcal Z$ where the resulting symbol vanishes and compute its measure, showing that this set depends only on the mesh, while the Jacobi parameters affect the values of the symbol on its support. When this measure is positive, the singular value distribution implies the discrete asymptotic ill-posedness of the problem, while a comparison with the Toeplitz setting suggests exponential growth of the spectral condition number. When the measure of $\mathcal Z$ is zero, the problem is well-posed and the observed conditioning is algebraic in the matrix size. Finally, numerical experiments supporting the theoretical findings are reported and discussed, while conclusions are given in the final section of the current work together with a few open problems.

We give checkable positive-definiteness criteria for the limiting neural tangent kernel (NTK) and high-probability finite-width gradient-descent guarantees for scalar- or vector-valued physics-informed neural networks (PINNs) with linear constraints. The constraints may be strong-form differential rows of any fixed finite order, including coupled systems, with any linear initial or boundary conditions; weak-form residual and boundary functionals; or finite-measure nonlocal observations such as integral, nonlocal-diffusion, and Dirac-data rows, all in any dimension. For each class, positive definiteness of the limiting NTK is equivalent to a rank condition on a finite coefficient, functional, or moment matrix of the fixed design: a certificate computed before training that detects structural zero modes. The model hypotheses are those of an ordinary two-layer PINN, a smooth nonpolynomial activation with bounded symmetric initialization, and are met by standard choices such as $\tanh$ with uniform initialization. Given a certificate, explicit width and step-size conditions ensure, with high probability, that the empirical NTK retains at least half the limiting gap and that full-batch gradient descent decreases the training loss geometrically at every iteration. Controlled experiments check the certificates and the finite-width mechanism.

In this paper we devise and analyze Lie Trotter splitting schemes for non-autonomous evolution equations, following two approaches: 1) A fully non-autonomous scheme, motivated by the potential use of explicit solutions of some subproblems. 2) A splitting scheme relying on the concept of Magnus integrators, where we freeze the operators at the beginning of each time interval as this suffices for first order convergence of the Lie-Trotter integrator. We derive first order convergence in a general setting and apply the analysis to the Bloch-Torrey equation of magnetic resonance imaging, where a decomposition into three subproblems is well-motivated by the physical role of the underlying operators. Numerical tests, complementing the time stepping scheme by a finite element space discretization, illustrate the theoretical findings.

We present a method for the data-driven inference of Hamiltonian and Lie-Poisson dynamical systems from vector-field samples, in which the Hamiltonian and the underlying Poisson tensor are estimated jointly rather than assuming the structure is known. The Hamiltonian is represented as a linear combination of random feature maps with fixed random internal weights, so its estimation is a linear regression problem for any given Poisson structure. We treat three cases of increasing generality: a known Poisson tensor, an unknown constant invertible (symplectic) structure matrix, and Lie-Poisson and affine Poisson structures, with structure matrix linear in the state. In the constant case, the Hamiltonian and structure matrix are recovered jointly from an eigenvalue problem that avoids the trivial zero solution due to the scaling symmetry of Hamiltonian systems. For Lie-Poisson and affine Poisson systems we alternate between an initial Hamiltonian estimate from a generalized eigenvalue problem for an approximate first integral, and a constrained nonlinear least-squares fit of the structure constants that enforces the Jacobi identity and, for affine structures, the cocycle condition, via continuation in a penalty parameter. Orthogonal matching pursuit selects a compact, well-conditioned subset of random features. Because the learned Hamiltonian is a sum of single-feature terms, the model admits an explicit Poisson/symplectic splitting integrator, which we use to validate the learned systems via long-time Poincare sections. We demonstrate the approach on the four- and five-dimensional Lorenz-86 model and a six-dimensional Kirchhoff rigid-body system, recovering Hamiltonians, structure constants and Casimir invariants to high accuracy, and showing the inferred Poincare sections converge to the correct invariant-tori topology as the number of random features increases.

We study the involution-preserving properties of the Active Flux method on unstructured triangular meshes. Two model problems are considered: the two-dimensional linear acoustic system, which satisfies a curl involution on the velocity field, and the two-dimensional Maxwell equations with an in-plane magnetic field, which satisfy a divergence involution. Both systems can be written in conservative form. The question we address is how a single discretization can satisfy the conservative form and the involution at the same time. Within a discrete Active Flux framework on simplices, we derive explicit formulas for a discrete curl operator that is exactly preserved by the scheme for acoustics, and for a discrete divergence operator preserved for the Maxwell system. The analysis applies to general simplex meshes and to any polynomial degree, including extensions of the method beyond third order. Numerical experiments with the third-order scheme confirm the theoretical results. Highlights: 1-Active Flux on unstructured triangular meshes preserves discrete involutions, 2-Explicit discrete curl operator preserved by the scheme for linear acoustics, 3- Explicit discrete divergence operator preserved for 2D Maxwell equations, 4- Results hold on general simplex meshes,% and for any polynomial degree 5- Third-order numerical experiments confirm the theoretical results. \end{itemize}