Papers

We prove positive definiteness and quantitative condition-number bounds for continuous Sobolev Gram matrices of deterministic shallow networks with Gaussian and tanh-difference activations. The prescribed feature families include the optimal Sobolev approximation constructions of Part I. On every fixed nonempty bounded Lipschitz domain in $\R^d$, $d\ge2$, and for every fixed integer $m\ge0$, the raw and column-normalized Gram matrices in the full $H^m$ inner product have spectral condition numbers bounded by $\exp(CN^{2/d}/\log(N+2))$ for sufficiently large $N$, where $N$ is the feature count and $C>0$ is independent of $N$. The upper bounds hold for midpoint and Gauss-Legendre offsets. For midpoint offsets and a suitable angular sequence, matching lower bounds on fixed ellipsoids establish sharpness of this exponent over the stated domain class for raw Gaussian matrices and both raw and normalized tanh-difference matrices. The proof combines homogeneous differential operators that isolate direction groups with one-dimensional coefficient-recovery estimates. Two-dimensional $L^2$ experiments on a square illustrate condition-number growth and normalization effects.

We construct shallow neural networks with fixed Gaussian and tanh-difference activations on bounded Lipschitz domains $\Omega\subset\mathbb{R}^d$, $d\ge2$. For each $s>0$, explicit deterministic rules prescribe all hidden parameters independently of the target $u\in H^s(\Omega)$. At frequency resolution $M$, the $N_M\asymp M^d$ prescribed features span a linear space $V_M$. For every sufficiently large integer $N$, we choose $M\asymp N^{1/d}$ with $N_M\le N$ and construct a bounded linear operator $A_N:H^s(\Omega)\to V_M$ such that \[ \|u-A_Nu\|_{H^m(\Omega)} \le C N^{-(s-m)/d}\|u\|_{H^s(\Omega)} \] for every $u\in H^s(\Omega)$ and every integer $0\le m<s$, with $C$ independent of $u$, $N$, and $m$. Sobolev width lower bounds show that this rate is optimal among linear spaces of dimension at most $N$ chosen independently of the target. We use explicit spherical cubature nodes for the directions and Gauss-Legendre nodes or equispaced midpoints for the offsets. We show that the exterior offset margin can shrink to zero as $M\to\infty$ while preserving the approximation rate. A margin beyond $\sup_{x\in\Omega}|x|$ of order $M^{-1}\log M$ for the Gaussian and $M^{-1}(\log M)^2$ for the tanh difference suffices to truncate the ridge integral and remove exterior midpoint nodes. Numerical experiments in two and three dimensions examine $L^2$ and $H^1$ approximation.

We propose an energy-stable arbitrary Lagrangian-Eulerian (ALE) finite element formulation for incompressible Navier-Stokes flows with surface-tension-driven free boundaries. The key ingredient is a multiplier-free compensated tangential velocity formulation, which redistributes boundary nodes via the discrete tangential curvature residual. The relaxation parameter governing the artificial tangential velocity is scaled according to the viscous-capillary balance, yielding the natural choice $\alpha= \gamma_0/ \mu $, where $\gamma_0$ is the surface tension coefficient and $\mu$ is the dynamic viscosity. We prove that the fully discrete scheme is well-posed, energy-dissipative, and temporally consistent with its semi-discrete counterpart. Numerical results confirm the robustness of the proposed method across a range of $\mu$ and $\gamma_0$ values. Furthermore, the method exhibits superior mesh quality preservation compared with the BGN method, particularly in the small-time-step regime.

