Papers

We establish a discrete $H^1$-Korn inequality for finite element spaces satisfying a second order weak continuity condition. This implies that classical $H^2$-nonconforming elements satisfy the discrete Korn inequality. Based on this result, we develop a stabilized Morley finite element method for linear strain gradient elasticity using only piecewise quadratic polynomials. Thanks to the discrete $H^1$-Korn inequality, this method does not require any penalty parameter. We prove that the method is robust with respect to both $\lambda$ and $\iota$. Numerical experiments confirm the theoretical results.

Autoencoders (AEs) have become a dominant approach to nonlinear latent-space construction in data-driven reduced-order modelling (ROM), with their decoders lifting latent representations back to the ambient state space. Their prominence, however, has overshadowed an established alternative: manifold-learning methods grounded in classical numerical analysis. We revisit this alternative using Parsimonious Diffusion Maps (PDMs), benchmarking them against Proper Orthogonal Decomposition (POD)-based ROMs and several convolutional AE architectures for the two-dimensional incompressible flow past a rotating cylinder —a bifurcating Navier-Stokes (NS) system organized by a codimension-2 Bogdanov-Takens point and its associated Hopf, saddle-node, and homoclinic bifurcations. Our approach uses PDMs to identify a parsimonious and interpretable set of intrinsic latent coordinates and to estimate their dimension directly from data. Gaussian process regression then learns the latent dynamics, while convex K-nearest-neighbor (K-NN) interpolation in PDMs space constructs the pre-image map, for which we establish pointwise consistency. The resulting nonlinear ROM substantially outperforms POD-based ROMs and achieves reconstruction and prediction accuracy comparable to —and, in some bifurcating regimes, better than— that of AE-based ROMs. At the same time, latent-variable learning with PDMs requires orders of magnitude less computational time than AE training.

Reduced-order models (ROMs) are widely used to accelerate the solution of parametrized partial differential equations, with projection-based methods relying on a linear representation of the governing equations in a low-dimensional space. Their efficiency, however, is fundamentally limited by the Kolmogorov N-width of the solution manifold: when the N-width decays slowly, accurate approximations require a large number of linear modes. This limitation has motivated nonlinear ROMs, which typically replace the linear representation of the solution field with nonlinear latent representations learned through autoencoder architectures. In this work, we propose a different perspective by introducing nonlinearity directly into the representation of the governing operators. Rather than compressing only the solution field, we learn nonlinear low-dimensional representations of the full-order operators while retaining an explicit reduced system of equations. Convolutional architectures are employed to exploit the locality and spatial structure of differential operators, and the resulting compressed operators are coupled with a nonlinear representation of the solution through a training strategy that accounts for both reconstruction accuracy and the error of the reduced-order solution. To address large sparse operators arising from fine and unstructured discretizations, continuous convolutions are used to operate directly on their non-zero entries. Finally, radial basis function interpolation enables the prediction of the compressed operators at unseen parameter configurations, providing an efficient online stage without assembling the corresponding full-order system. Numerical experiments on three parametrized problems demonstrate the potential of nonlinear operator compression as an alternative strategy for equation-based ROMs.

This paper proposes a novel discretely entropy conserving/stable flux reconstruction (FR) scheme for nonlinear conservation laws. The scheme combines a hybridized split-form differential operator with a generalized entropy projection, enabling discrete entropy conservation for arbitrary nodal distributions and arbitrary correction functions without modifying the original quadrature weights. The formulation is built directly upon the classical FR framework and is shown to preserve conservation properties while achieving entropy conservation in the semi-discrete sense. Numerical experiments using the isentropic vortex problem confirm the theoretically predicted accuracy and demonstrate robust entropy behavior with entropy-dissipative interface fluxes.

This paper develops a staggered discontinuous Galerkin (SDG) scheme based on Cartesian grids for Stokes equations that is simple to implement, intrinsically pressure-robust, and superconvergent for all variables. Instead of the composite meshes used in standard SDG, we construct staggered quadrilateral meshes from Cartesian grids. The scheme takes velocity, pressure, and velocity gradient as unknowns, and employs piecewise-constant spaces with carefully designed staggered continuity. The additional gradient unknowns are locally eliminated via static condensation and can be further removed by mass lumping without loss of accuracy. An explicit pointwise formulation of the scheme is derived, which facilitates implementation and enables a detailed pointwise analysis. We rigorously prove pressure robustness and second-order superconvergence, which holds on general non-uniform Cartesian grids. The scheme is further extended to Navier-Stokes equations by introducing a novel discrete convection term with second-order consistency. Combined with the scalar auxiliary variable (SAV) approach and the Crank-Nicolson (CN) scheme, this yields an unconditionally energy-stable and second-order accurate scheme. Numerical experiments validate the theory and demonstrate accuracy and robustness.