We derive error bounds from empirical squared residuals for two non-monotone, mesh-free approximations of Hamilton–Jacobi–Bellman equations on unbounded spatial domains: collocation with Wendland kernels under a native space norm constraint, and physics-informed neural networks with smooth activations under a Sobolev-norm constraint. The main step is to recover the pointwise residual control needed for viscosity-solution error estimates from the empirical loss, while quantifying the discretization or sampling error. For the kernel method, we prove uniform convergence on compact sets under classical and native space regularity and derive an a posteriori quantitative error bound for the unique bounded continuous viscosity solution. For the neural-network method, we prove uniform convergence on compact sets in probability and give a quantitative high-probability error bound under classical Sobolev regularity. The bounds are governed by the fill distance for kernels and by the network and sample sizes for neural networks. The proofs combine empirical-to-population estimates, an interpolation inequality that turns an $L^2$ residual into a sup norm bound, and a controlled diffusion estimate for the error caused by truncating the domain. The regularity assumptions used to obtain convergence from approximation properties are additional to viscosity wellposedness. Numerical experiments in one and two space dimensions illustrate the methods and the estimates.

We study explicit lightning-plus-polynomial (LP) approximation of $g(z)z^\alpha(\log z)^m$ on a closed unit-radius sector, where $0<\alpha<1$, $m$ is fixed, and $g$ is analytic nearby. Since the basic sector root-exponential rate is already known, we address three finer questions. First, a one-sided Abel–Poisson formula and a pole-free notched contour give a branchwise error bound that is uniform for grid shifts in compact subsets of $(0,1)$ and valid for every fixed logarithmic power. Second, for $z^\alpha$ in the dense and balanced regimes, an exact vertex-scale rescaling identifies a positive finite leading constant for a specified shifted sequence. Third, a rate-matched band of additional representation poles followed by polynomial compression gives the explicit representation dimension $N+2\sqrt{2\alpha\lambda_\sigma N}+O(1)$ while leaving the first $N$ tapered poles unchanged. A slightly overresolved hybrid sequence also preserves the exact leading constant. Matched-dimension experiments, including reentrant sectors, illustrate the allocation. These exact-asymptotic and dimension statements concern the stated constructions, not sector minimax approximation or total computational complexity.

For $h>0$, degree $m\ge0$, $0\le\beta<\pi$, and $q>0$, let $\Sigma_\beta=\{re^{i\theta}:r\ge0,\ |\theta|\le\beta/2\}$ and let $E_{m,\beta}(q)$ denote the best uniform error for approximating $e^{-hz}$ on $\Sigma_\beta$ by rational functions whose only possible finite pole is $s=-q/h$, with multiplicity at most $m$. For the linear scaling $q=\alpha m$ with a fixed $\alpha>0$, a Möbius map reduces the problem to polynomial approximation on a compact lens. Write $\gamma_\beta=2-\beta/\pi$ and let $\alpha_c(\gamma_\beta)$ denote the saddle-coalescence threshold. For $0<\alpha<\alpha_c(\gamma_\beta)$, let $H_\beta(\alpha)\in(0,1)$ be the exponential of the real saddle action. An exact boundary Faber integral, a global fractional-power contour deformation, and a two-saddle expansion yield the sharp two-sided order $E_{m,\beta}(\alpha m)\asymp m^{-1/2}H_\beta(\alpha)^m$, uniformly on compact subsets of $0<\alpha<\alpha_c(\gamma_\beta)$. On this noncoalescing interval, the exponential factor $H_\beta(\alpha)$ has the unique minimizer $\alpha_\beta^*=1/[\gamma_\beta\sin(\pi/(2\gamma_\beta))]$ and minimum $\tan(\pi/(4\gamma_\beta))$, recovering the classical half-line factor $\tan(\pi/8)=\sqrt{2}-1$ at $\beta=0$. Numerical-range estimates give matrix dimension-independent upper bounds for projected shift-and-invert Arnoldi, while finite-degree pole searches and nonnormal experiments illustrate the distinction between the sector-uniform law and individual matrix computations.