In this paper, we develop a preconditioner that is robust with respect to both the mesh size and the interaction horizon for finite element discretizations of volume-constrained nonlocal diffusion problems. Motivated by the Fourier multiplier of the nonlocal operator, we construct the preconditioner $B_h=G_h+\beta\delta^2D_h^{-1}$, where $G_h$ is a symmetric local Poisson solver, $D_h$ is the lumped mass matrix, and $\beta$ is a kernel-dependent constant. We rigorously prove that the spectrum of the $B_h$-preconditioned system is bounded above and away from zero uniformly with respect to the mesh size and horizon. Moreover, we develop a pair-cell assembly for the finite element matrix that preserves the symmetry and positive definiteness required by the preconditioner and PCG, and establish an $L^2$-error estimate for the finite element discretization. Finally, we use a single symmetric geometric multigrid V-cycle for $G_h$ in the numerical experiments, and the corresponding numerical results confirm the predicted convergence rate and show bounded condition numbers and PCG iteration counts over the tested meshes, horizons, and kernels.

Point cloud-based CFD enables flow analysis directly on discrete points obtained from 3D scanning and medical imaging, bypassing surface reconstruction, geometry cleanup, and boundary-fitted mesh generation. Derived from immersogeometric analysis, the method immerses the point cloud in a background mesh and enforces no-slip conditions on discrete points through a Nitsche-based weak boundary condition (BC). The framework delivers accurate velocity fields, pressure distribution, and integrated loads; however, accurate prediction of the local wall shear stress (WSS) has remained a critical challenge. The geometry intersects the background mesh arbitrarily, producing cut elements that lack the regularity required for consistent gradient evaluation. The issue is compounded by the stabilization term of the weak BC, whose parameter estimation in the symmetric Nitsche formulation is dependent on the cut configuration and affects the variationally consistent definition of traction from which the WSS is computed. In this work, we propose a new method to obtain accurate wall shear stress in immersed flow analysis with application to point cloud-based CFD, using a non-symmetric Nitsche's formulation with near-wall modeling and a patch-based stress recovery approach with traction compatibility. The method is validated on canonical benchmarks and applied to turbulent flow past a sphere and to a patient-specific aorta, showcasing excellent agreement with reference results.

Neural operators are data-driven models that learn mappings from inputs that parameterize partial differential equations, such as spatially varying coefficients, initial conditions, forcing terms, boundary conditions, or geometries, to solution fields or quantities of interest. Once trained, they can serve as surrogates for classical numerical solvers in many-query settings that require repeated evaluations for varying inputs. We address the question of when, and then why, neural operator surrogates outperform classical numerical solvers, in terms of cost for a given accuracy. We focus on the post-training, many-query limit, in which data-acquisition and training costs are treated as fixed and fully amortized. Even in this deliberately favorable regime for neural operators, there are regimes in which classical solvers outperform the surrogate models. We compare the cost-accuracy performance of neural operator surrogates and classical numerical solvers through a reproducible benchmark study comparing neural operators with problem-matched classical solvers on representative problems in computational science and engineering, focusing on prediction error, per-query floating-point cost, and wall-clock runtime. Neural operators are most competitive at low-to-moderate accuracy requirements. Their floating-point cost advantage depends strongly on the problem structure, arising when they avoid temporal or nonlinear iterations or predict a reduced quantity of interest rather than a full solution field. Additional wall-clock speedups result from dense tensor operations that are well suited to modern hardware. As the target accuracy is tightened, achieving the required accuracy with neural operators becomes increasingly challenging, and classical solvers outperform surrogates in this regime; thus classical solvers will remain important for verification and high-accuracy computation.