Entropy stable boundary conditions are critical for ensuring that the corresponding numerical scheme satisfies the discrete entropy inequality, thus, mimicking the second law of thermodynamics for the compressible Navier-Stokes equations discretely. Provably entropy-stable adiabatic wall boundary conditions were introduced in Parsani 2014~\cite{parsani2014entropy}. These discrete wall boundary conditions were further generalized in Dalcin 2019~\cite{dalcin2019conservative} to a moving adiabatic solid wall or a wall with a prescribed heat flux for the compressible Navier–Stokes equations discretized by using summation-by-parts (SBP) and simultaneous-approximation-term (SAT) operators. Currently no entropy-stable formulation using SBP-SAT operators for isothermal wall boundary conditions is available in the literature. This paper presents a new isothermal no-slip wall boundary conditions that: 1) enforces the $T|_{y=0}=T^{wall}$ and $\vec{v}|_{y=0}=\vec{v}^{wall}$ conditions with the design order of accuracy while maintaining stability, 2) provides the correct sign of entropy production at the wall, and 3) mimics the entropy balance of the Navier-Stokes equations at the discrete level. The structure-preserving and design-order properties of the proposed methodology are demonstrated and verified on standard benchmark problems for compressible flows.

This work analyzes a framework for learning the solution operator, or equivalently, the Green's function, of a general self-adjoint second-order elliptic operator $\cL$. Subtracting the singular part of the Green's function, either known analytically or supplied by an explicit Laplace parametrix, recasts the learning of $\cL^{-1}$ as learning an $\cL$-harmonic extension operator on the boundary. We formulate the process as operator learning between reproducing kernel Hilbert spaces (RKHSs) to show that the kernel decomposition strategy raises the algebraic input rate by two orders and, under uniform analytic regularity of the induced boundary data together with compatible output spaces, upgrades the output rate from algebraic to exponential. We extend the framework further to semilinear problems via an iterative decomposition in which the only learned object remains a linear operator, ensuring that the improved linear rates translate directly to the nonlinear iteration. Numerical experiments are provided to validate the theoretical development.

In this paper, we present JAXIGA, an open-source library for isogeometric analysis based on JAX. The library separates the construction of the mesh from the computations which are traced for compilation and differentiation. This allows geometry evaluation, assembly and the supported solvers to be combined in a just-in-time compiled program, for which first-order reverse-mode derivatives can be computed. We consider three numerical methods obtained from a single pointwise energy density. The Galerkin and energy-minimization methods use the same discrete potential and have the same stationary points, while collocation enforces the corresponding strong form. The sensitivities of a converged solution are obtained by solving one adjoint problem, and their accuracy depends on the tolerances of both the primal and adjoint problems. We also describe hierarchical extraction for local refinement of the variational methods and the use of smooth spline spaces for fourth-order problems. The numerical results are compared with analytical solutions and independent FEniCSx and JAX-FEM computations. For the reported GPU benchmark with 36,448 unknowns, assembly is twenty-four times faster than on the host CPU. When the conjugate-gradient solver is also run on the GPU, the total solution time is reduced from 31.8 to 1.38 seconds, and the result agrees with the direct solution to $5\times10^{-11}$. These results show that the choice of solver and the reuse of compilation are important for obtaining an overall reduction in computational cost.

Toeplitz matrices are among the simplest structured models of translation invariance, arising naturally in convolution, stationary covariance models, and discretizations of translation-invariant operators. For a pair of real matrices, the Böttcher–Wenzel form (BW form) is a nonnegative quartic measuring the gap in the corresponding commutator inequality. László conjectured that, when both matrices are Toeplitz, this quartic is a sum of squares for every matrix order. We show that the conjecture is false. More precisely, we identify an explicit family of corner-supported Toeplitz pairs, retaining only the outermost diagonals on the two sides, for which the associated BW form is not a sum of squares at all sufficiently large orders. The proof reduces the complete family of Gram representations to a finite complex test: its Gram-invariant part develops a strictly negative limiting direction, while the remaining Gram-dependent contribution vanishes. We also give an explicit sufficient threshold for the matrix order. The negative result is complemented by two positive regimes. If either Toeplitz factor is symmetric or skew-symmetric, then the corresponding BW form is SoS for every order, with the other factor allowed to be an arbitrary real Toeplitz matrix. For unrestricted Toeplitz pairs, we further prove SoS representability for every $2\le N\le50$, with formal Lean~4 verification through order 20. The large-order non-SoS theorem, including the explicit sufficient threshold, has also been formally verified in Lean 4.

Experience with invariant finite-difference schemes for one-dimensional hydro\-dynamic-type equations for a function $u(t,x)$ suggests that nonlinear polynomial schemes on orthogonal uniform meshes with the 9-point stencil cannot simultaneously preserve difference analogues of the energy and momentum conservation laws, with the multipliers given by the standard difference approximations of $u_t$ and $u_x$. For the 9-point stencil this hypothesis is proved for polynomials of arbitrary degree. On larger stencils it is investigated by means of the difference analogue of the direct method: for polynomials of bounded degree, such schemes are found to exist only on sufficiently large stencils, and the second-order equations they approximate are of Monge–Ampère type. As a by-product, a scheme for the homogeneous Monge–Ampère equation is obtained that has five conservation laws, admits most of the symmetries of the differential equation, and is invariant under the discrete dihedral group $D_4$. Numerical tests on exact solutions confirm its second order of approximation and show that the conservation laws are satisfied on these solutions with the same order of accuracy.

We prove the convergence of the isoparametric Barrett–Garcke–Nürnberg (BGN) method for mean curvature flow of closed surfaces. For finite elements of degree $k\ge6$ and time steps $\tau=O(h^2)$, the numerical surfaces converge to the exact surfaces in the manifold distance with order $h^2$. The proof relies on two ingredients. First, we construct approximation flows by an expansion in time, whose interpolations satisfy the BGN method with an improved defect. Second, we derive a refined estimate for the tangential error from an algebraic identity for tangential test functions.

We propose an algorithm for the numerical solution of the time-dependent Schrödinger equation for N >> 1bosons interacting via a two-body interaction potential in the mean-field scaling regime. Unlike existing numerical algorithms for quantum dynamics, the proposed algorithm has a computational complexity that is independent of the number N of particles. It is based on an ansatz for the wave function inspired by the Bogoliubov theory, which is an asymptotic theory for large N. In a first part, we follow Bogoliubov theory in evolving the condensate function by the Hartree equation, but we use the Dirac–Frenkel time-dependent variational principle to derive nonlinear equations of motion for the excitation functions in the Bogoliubov ansatz. By carefully exploiting the bosonic symmetry, the so obtained approximate wave functions are shown to have derivatives bounded independently of N apart from a time-dependent global phase factor which rotates with a frequency proportional to N. In the second part, we numerically approximate the wave function (modulo the phase) by a composite algorithm, which employs an appropriate splitting of the equations of motion and uses a dynamical low-rank approximation of the fundamental two-particle excitation function by linear combinations of products of single-particle functions. The low-rank excitation function is evolved in time by a basis update and Galerkin (BUG) integrator that is robust to small singular values of the coefficient matrix and allows for rank adaptivity. Numerical results are provided to assess the performance of the proposed algorithm.

We prove that the probability that Gaussian elimination with partial pivoting on $n \times n$ random matrices has growth $\rho$ is at least inverse quasi-polynomial: $\Omega(\exp(-c \log^2 (\rho)\log(n)))$ for some constant $c > 0$. This lower bound breaks standard conjectures in the smoothed and average-case analysis of Gaussian elimination. To the best of our knowledge, it is the first non-trivial lower bound on the probability of large growth for Gaussian random matrices.

In this paper we consider matrices whose entries are rational functions of several parameters over a Moore-Penrose field, that is, a field with an involutory automorphism over which every matrix has Moore-Penrose inverse. We prove that over any such field the Penrose conditions can be solved with linear algebra alone, even though they form a polynomial system of degree two. We then bound the degrees of the numerator and of the denominator of the pseudoinverse in terms of the degree of the entries and of the rank, and not of the dimensions of the matrix. Specialization is the thread: the closed form is computed once over the field of rational functions, and the question is at which parameter values it still returns the pseudoinverse of the specialized matrix. We determine the values where it does not, from the matrix alone and before computing the pseudoinverse. They are the zeros of a polynomial built from the maximal minors, and they are the values at which the rank decreases. This turns previous sufficient conditions into an exact characterization. The theory also gives symbolic algorithms for the real, the complex and the parametric case, which we implement in Maple and test on 972 timed replications. We also apply them to a Leontief economic model, where the excluded value is the point at which the economy stops being viable.