Tracking neurite morphology over time helps characterize structural changes during neuronal development and deterioration, but long-term time-lapse imaging is resource-intensive and difficult to scale. Forecasting future morphology could reduce this burden. Existing neurite digital-twin models such as gated spatiotemporal attention (gSTA) produce a single deterministic forecast without representing variability among plausible futures. We introduce Morphology-Gated Residual Diffusion (MGRD), a compact stochastic surrogate that jointly forecasts twenty future neurite-morphology frames from ten observed frames while conditioning on morphology features derived from the latest observation. On controlled phase-field trajectories, MGRD reduces trajectory-wise mean MAE by 9.7% relative to a matched control while updating 4.46 times fewer parameters. On human iPSC-derived neuron microscopy, MGRD improves all four reported metrics over gSTA, including a 39.6% reduction in trajectory-wise mean MAE and a 45.3% increase in skeleton F1. Without mouse-domain retraining or fine-tuning, MGRD also improves MAE and skeleton F1 on mouse cortical-neurosphere microscopy across 10-40-min sampling intervals and forecast horizons beyond 13 hours. Repeated sampling provides a case-level variance score for ranking forecast difficulty. Retaining approximately 60% of the lowest-variance cases reduces mean MAE by 17.6% on iPSC microscopy and 16.8% on simulation data. MGRD uses 1.01% of gSTA's parameters, requires less than one tenth of its training-update time, and generates a 50-step DDIM trajectory 7.9% faster when morphology features are cached. These results establish MGRD as a compact stochastic surrogate for neurite-morphology forecasting and case prioritization across simulation and microscopy datasets.

This paper establishes time-uniform accuracy guarantees for deterministic square-root ensemble Kalman filters in linear–Gaussian state-space models under perfect-model dynamics. For hyperbolic, detectable systems, we prove that the ensemble means and covariances approximate their Kalman filter counterparts uniformly in time, with high probability and without inflation or resampling. The required ensemble size depends on the effective rank of the initial covariance and the dimension of the unstable subspace, rather than on the ambient state dimension. We also analyze a localized square-root ensemble Kalman filter for weakly coupled spatial systems. We show that the stabilizing Kalman covariance inherits spatial decay from the dynamics and derive accuracy bounds that separate sampling error from localization bias. For the localized filter, the required ensemble size depends on the local intrinsic dimension and only logarithmically on the number of spatial blocks, while the localization bias scales with the interaction strength. Our analysis combines stability and forgetting for possibly singular Riccati recursions, nonasymptotic covariance concentration, and perturbation estimates for localized covariance dynamics.

We consider the problem of numerical differentiation of functions defined on the half-line: the derivative $f^{(r)}$, $ r=1,2,\ldots$, is recovered from a finite set of perturbed Fourier–Laguerre coefficients of a function $f$, whose error is measured in $\ell_p$, while the accuracy of approximation is measured in the uniform metric with the weight $t^{\alpha}e^{-t}$. It is established that the Wiener class $W^{\mu}_{s}$ admits such a setting for $\mu>r+1-1/s$. A class of methods $\mathcal{S}^{\theta}$ is constructed that realize the optimal order of accuracy of numerical differentiation on $W^{\mu}_{s}$ and use the smallest amount of input information in order.

Low-rank approximation is a fundamental tool for scalable matrix computations. While such approximations have classically been formed via the truncated SVD, recent advances in randomized numerical linear algebra have produced methods of comparable accuracy at a fraction of the cost. A key advantage of these approaches is their ability to operate with limited access to the full matrix. A prominent example is the CUR decomposition, which builds an approximation from only a subset of columns and rows, making it particularly well-suited to settings where global matrix access is unavailable, such as in the design of spectro-microscopy experiments where even matrix-vector products cannot be computed. However, the randomness inherent to these methods introduces a key challenge: efficiently assessing approximation accuracy under the stringent constraint that only a subset of matrix entries can be observed. We derive a "blind" estimator of the Frobenius error of a CUR decomposition. For this estimator, we derive a worst-case lower bound on the number of required column queries, and provide techniques for quantifying its uncertainty. Finally, we demonstrate the performance of the estimator and the uncertainty sets on a mix of synthetic problems and real spectro-microscopy examples.