Multiscale, parameter-dependent partial differential equations (PDEs) pose severe computational challenges due to strong coefficient heterogeneity and high-dimensional parameter spaces. We develop a geometric interpolation approach within the multiscale generalized finite element method (MS-GFEM) that targets the most expensive component: computing parameter-dependent optimal local approximation spaces. Leveraging the spatial localization of MS-GFEM and assuming local parameter dependence, we decompose the global problem into parametrically low-dimensional local subproblems. The optimal subspaces for each parameter are identified as points on a Grassmann manifold and approximated via Grassmann interpolation on sparse grids, which preserves the geometric structure of these spaces while efficiently handling high-dimensional parameter spaces. The resulting localized model reduction method inherits the nearly exponential spatial convergence of MS-GFEM and the parametric convergence rates of sparse grids. Numerical experiments for elliptic problems confirm the theoretical convergence results.

We address the Subset selection problem for matrices, where the goal is to select a subset $\mathcal{S}$ of $k$ column indices from a \enquote{short-and-fat} matrix $X \in \mathbb{R}^{m \times n}$, such that the sampled submatrix $X_{\mathcal{S}}$ has $\|X_{\mathcal{S}}^\dag X\|_F$ as small as possible. Our approach is centered on volume sampling, which attains the tightest known bound on this objective in expectation. As our primary contribution, we propose a new deterministic algorithm, Forward derandomized volume sampling (FDVS), which provably attains this bound and has asymptotic complexity $O(nkm)$. In contrast, the complexity of all previously known algorithms with this guarantee scales at least quadratically in $n$ when $k \ll n$, making FDVS particularly attractive for very wide matrices. In addition, we systematically structure the landscape of volume sampling methods: we propose a simple $O(nm^2)$ algorithm for exact forward volume sampling, show that a known deterministic method is in fact a derandomization of Reverse iterative volume sampling, and derive a modification of the latter that avoids repeated SVD downdating, along with a fast greedy forward variant. The algorithms are verified and compared in numerical experiments involving optimal experimental design, sensor placement, and DLRA-DEIM.

Advanced electronic packaging structures represent a key technological approach for extending Moore's Law. An electronic packaging structure consists of multiple materials with distinct properties, constituting a composite structure with complex spatial architecture. This paper develops an unfitted cut finite element method (CutFEM) for transient heat conduction simulation in multi-material electronic packaging structures with complex material interfaces. In the proposed computational framework, the heterogeneous spatial geometric features of electronic packaging structures is embedded and characterized in a regular background mesh, thereby avoiding the generation of body-fitted meshes on complicated interfaces. Interface temperature and heat-flux transmission conditions are imposed by a coefficient-weighted symmetric Nitsche formulation. In addition, a ghost-penalty stabilization controls arbitrarily small intersections between the physical materials and background meshes. Furthermore, the semi-discrete and fully-discrete numerical schemes are proposed, and an explicit energy estimate is derived in detail. Finally, two-dimensional and three-dimensional electronic packaging structures containing substrate, molding compound, die, and solder balls are designed to validate its accuracy, efficiency, and scalability in simulating challenging transient thermal problems.

We introduce a continuous-time counterpart of the median filter that gradually denoises a given image. In the limit of vanishing stencil size, our results show that the evolution converges to level-set mean curvature flow. Surprisingly, and for the first time for any median-type filter scheme, we can prove the convergence of energies in this limit. Such strong convergence results are of major interest as they often have to be assumed in the literature to prove convergence of related schemes. Moreover, this result carries crucial geometric information, namely the unit multiplicity of interfaces. The proof relies on a compensated compactness argument inspired by the work of Evans and Spruck.

In many-query scenarios, data-driven surrogate models provide an efficient alternative to high-fidelity solvers for simulating physical systems governed by Partial Differential Equations (PDEs). In this context, the Latent Dynamics Network (LDNet) has recently demonstrated remarkable performance in predicting the response of spatio-temporal systems, combining Neural Ordinary Differential Equations with nonlinear dimensionality reduction. However, the original formulation assumes a fixed initial condition, limiting its applicability to many real-world applications where a system evolves from varying starting states. In this work, we overcome this limitation while keeping the end-to-end training procedure of the original LDNet and its encoder-free nature, which preserves its intrinsic independence from spatial resolution and grid topology. We infer the initial latent state directly from a small set of early-time observations, treating latent-state initialization as an adaptation problem, and investigate two strategies: an auto-decoding formulation and a meta-learning approach in which the initial latent state acts as a task-specific context variable. We demonstrate the accuracy of the proposed methods across diverse physical phenomena, spanning advection-diffusion, fluid dynamics, and solid mechanics. Meta-learning markedly accelerates latent-state inference and induces smoother, better-conditioned optimization landscapes, and spontaneously organizes the latent space into a structured representation that reflects physically meaningful features of the underlying dynamics. The coordinate-based decoder enables training from spatially subsampled data while recovering high-resolution solution fields at inference. The resulting approach provides an efficient and resolution-independent surrogate modeling framework for many-query simulations of time-dependent PDEs with varying initial conditions.

We exhibit a triangulation of an equilateral triangle for which the space of bivariate $C^1$ cubic splines has dimension larger than the dimension formula conjectured by Strang. Notably, the triangulation is such that no two edges sharing a vertex are collinear, and the triangulation is invariant under the action of the isometry group ($D_3$) of the equilateral triangle which is triangulated.

Adaptively Implicit-Explicit (AdImEx) time stepping provides stability for large time steps for mass-conservative transport, making it attractive for the numerical representation of advection in weather and climate prediction. Existing AdImEx transport schemes become first-order accurate in the large-time-step, implicit regime and may require multiple sparse matrix solutions per time step. This work introduces AdHImEx, a new AdImEx scheme that addresses both limitations simultaneously. AdHImEx provides second-order accuracy for large Courant numbers while requiring only a single matrix solution per time step. It is a Runge-Kutta scheme that blends third-order accurate explicit time stepping and Crank-Nicolson implicit time stepping to provide numerically verified stability for Courant numbers up to 100. AdHImEx is designed for flows that are predominantly explicit, activating implicit time stepping only where large Courant numbers occur. It retains the full efficiency and third-order accuracy of the explicit time stepping in regions with small Courant numbers. In this work, it is combined with a fifth-order accurate finite-volume discretisation in space for mass conservation. A second contribution of this work is a novel stage-dependent treatment of the divergence operator that preserves constancy despite different implicit and explicit Runge-Kutta stage time steps and spatially varying implicitness, removing spurious divergence. This innovation is applicable to other AdImEx time stepping schemes. Convergence, transport, and efficiency tests demonstrate substantial improvements over previous first-order AdImEx schemes, little phase and amplitude error, accurate transport on locally refined meshes with large Courant numbers, and promising efficiency characteristics for atmospheric transport.

Phase-field equations are fundamental examples of gradient-flow systems and require numerical solvers that preserve intrinsic energy dissipation. In this work, we propose a neural network-based method that combines an SAV time discretization with a mesh-free neural representation. We establish the unconditional energy stability of the underlying SAV semi-discretization, quantify the energy defect caused by solver residuals, and derive a conditional high-probability error estimate for the Allen–Cahn equation. Numerical experiments show that our approach consistently exhibits energy-stable behavior, while outperforming existing neural network-based solvers.

We propose a one-step implicit-explicit (IMEX) Runge-Kutta framework for high-order time integration of geometric flows. Starting from the Dziuk and Deckelnick-Dziuk parametric finite element formulations of curve shortening flow, we recast the spatially semidiscrete equations as systems of ordinary differential equations and construct second- and third-order schemes. Each stage requires only a linear solve. The key ingredient is a reduced-coefficient structure that allows the variational problem at each stage to be posed on the geometry generated by the preceding stage. The Runge-Kutta stages thus supply the required prediction geometry without a separate prediction procedure. The resulting schemes require no additional starting values. The framework also extends directly to the Barrett-Garcke-Nürnberg formulations of mean curvature flow and surface diffusion for curves and surfaces. Numerical experiments support the designed orders of temporal convergence and indicate that suitable coefficient choices yield mesh quality comparable to that of the corresponding first-order methods.