In undersampled inverse problems (such as sparse-view computed tomography), only a small number of measurements are collected, which reduces radiation exposure and acquisition time and cost, and can also address the inaccessibility of certain acquisition arrangements. Most learning-based methods for such problems require supervision in the form of fully sampled measurements and ground-truth images, which are costly or even infeasible to acquire. To overcome this issue, we propose Exact Learning for Inverse Problems from Partial Self-supervision (ELIPPS), an incomplete self-supervised training paradigm for undersampled inverse problems that requires neither ground-truth images nor fully sampled measurements. ELIPPS learns solely from incomplete forward measurements on a fixed incomplete supervision set. Our theory shows that when the data distribution is invariant under certain transformations, minimizing a masked empirical risk is equivalent to minimizing the full self-supervised risk, i.e. the risk against the complete, noise-free measurement. The equivalence is exact for arbitrary equivariant hypothesis classes and holds for signal-dependent noise such as the pre-log Poisson statistics of low-dose tomography, and for the log-transformed count model up to a quantifiable bias. When the underlying coverage condition is only approximately satisfied on a discrete grid, we give a stability estimate in terms of the associated frame constants and the irreducible error of the inverse problem. We realize ELIPPS for computed tomography by exploiting rotation and reflection invariance. In our experiments, ELIPPS substantially outperforms naive masked supervision and reaches the same order of accuracy as a reference model trained with full clean measurements.

Rigorous error analysis for numerical approximations of parabolic equations with hysteresis remains limited, even for the widely used Prandtl–Ishlinskii model of play type. In this work, we establish an $O(h+\tau)$ error bound for an implicit Euler $P_1$ finite element discretization. The analysis requires neither higher-order temporal regularity nor additional spatial regularity of the hysteresis variables. For the temporal discretization, we exploit a convex subgradient-flow structure in a weighted Hilbert space and estimates of residual monotone terms to obtain first-order convergence. For the spatial discretization, only the diffusive field is restricted to the finite element space, and a constraint-preserving comparison yields an $O(h)$ semidiscrete estimate. The analysis is developed for play-type Prandtl–Ishlinskii operators formulated directly on a spatial Hilbert space, encompassing the canonical pointwise extension as well as more general spatially structured constraints.

One novel monotonous finite volume element (MFVE) scheme is put forward for convection dominant diffusion problem. The main contributions of this paper include four aspects. Firstly, one upwind sub-control volume, called upwind volume, is introduced for the discretization of the convection term similar to the upwind element, which leads to the upwind property. Secondly, one first-order linear discrete operator for the convection term in the balance equation is directly obtained by numerical integration in upwind volume, which is different from the classical finite volume methods.Thirdly, one asymptotic expansion of gradient functions in the dual element of each node is derived by one extension interpolation function. Then, one second-order nonlinear discrete operator for the convection term is constructed by the asymptotic expansion and error estimation of its FVE solution, where some perturbed coefficients are skillfully put into the expansion, which serves as a high-order correction. Fourthly, one novel two-order MFVE scheme, also accompanied by another one-order MFVE, is designed for convection-dominated diffusion problems together with one positivity-preserving finite volume element (PFVE) discrete diffusive operator. The L2 norm of the error of the approximate solution is derived. Finally, numerical results validate theoretical conclusions.

We investigate an inverse coefficient problem for an augmented Keller–Segel model arising in the description of phototactic and chemotactic population dynamics. The spatiotemporal evolution of the population density is governed by an advection-diffusion equation in which both the diffusion coefficient and the advection coefficient depend on space and time. Due to the underlying intracellular adaptation mechanism, even in a time-independent external environment, the diffusion coefficient \(D(x,t)\) and the scalar advection coefficient \(K(x,t)\) remain time-dependent and relax exponentially toward their respective steady-state profiles. Importantly, their relaxation is governed by the same rate, determined by the intracellular adaptation dynamics. From internal measurements of the density, we establish conditional Lipschitz stability for recovering either coefficient when the other is known, and derive cross-sensitivity estimates quantifying the compensation between diffusion and advection perturbations that produce the same density data. We further identify a finite-time separation between the transient and steady regimes: the full dynamics admit an exponentially accurate frozen-coefficient approximation at late times. When the steady-state coefficients are known, we also prove unique identification of the relaxation rate from the long-time behavior of the solution. Finally, we derive an explicit linearization of the forward map and, motivated by the temporal separation above, develop a two-stage gradient-based reconstruction strategy. Numerical experiments demonstrate the feasibility of the proposed method.

A new approach is developed to study the convergence of parametric finite element approximations to the mean curvature flow of closed surfaces in three-dimensional space. In this approach, the error analysis is conducted by comparing the numerical solution to a dynamic Ritz projection of the mean curvature flow introduced in this paper, rather than an interpolation of the mean curvature flow, as commonly used in the literature. The errors associated with the dynamic Ritz projection in approximating the mean curvature flow are established in the $L^2$ and $W^{1,p}$ norms. Leveraging these results, optimal-order convergence of parametric finite element methods for mean curvature flow of closed surfaces in the $L^\infty(0,T;L^2)$ norm is proved, including the convergence of parametric finite element methods with piecewise linear finite elements.