We construct neural approximation classes native to the graph spaces H(div) and H(curl), in two and three dimensions and, for H(div), in any dimension. Every realization lies in the space for all parameter values, and with kinked potentials, such as ReLU networks, the admissible jumps appear at finite width. The classes are images of scalar and componentwise networks under fixed operators of the de Rham complex, and do not involve a mesh or finite element emulation. For H(div) in R^n two native classes are given on an equal footing, with a skew-symmetric potential $A$: $\mathrm{Div}\,A+R_nq+\mathbf{h}$, with the divergence $q$ as an explicit unknown, and $\mathrm{Div}\,A+\mathbf{z}$ with an $H^1$ field $\mathbf{z}$; for H(curl) the analogous classes are $\mathrm{grad}\,\phi+Sr+\mathbf{h}$ in two dimensions and $\mathrm{grad}\,\phi+\mathbf{z}$ in two and three dimensions. In all of them every interface jump of the field is carried by the potential term, $\mathrm{Div}\,A$ or $\mathrm{grad}\,\phi$, while the remaining part has no interface jump (it is an $H^1$ field in the regular-decomposition classes); the classes with $\mathbf{z}$ are the componentwise approach enriched by this term. Known or learned interface geometry enters the potential through factors with trainable amplitudes, and the remaining part if the divergence jumps. The classes lead to the FOSLS-deRhaNN method, first-order system least squares with de Rham neural networks, whose loss is the least-squares functional posed in the natural spaces of the weak formulation; for elliptic equations this includes $H^{-1}$ right-hand sides and $H^{1/2}$ Dirichlet data. Elliptic equations with discontinuous coefficients and curl-curl problems are treated as instances, with the functional equivalent to the error; linear transport with discontinuous solutions and conservation laws with shocks use the same flux classes.

Splitting the exponential and exponential-like functions can significantly reduce the overall computational cost of exponential integrators, providing a strong motivation for their rigorous analysis. In this paper, we analyze splitting errors for two noncommuting unbounded operators and derive new representations that include the popular Lie–Trotter and Strang cases. Moreover, within the abstract framework of strongly continuous semigroups, we prove the convergence of a split version of the exponential Euler method and a second-order exponential Runge–Kutta method for semilinear problems. Numerical experiments on a nonlinear dispersive equation and a diffusion–reaction equation demonstrate the theoretical findings, stability properties, and efficiency of the proposed split exponential integrators compared to the state of the art.

The Multilevel Markov chain Monte Carlo algorithm accelerates Bayesian inversion by exploiting correlations between discretized quantities of interest across a hierarchy of model resolutions. A key practical parameter is the subsampling rate used to generate coarse-level proposals for the finer-level chains. Existing practice sets this rate equal to an integrated autocorrelation time as a proxy for proposal decorrelation, but this choice can be overly conservative and lead to substantial computational overhead. In this work we treat finite subsampling as an error-budget question. We derive a cost-error criterion for selecting a computationally efficient subsampling rate using estimates of sample variances, integrated autocorrelation times, and per-level costs. The criterion balances multilevel sampling errors across levels rather than enforcing approximate independence of the coarse proposals. We clarify the theoretical role of independence through an ideal reset-reference kernel, for which independent coarse proposals may be replaced by proposals generated from a reversible coarse-level kernel. We derive a Poisson-equation representation of the perturbation bias resulting from finite subsampling. Under a state-averaged lower-level forgetting assumption, this yields a conditional geometric perturbation-bias bound; in computation, the corresponding projected diagnostic is used as an empirical error-budget check. Numerical experiments on a highly autocorrelated elasticity beam benchmark show substantial reductions in computational cost compared with IAT-based subsampling. A second Darcy-flow inverse problem gives supporting evidence of qualitatively similar finite-subsampling behaviour in another elliptic benchmark.