This paper introduces a novel formulation and an associated space-time finite element method for simulating solutions to the nonlinear Schrödinger equation. A major advantage of the proposed algorithm is its intrinsic ability to preserve the conservation of mass, energy, and momentum at the discrete level. This is proved for the numerical solutions determined by the fully discrete implicit scheme. An effective iterative scheme is proposed for solving the nonlinear system based on an equivalent formulation which suggests using Newton's iteration for the solution and no iteration for the Lagrange multipliers in the nonlinear system. Extensive numerical examples are provided to demonstrate the high-order convergence and effectiveness of the proposed algorithm in conserving mass, energy, and momentum in the simulation of one-dimensional Ma-solitons and bi-solitons, as well as of two-dimensional solitons governed by the nonlinear Schrödinger equation. The numerical results show that the mass-, energy- and momentum-conserving method designed in this paper also significantly reduces the errors of the numerical solutions in long-time simulations compared with methods which do not conserve these quantities.

We develop a bihomotopy argument to estimate the discrete material time derivatives of consistency functionals. Combined with summation-by-parts, this argument yields a directional super-convergence estimate for the consistency contribution when tested with the normal error velocity. Together with the super-approximation of the reversely weighted discrete normal, these techniques establish $L_t^\infty H_x^1$ super-convergence of the projection error and optimal $L_t^\infty L_x^2$ convergence of the trajectory error of a stabilized BGN method for curve-shortening flow with piecewise linear parametric finite elements.

We construct a mesh-induced sparsely connected neural network framework that exactly reproduces arbitrary-order Lagrange finite element spaces over simplicial meshes. Unlike conventional black-box neural surrogates, the proposed network architecture is fully dictated by finite element discretization rules: local computations stem from simplex geometry and barycentric coordinate transformations, while global consistency is enforced through shared degrees of freedom. For linear Lagrange elements, local basis functions are directly implemented via affine barycentric layers, and high-order polynomial bases are explicitly decomposed into barycentric product compositions realized by specially designed $\mathrm{ReLU}^p$ modules. Equipped with element indicator branches and multiplication units, these modular local components are globally assembled into a sparsely connected neural network whose function space coincides exactly with the target finite element space, thereby inheriting the complete classical finite element approximation theory. By prescribing customized backward differentiation rules for piecewise activations, function values and their spatial gradients can be simultaneously extracted via automatic differentiation within a unified computational graph, eliminating the separate gradient calculation subroutines required in standard finite element implementations. Numerical experiments verify the accuracy of the neural network representation of Lagrange finite elements. Furthermore, by virtue of the intrinsic mesh-free nature of this neural network representation, finite element functions can be interpolated between non-matching meshes, and the proposed scheme can be applied to adaptive finite element methods for solving parabolic partial differential equations. An open-source code implementation of the proposed architecture is made publicly available.

High-order isogeometric discretizations produce a small number of non-physical outlier modes at the top of the discrete spectrum. These modes can pollute the numerical solution in wave propagation and structural dynamics. Most existing remedies modify the spatial discretization through outlier-free bases, mass lumping, subspace projection, and related constructions. We introduce MIISO (Modal Integrators for Isogeometric Stabilization of Outliers), a family of exponential Rosenbrock-Krylov methods of orders two through six that suppress the outliers in the time integrator and leave the spatial discretization unchanged. We construct a non-holomorphic amplification map by blending a Padé approximation that preserves physical oscillatory modes with a Pade approximation that damps high frequencies, with the blend placed at the known spectral location of the outliers. For orders three and higher the damping approximant is repaired so that the method is unconditionally A-stable. Physical modes receive no numerical dissipation and the outliers are driven to zero. The outlier handling and the method steps are fully matrix-free, and the method extends to nonlinear problems. Numerical experiments recover the predicted orders, match spatial outlier removal without changing the discretization, show lower cost than generalized-alpha at matched accuracy and in scaling, and demonstrate the method on nonlinear and three-dimensional problems.

Data assimilation aims to infer the state of complex dynamical systems based on observational data. However, accurate inference of the multimodal posteriors induced by nonlinear or non-injective observation operators remains a key challenge under high-dimensional and sparse observation conditions. Ensemble filters scale to high dimensions but are confined by restrictive distributional assumptions, while training-free generative filters (e.g., EnSF, EnFF) alleviate this limitation but may introduce structural errors and hinder information propagation under sparse observations. To address these issues, we propose a training-free, asymptotically exact posterior transport method. Firstly, a covariance-aware posterior transport scheme is designed, which embeds the forecast cross-covariance into flow-based transport and accurately recovers unobserved states while preserving the non-Gaussian posterior structure. Furthermore, the method combines a tractable observation-adaptive proposal with posterior correction, ensuring accurate approximation of the nonlinear posterior distribution. Finally, we establish the corresponding posterior flow theory, from which the asymptotic exactness of the proposed method relative to finite-ensemble surrogates and the Wasserstein error bound are derived. Experiments on Double-Well, Lorenz-96, and Kolmogorov flow show that the proposed method captures complex posterior structure and remains accurate under sparse, nonlinear, and non-injective observations. In the sparse non-injective setting, it reduces RMSE by 56% relative to the best baseline.

We study numerical observables relevant to delocalization in the three-dimensional Manhattan mirror model, also known as the K-lattice. Our principal experiment separates the locations of transmitting boundary channels from the routing between them. We derive an exact conditional comparison that preserves all transmission and reflection aperture sets, and find a statistically resolved deficit of physical transmission relative to this comparison. Nearby transmitting channels are themselves suppressed relative to a uniform subset with the same cardinality. A second experiment resolves the contribution of repeated vertices by their return age. The squared norm of the signed old-memory remainder is consistent with an inverse-square-root decay over the measured range. Annular experiments quantify decreasing continuation losses, while an initially apparent history dependence fails independent replication. Finally, a periodic configuration has a connected auxiliary contact graph but only finite routing cycles. These results identify estimates on reflection routing and adaptively selected revisit charges that may assist a rigorous proof. They do not establish an infinite-volume delocalized phase.

We propose a fast deterministic particle method for the spatially homogeneous Landau equation. The method combines a particle representation of the solution with a Fourier approximation of the nonlinear collision flux. Nonuniform fast Fourier transforms are used to reconstruct the density from the particles and to evaluate the flux and density at the particle locations, while the convolutional structure of the flux enables its efficient computation by FFTs. For a fixed transform tolerance, the cost per time step is $\mathcal{O}(N+M^d\log M)$, where $N$ is the number of particles and $M$ the number of Fourier modes per velocity dimension. We establish a consistency estimate for the reconstructed velocity field on regions where the reference density is bounded away from zero. The estimate separates the spectral truncation error from the particle-density reconstruction error, showing how the latter can dominate for sufficiently smooth densities. Two-dimensional numerical experiments with Maxwellian interactions illustrate the accuracy and efficiency of the method, its sensitivity to particle and spectral resolution, and the effects of filtering. Comparisons with a direct blob implementation demonstrate improved accuracy at a strongly reduced computational cost.

The Time-Evolving Natural Gradient (TENG) method evolves neural approximations of time-dependent partial differential equations by projecting time-discrete target states onto the neural-network manifold through local natural-gradient iterations. For stiff diffusion-dominated problems, however, explicit Runge–Kutta methods face severe stability restrictions. High-frequency diffusive components can drive stage targets beyond the locally reachable region of the neural manifold, worsening the conditioning of local least-squares projections and destabilizing parameter updates. We propose TASE–TENG, which applies the highly stable explicit operator $T_p(hW_0)$ to each Runge–Kutta stage increment before projection, where $W_0$ approximates the dominant stiff diffusion operator. To construct this surrogate in a mesh-free setting, we develop a learnable operator whose local interaction coefficients are generated by a shared geometry-dependent neural rule. A structured factorization enforces discrete self-adjointness, dissipativity, and constant preservation. The operator is trained offline and remains fixed during online integration, where it is used only within the TASE stabilization operator. Numerical experiments on stiff diffusion and reaction–diffusion problems demonstrate improved stability and accuracy at the same time-step size, with reduced errors from high-frequency diffusive modes. The learned rule also allows the stabilization operator to be reconstructed at different spatial point resolutions without retraining.