Navier-Stokes (NS) equations lose accuracy in the transition regime when the Knudsen number is large, while the Boltzmann equation provides an adequate kinetic description with a substantially higher computational cost. In this work, a data-driven generalized Navier-Stokes (GNS) model is proposed to extend the applicability of the NS equations to the transition regime. A new relationship between the non-equilibrium variables, such as the stress tensor and the heat flux, and the equilibrium variables, including the density, macroscopic velocity, and the temperature, together with the Knudsen number, is learned from the data formed by the numerical solution to the Boltzmann equation. The standard end-to-end learning is utilized to construct the loss function, based on which, a long-term end-to-end learning method is proposed for the loss function to reduce accumulated error and improve long-term predictive accuracy. Several numerical examples, including one-dimensional wave and Riemann problems as well as two-dimensional isentropic vortex and Taylor-Green vortex problems, are studied to validate the efficiency of this new GNS model.

Adaptive multiresolution methods reduce representation cost by concentrating fine-scale degrees of freedom where needed, but their tree updates are usually governed by fixed local criteria. We introduce Learned Adaptive Multiresolution Diffusion Imaging (Learned AMDI), which preserves the AMDI fixed-tree propagator and hierarchy constraints while replacing the post-propagation selector with a shared local policy trained by proximal policy optimization. Regression tests reproduce deterministic AMDI trajectories to machine precision when identical trees are used. In the Haar implementation studied here, the deterministic one-step selector accepts no refinements in 54 decisions. Across nine held-out cases, Learned AMDI executes 393 refinements and reduces the mean terminal reference discrepancy from $0.17496$ to $0.13657$, while occupancy rises from $0.13737$ to $0.26660$. Step-resolved diagnostics reveal occasional small adaptation-energy increases; fixed-tree energy stability therefore does not guarantee monotonicity of the learned outer iteration. At comparable occupancy, a validation-tuned observed-detail threshold reaches a discrepancy of $0.13792$ with slightly better RMSE and SSIM, placing both methods on essentially the same accuracy–occupancy tradeoff. A decision-1-only control reaches $0.13742$, indicating that most of the improvement on this static benchmark arises from the initial allocation. The shared actor transfers without retraining to $64\times64$ and $128\times128$ images, improving reference discrepancy, RMSE, and SSIM relative to deterministic AMDI, while the frozen threshold rule remains competitive. Learned AMDI thus provides a hierarchy-constrained, resolution-transferable mechanism for adaptive allocation and clarifies the contribution of sequential decisions.

We study sampling discretization of the uniform norm for trigonometric polynomials with frequencies in a hyperbolic cross of level $N$ on $\T^d$. For every $d,N\ge2$ and $\varepsilon\in(0,1)$, we construct an explicit norming set with at most $\bigl(C_1(1+\varepsilon^{-1})\bigr)^{d-1}N^{1+\varepsilon}$ points and norming constant at most $\bigl(C_2(1+\varepsilon^{-1})\bigr)^{d-1}$, where $C_1,C_2$ are absolute constants. We also prove that independent Haar-distributed points achieve the same exponent $1+\varepsilon$: a sample of size at least $C(d,\varepsilon)N^{1+\varepsilon}\log(2/\eta)$ norms the entire space simultaneously with probability at least $1-\eta$, with a norming constant bounded by $\exp(C_d\varepsilon^{-d})$ and independent of $N$. The deterministic construction combines uniformly stable de la Vallée Poussin sampling operators with a Smolyak-type combination identity. The random result follows from an abstract norming theorem for sums of spaces associated with commuting regular partitions, together with a multiscale approximation of hyperbolic-cross polynomials. Finally, we show that every norming set of constant $B>1$ for the trigonometric polynomials with frequencies in $\{-1,0,1\}^d$ has at least $\bigl(\pi d/(2e\log B)\bigr)^{d/2}$ points. This rules out quadratic bounds in the dimension of the polynomial space with a prefactor growing only polynomially in $d$ when the norming constant is fixed.