We establish the asymptotic-preserving property of the exponential Euler method for stochastic reaction-diffusion-advection equations in $\mathbb R^2$ under fast Hamiltonian advection with multiple critical points. In the fast-advection limit, the dynamics reduces to a stochastic partial differential equation on a noncompact metric graph. A key ingredient is the strong convergence analysis of the limiting exponential Euler scheme, which is complicated by the graph's noncompactness, its multiple edges and vertices, and the nonuniform ellipticity and vertex degeneracy of the graph operator. Using a weighted $L^2$ framework and analytic-semigroup smoothing estimates, we prove a temporal convergence order arbitrarily close to $1/2$ for weighted $L^2$ initial data and exactly $1/2$ under a half-order fractional-domain condition. Combining these estimates with the fast-advection limits yields the asymptotic-preserving property in the multiple-critical-point setting. Numerical experiments illustrate the asymptotic behavior and confirm the predicted temporal convergence rate.

We introduce Neural $S$-Transform Reconstruction (NeSTR) for Gaussian stochastic partial differential equations (SPDEs). The learned object is not the random solution itself, but its $S$-Transform on a finite dimensional subspace of the white noise test space. For Wick type SPDEs, NeSTR yields a deterministic parametric PDE for the restricted $S$-Transform, and the stochastic solution is recovered from the Taylor coefficients of this transform at the origin. Thus the neural approximation is tied directly to the Wiener chaos expansion of the solution. The main novelty of NeSTR is an inverse $S$-Transform reconstruction principle combined with finite mode learning: we prove consistency of the recovered chaos coefficients, separate noise mode, chaos truncation, deterministic solver, and neural approximation errors, and give a deterministic learning formulation for the restricted transform. Numerical tests on an additive stochastic heat equation and a multiplicative Wick heat equation show that polynomial spectral representations in the transform variables provide stable access to higher order chaos coefficients. The resulting reconstruction separates the equation dependent deterministic coefficient fields from a reusable Gaussian–Hermite basis. Consequently, NeSTR learns an analytic generating function whose derivatives recover the Wiener chaos coefficients, rather than learning stochastic trajectories directly.

We extend the high-order Zipped Finite Element Method to three-dimensional polyhedral elements. The construction relies on a local sub-tetrahedralization obtained by triangulating the polygonal faces and connecting the resulting vertices to a suitable interior point. High-order shape functions are constructed as linear combinations of standard finite element basis functions, with coefficients chosen to ensure exact polynomial reproduction on the faces and in the element interior while preserving global C0-conformity. We further introduce a QR-based strategy as a systematic alternative to the heuristic procedure that was adopted for the two-dimensional case. The QR-strategy is more systematic and robust than the heuristic strategy, while being independent of the spatial dimension and avoiding geometric alignment tests and problem-dependent tolerances. The resulting method is formulated for a diffusion-reaction problem and numerically assessed through three-dimensional convergence tests, confirming the expected optimal polynomial order of accuracy.

We study the application of exponential time integration methods, coupled with Fourier pseudospectral space discretization, to the hyperbolic Korteweg-de Vries (KdVH) system. We investigate the asymptotic preserving (AP) properties of such discretizations, showing that, in general, Lawson methods are not asymptotic preserving, because the auxiliary (derivative-approximating) variables do not satisfy the limit equilibrium manifold, whereas exponential time differencing (ETD) methods are AP for all components, including both the solution variable and the auxiliary variables. We also present an efficient numerical implementation based on an exact formula for the matrix exponential of this system, and compare it with previously-proposed ImEx Runge-Kutta time integration, showing that exponential methods can be competitive in this context.

We study the global continuation of solution manifolds arising in dynamical systems. We present a rigorous continuation method based on a Chebyshev series expansion of the solution manifold. The branch is first approximated by a high-order Chebyshev interpolation polynomial, and an explicit error bound is then obtained by verifying the contraction of a quasi-Newton operator near this approximation. The contraction is formulated on a weighted $\ell^1$ space, giving a finer control than the typical $C^0$-error bound obtained from the uniform contraction theorem. In fact, the latter follows directly from our contraction operator. Furthermore, we discuss how our strategy applies naturally to pseudo-arclength continuation, where the continuation parameter fails to provide a valid local coordinate, and extends to multi-parameter continuation. Lastly, we detail two applications in which we compute a two-parameter family of steady-states for the Cahn–Hilliard equation, and a one-parameter family of steady-states undergoing saddle-node bifurcations for the Shigesada–Kawasaki–Teramoto system.