Papers

Direct numerical simulation (DNS) of compressible transitional and turbulent flows requires numerical methods that combine high-order accuracy, robustness, and computational efficiency to resolve a broad range of spatial and temporal scales. This paper presents a massively parallel hybridizable discontinuous Galerkin (HDG) solver for DNS of the compressible Navier-Stokes equations on GPU-accelerated high-performance computing systems. The proposed solver combines high-order HDG discretization with robust shock capturing, diagonally implicit Runge-Kutta (DIRK) time integration, and an efficient Newton-GMRES solution strategy accelerated by additive Schwarz preconditioning and reduced-basis approximation. A distributed implementation of these methods based on GPU-aware MPI, Kokkos, and CUDA/HIP libraries enables scalable execution on heterogeneous computing platforms. The solver is demonstrated on three canonical benchmark problems covering a wide range of Mach-number flow regimes: subsonic transitional flow over the Eppler 387 airfoil, the supersonic Taylor-Green vortex, and hypersonic boundary-layer transition. Numerical results are compared with available experimental measurements and published DNS data, showing good agreement across distinct flow regimes. The results demonstrate the ability of the proposed solver to resolve laminar-turbulent transition, strong compressibility effects, shock-associated flow structures, and fully three-dimensional turbulent dynamics.

We study a PDE-constrained optimization problem for an elliptic equation with the spectral fractional Laplacian and a linear combination of Dirac measures as the forcing term; the controls are the amplitudes of these singular sources. We prove existence and uniqueness of an optimal solution and derive first-order optimality conditions. We then propose a discretization based on finite elements. Since the set of admissible controls is finite dimensional, the control variable itself does not require discretization. We conclude by deriving a priori error bounds

Dynamical low-rank approximation has become a widely used numerical method in diverse disciplines. Its main idea is to represent the matrix or tensor-valued solution to a time-dependent differential equation as a low-rank factorization. The evolution of the factorization leads to highly stiff dynamics and requires the derivation of novel time integration methods that are not prone to this stiffness. A promising family of integrators are basis-update & Galerkin (BUG) integrators as they enable implicit time integration and structure–preservation. However, current BUG integrators are limited to second–order accuracy while general-order BUG integrators are designed as projections of explicit time integration methods, thus severely limiting their use. In this work, we propose general-order BUG integrators that do not rely on an explicit time integration scheme while requiring a smaller number of basis functions to achieve high-order accuracy. We prove a general order error bound for the proposed augmented and parallel BUG integrators and demonstrate their behaviour for a series of stiff and non-stiff numerical benchmarks in which they significantly outperform previous BUG versions.

We extend the two-dimensional singular Zienkiewicz element to three dimensions, leading to a novel $H^2$-conforming finite element on tetrahedral meshes based on rational shape functions. Besides the finite element construction and its conformity analysis, we develop an exact iterative integration procedure for the associated class of rational functions. The resulting formulae allow for the exact integration of the basis functions, their derivatives, and the products occurring in finite element assembly.

Quantum time-marching algorithms for transport PDEs often represent variable coefficients and forcing through register-expanding dilations, block-encoding oracles, or repeated postselection. We present an alternative algorithm for a forced variable-coefficient advection-diffusion equation in flow-inspired skew-symmetric form that incorporates spatially varying velocity, viscous dissipation, and persistent source injection with a peak logical requirement of $n_q+1$ qubits. A centered skew-symmetric discretization makes the advection operator strictly skew-Hermitian for arbitrary velocity profiles, enabling an ancilla-free unitary realization using a Gray-code Trotter sequence of controlled-$R_y$ rotations. Diffusion is applied in the Fourier basis through a uniformly controlled rotation on one postselected ancilla, which is measured, reset, and reused between the two diffusion half-steps, while the source is incorporated classically through second-order Strang splitting. Statevector simulations for $N=16$ and $32$ recover second-order temporal convergence against high-accuracy classical solutions, while Richardson extrapolation gives fourth-order accuracy and reduces kernel calls by factors of four to fourteen. Independent tests through $N=256$ confirm second-order spatial consistency. We further show that the per-step ancilla failure probability is proportional to the instantaneous viscous dissipation rate, making postselection cost self-regulating over a fifty-fold viscosity range. Stable evolution is demonstrated for $5\times10^4$ time steps without observable secular error growth, while Gray-code advection accounts for $71$–$95\%$ of transpiled controlled-NOT gates. The fixed-width kernel provides a qubit-efficient building block for near-term hardware studies, although classical readout and state re-preparation remain the main obstacles to coherent multistep evolution.

Continuous data assimilation improves flow predictions by continually nudging a model toward available observational data. For slightly compressible flow, a recent model addresses the limitations of velocity-only nudging by assimilating both velocity and pressure data and nudging both quantities into the incompressible Navier-Stokes equations [5]; continuous-in-time error estimates and preliminary experiments show that this joint nudging is effective and substantially reduces the model error relative to velocity-only nudging. Motivated by these results, we carry out the numerical analysis of the model and its finite element discretizations. We establish stability and error estimates for the semi-discrete scheme and for the fully discrete, linearized backward Euler scheme. The analysis shows an infinite predictability horizon: the effect of the initial error decays exponentially in time, and the model error is first order in the observation resolution H and of order $\mu_1^{-1/2}$ in the pressure nudging parameter $\mu_1$. Balancing these two error terms, we choose $\mu_1=\mathcal{O}(H^{-2})$, which yields the optimal convergence rate. Numerical experiments confirm the predicted rates.

We present a finite element framework to simulate rotating mixing devices using the Mortar Element Method as a domain decomposition strategy. The model is implemented within a matrix-free Navier-Stokes framework which uses a high-order Continuous Galerkin method. The discretized domain is subdivided into rotor and stator parts. An Arbitrary Lagrangian-Eulerian approach accounts for the relative rotor-stator motion, and stabilization is ensured through the Streamline-Upwind/Petrov-Galerkin and Pressure-Stabilizing Petrov-Galerkin methods. The rotor-stator domains are connected by an interface composed of mortar cells, and continuity is weakly enforced in a Discontinuous Galerkin fashion by accounting for boundary integrals at the rotor-stator interface. Verifications of the convergence order in two-dimensional steady and transient examples report optimal rates. The geometric non-conformity created at the mortar interface due to rotor rotation does not introduce significant error in the solution. A three-dimensional example is used to investigate the model's scalability, which yields ideal strong scaling for large problems. A two-dimensional Rushton impeller example uses a torque analysis to showcase the mesh convergence, and the corresponding velocity profile is in agreement with existing numerical results. In a three-dimensional pitched blade turbine case, the power number curve (Np vs Re) shows good agreement with experimental data for Reynolds number values from 1 to 2000. An energy balance analysis reports a numerical dissipation of 1% for Re=200 and of 10% for Re=2000. By exploiting modern hardware capabilities through matrix-free methods, the proposed model is a robust, accurate, and efficient framework suitable for simulating flows with rotating geometries.

We propose a novel method that combines a multilevel decomposition with the framework of Interacting Particle Systems to compute the probability that given quantities of interest take values above or below a certain pre-specified value. Through the application of a sequential sampling scheme, the algorithm is able to achieve convergence rates faster than those achieved by the classical Multilevel Monte Carlo algorithm, and competes with the current state of the art techniques, while also being dimension independent. Furthermore, the sample points generated by the algorithm form a sequence which concentrates near the critical value of the quantity of interest. Our results are rigorously established and later verified through numerical experiments.

Operator learning in scientific machine learning is concerned with approximation of maps between infinite-dimensional function spaces; such maps frequently arise as the solution operators of partial differential equations (PDEs). Neural operators have demonstrated broad empirical success at approximating such maps from data. However, most existing neural operator architectures enforce boundary conditions indirectly through training from data even though the boundary condition is often known exactly. Furthermore, existing modifications and approaches that do enforce boundary conditions explicitly suffer from impractical restrictions, including boundary smoothness, uniform grids, and separable, box-like domains. In this work, we propose an architecture which, independently of training, satisfies homogeneous Dirichlet boundary conditions, whilst simultaneously retaining the expressivity of existing kernel-integral neural operator architectures. This is achieved by enforcing the property that the output of each layer is contained in the span of a subset of the homogeneous Dirichlet eigenfunctions of the Laplacian on the output domain. The method requires only that the output domain be bounded with Lipschitz boundary and places no restriction on the choice of discretization, making it applicable to arbitrary mesh data and general geometries. We prove universal approximation for the resulting architecture; furthermore the approach we adopt in the analysis proves universality for a broad class of kernel-integral neural operators thereby uniting existing theory for a variety of operator learning methods. We validate the proposed method on maps defined by the coefficient to solution map in 2D PDEs: Darcy flow on a square domain and the Helmholtz equation on a circular domain. Comparisons are made with alternative methods.

In this paper, we develop a deterministic, physically constrained generative framework based on the Madelung formulation of the free-particle Schrödinger equation. A reference Born probability density and a controllable initial phase function serve as initial data for the free Madelung system, which couples the Born probability density and phase function through the Bohm quantum potential, while the phase function determines the hydrodynamic velocity field. Provided that the Born probability density remains positive and the hydrodynamic velocity field generates a unique global characteristic flow, samples drawn from the reference density and transported along the characteristic flow are distributed according to the evolving Born probability density at every time. As a consequence, randomness enters only through the initial sampling; the subsequent generation is deterministic and involves neither stochastic dynamics nor an independently parameterized time-dependent velocity field. We formulate terminal-time distribution matching as a PDE-constrained phase-identification problem and derive the underlying Hamiltonian and Fisher-information structure. For isotropic Gaussian wave packets, we obtain explicit dynamics and a necessary and sufficient condition for exact reachability of isotropic Gaussian targets by quadratic initial phase functions, together with the corresponding sampling map. For a smooth prescribed potential initial velocity field, we further establish that the characteristic flow approximates the associated first-order transport map with an O(T^2) error, both uniformly and in the 1- and 2-Wasserstein distances. A numerical Gaussian benchmark validates the fully discrete forward solver, while full-grid PDE-constrained phase identification is demonstrated for asymmetric bimodal targets.

In Abner J. Salgado and Ignacio Tomas. Diagonally implicit Runge-Kutta schemes: discrete energy-balance laws and compactness properties. J. Number. Math., 31(4):313-341, 2023, the notion of $U$-stability for Diagonally Implicit Runge-Kutta (DIRK) schemes was introduced. Here we establish the equivalence between $U$- and $B$- (algebraic) stability for two- and three-stage DIRK schemes, which then} provides suitable Bochner norm a priori estimates for $B$-stable methods. As applications, we first prove the convergence of $U$-stable DIRK schemes for linear coercive evolution problems under minimal regularity via energy estimates and compactness. Second, we apply these discretizations to gradient flows, which allow us to derive discrete local energy dissipation inequalities and provide counterexamples that demonstrate the limitations of stagewise energy monotonicity.

This paper introduces a multi-class extension of a discrete-velocity kinetic traffic flow model based on a non-local Prigogine-Herman framework. We derive a hyperbolically scaled system of equations from a continuous kinetic formulation describing interactions between different vehicle classes through braking and relaxation terms. The model is then discretized with respect to the velocity variable for an arbitrary number of vehicle classes, and the structural properties of the resulting formulation are analyzed. In particular, we prove hyperbolicity and total linear degeneracy. Due to the non-conservative structure of the model, we employ a path-conservative finite volume scheme for the numerical approximation of the system. Finally, we derive the corresponding diffusively-corrected macroscopic multi-class model, investigate its stability and present numerical simulations on a single-lane road to illustrate the theoretical findings.

This work develops a regularization analysis of non-accelerated and accelerated primal-dual methods for solving linear inverse problems in the presence of noisy data. We investigate a Condat-Vũ algorithm and an accelerated primal-dual hybrid gradient method in Hilbert spaces, with focus on quantifying the effect of data perturbations on the reconstruction error. For the non-accelerated scheme, we derive error estimates in terms of Bregman distances, whereas for the accelerated scheme we establish error estimates in norm. The study accommodates a general class of convex data fidelities satisfying suitable perturbation conditions, which are verified explicitly for equality constrained and Morozov regularization. For the non-accelerated method, the analysis is further extended to Banach spaces, taking into account non-Euclidean geometries and including a particular non-reflexive setting tailored to nonnegative solution reconstruction. The results recover known behavior in classical settings while extending the regularization analysis to these more general frameworks. Numerical experiments with representative regularizers, including sparsity and entropic models, support the theoretical findings and illustrate practical performance under noise.

This paper considers the local convergence for solving general nonlinear equations. We establish randomized quasi-Gauss–Newton methods based on the approximation of the Gram matrix, addressing both the underdetermined and the overdetermined settings. For the underdetermined case, we show that our methods achieve local superlinear convergence to the optimal solution. For the overdetermined case, we show that our methods achieve local condition-number-free convergence to the stationary point of the nonlinear least-square formulation of the problem. In contrast, existing quasi-Newton methods mostly focus on square systems and additionally require the initial Jacobian estimate to be sufficiently close to the exact one.

We study exact low-dimensional reformulations for a regularized spectral approximation problem and its weighted variant. In the unweighted setting, weak-majorization monotonicity of the fidelity term enables an exact reformulation over singular values. In the weighted setting, we establish such a reformulation under compatibility and simultaneous diagonalization conditions. We further show, via a counterexample, that this exact reformulation can fail in the absence of simultaneous diagonalization.

Quantum computing offers potential advantages for solving partial differential equations (PDEs). However, most existing quantum PDE solvers primarily focus on preparing quantum states for solutions, while the efficient recovery of classical information from these states remains less explored. Motivated by the readout limitation, we propose a quantum-classical hybrid framework for Bayesian PDE inversion problems: The quantum processor evolves the PDE and evaluate the loss function with sampling noises, while the classical computer tunes the hyper-parameters in the Gaussian Process Regression to explore the next trial candidate. To match the quantum solvers for linear and semi-linear autonomous evolution PDEs, we suggest to use a normalized quantum-state loss as the data-misfit function and evaluate the new misfit by combining quantum PDE solvers with the Hadamard test, thereby allowing us to extract useful classical information using only a limited number of quantum state copies without reconstructing the full solution vector. The analysis of error propagation and overall complexity of loss evaluation under a prescribed accuracy shows that the new data-misfit function outperforms the conventional L2-loss under quantum measurements. Quantum circuit simulations of 1D and 2D linear convection diffusion equations under approximate and finite sampling loss evaluations, together with classical numerical experiments on a nonlinear forced viscous Burgers equation, demonstrate the feasibility of the proposed approach for parameter inversion even when the loss evaluations are affected by sampling noise. This framework may provide a viable quantum-assisted scheme for PDE-based inverse problems and elucidate the potential of quantum PDE algorithms in addressing a complete quantum-to-end optimization stack.

Kernels measure similarity or correlation in tasks such as regression and classification. The Gaussian kernel, other names of which include squared exponential and radial basis function kernel, is one of the most popular in Gaussian process regression. We argue that the Gaussian kernel is best avoided and should never be used as a default. The argument rests on two results demonstrating that the Gaussian kernel is extremely brittle. First, the Gaussian kernel gives rise to a conditional variance that is unrealistically small. If the variance is used to quantify predictive uncertainty, catastrophic overconfidence is almost inevitable. Second, a small variance goes hand in hand with numerical ill-conditioning, so that to use the Gaussian kernel in practice requires tricks such as nugget terms that effectively modify the underlying regression or classification model. These problems are caused by the unnatural smoothness of the Gaussian kernel, a fact we are far from the first to take notice of. The problem is not the Gaussian form itself but the analyticity of the kernel: Our argument is more broadly that analytic kernels are best avoided. For stationary kernels analyticity is essentially equivalent to an exponential decay of the spectral density.

Inverse inequalities are an important tool in the stability and convergence analysis of kernel approximation methods. In a multiscale setting, however, the trial space changes with the kernel scale $\delta$, and inverse estimates for a fixed kernel are not sufficient. The constants must remain controlled as $\delta\to0$. In this paper, we study scale-uniform inverse inequalities for spaces generated by scaled positive definite kernels whose native spaces are Sobolev spaces. We first establish inverse estimates in scale-dependent Sobolev spaces. On bounded domains, this leads to a scale-uniform inverse inequality in standard Sobolev norms with $L^2$ as the weaker norm. The estimate requires only that the separation distance of the centers be bounded above by a fixed multiple of the kernel scale. This allows the kernel scale to decrease more slowly than the separation distance, as is relevant in multiscale refinement. We then establish more general scale-uniform Bernstein inequalities on the whole space. Our result covers a broad range of weaker and stronger Sobolev indices under the same relation between the separation distance and the kernel scale. The proof works directly with the Fourier representation of functions in the kernel space and combines a low-high frequency decomposition with frame estimates for separated exponential polynomials. This avoids the additional scale factors that arise from separate comparisons of scaled and standard Sobolev norms.

Peer two-step methods are attractive for the numerical solution of stiff initial value problems because they combine favorable features of Runge-Kutta and multi-step methods, including high-order accuracy, strong stability properties, and the avoidance of order reduction. This paper develops rigorous norm estimates for the propagation operators of variable-stepsize diagonally-implicit $L(\alpha)$-stable Peer methods applied to stiff semi-linear systems. The analysis relies on a class of Peer methods whose propagation matrices possess a grid-independent left eigenvector, enabling the construction of suitable matrix norms for entire coefficient families. We show that this class contains $L(\alpha)$-stable methods of stage order $q=s$ and super-convergence of order $s+1$ with $s$ stages, exceeding the classical stage-order barrier $q\le 2$ of irreducible $s$-stage diagonally-implicit Runge-Kutta methods, and construct such methods up to order 5. Achieving high stage order is one of the principal advantages of Peer methods, especially for stiff differential equations where it helps mitigate order reduction. The resulting theory provides uniform stability bounds for stiff semi-linear problems on variable grids and establishes a foundation for further applications to optimal control and related boundary value problems. Comparative numerical results are presented for three classical stiff benchmark problems, including three novel Peer two-step methods of order 3, 4, 5, and four singly diagonally implicit Runge-Kutta methods with a first explicit step of stage order 2 and order 3, 4, 5, and 6.

Lossy compression is essential for managing massive scientific data, but per-element error bounds do not translate into bounds on downstream quantities of interest (QoIs) such as regional averages, neural network predictions, or multi-field derived quantities. We present TOPIQ, a statistical error-propagation framework that predicts QoI-level bias and uncertainty from compact compression metadata (less than 0.1% of original data). TOPIQ decomposes QoIs into primitive operators with closed-form propagation rules accounting for spatial error correlation and data-error coupling; new QoIs are supported by composition at runtime with no per-QoI derivation or retraining. Across 552 evaluations spanning 4 datasets, 3 compressors, 4 QoI families, and 8 error bounds, 93.1% of configurations achieve well-calibrated predictions. Pre-computed metadata enables post-hoc uncertainty quantification for arbitrary query regions at 56x-402x speedup over direct computation. A case study demonstrates integration into an AI-driven analysis pipeline with end-to-end confidence intervals for dynamically composed queries.

In this work, we study guided modes in systems of high-contrast resonators near nonzero interior Neumann frequencies, beyond the subwavelength regime. In the regular exterior regime, where the exterior Dirichlet problem is well-posed at the reference wavenumber, we introduce an infinite-dimensional frequency-dependent capacitance operator obtained by compressing the exterior Helmholtz Dirichlet-to-Neumann map to the traces of the interior resonant Neumann eigenspaces. We prove the norm-resolvent convergence of the continuous problem to this discrete effective operator as the contrast $\delta\to0$, and derive first-order asymptotic formulas for compact-defect frequencies and line-defect band functions. We then establish exponential off-diagonal decay of the capacitance coefficients by a Combes–Thomas argument, yielding an exponentially accurate truncation of the discrete operator, and show that its retained coefficients can be computed from local Helmholtz problems. At the physical frequency, this local approximation converges exponentially under a uniform stability assumption for the growing finite-cluster problems. The stability assumption can be removed by introducing a vanishing complex absorption together with a Hermitian symmetrization. In particular, an absorption parameter of order $\sqrt\delta$, together with interaction truncation and patch radii of order $|\log\delta|$, suffices to preserve the $O(\delta^2)$ accuracy of the first-order high-contrast expansion of the defect eigenfrequencies, yielding a fast computational method. Numerical experiments for dipole and quadrupole resonances illustrate the accuracy, exponential locality, and applicability of the discrete model to straight and bent waveguides generated by material or geometric detuning.

Multiple scattering is a fundamental wave interaction phenomenon in acoustics and electromagnetics. The multipole expansion method (MEM) is the basis of many fast algorithms, such as the fast multipole method (FMM), for such problems. However, due to the infinitely many wave reflections involved, its convergence in three dimensions remains unexplored. In this paper, we prove spectral convergence of the MEM for time-harmonic acoustic scattering by finitely many well-separated spheres in three dimensions. Using a diagonally preconditioned single-layer formulation, we analyze the degree-$N$ truncated system in a natural spherical harmonic energy space. We split the interaction truncation into target-side and source-side high-degree parts and choose a different representation for each: a projected Green-kernel representation for the former and degree-wise estimates of a translated spherical wave family for the latter. The resulting argument, based on Parseval's identity and the spherical harmonic addition theorem, provides a general framework for the convergence analysis of MEM and reveals the geometric and physical origins of the convergence factors. We also obtain a sharper estimate through the first-transfer analysis. Numerical experiments confirm the predicted spectral decay and the geometric convergence factor. This paves the way for the convergence analysis of a large class of fast algorithms for multiple scattering.

This paper proposes a fast solver for time-harmonic elastic scattering by multiple particles embedded in layered media, with either Dirichlet or Neumann boundary conditions imposed on the particle surfaces. Such problems arise in many important applications, including composite material optimization, nondestructive testing, and subsurface imaging. They are computationally challenging because of strong multiple scattering interactions among the particles and the layered interface. The proposed method represents the layered medium contribution through Sommerfeld integrals and couples this representation with a well-posed boundary integral formulation for the particle scattering problem. High-order integral equation discretization and scattering matrix are used to handle particles of general shape, while multiple scattering theory provides an efficient description of particle interactions. To reduce the cost for large particle systems, the resulting multiple scattering computation is further accelerated by the fast multipole method. The main formulation is developed in both two and three dimensions. Numerical experiments for rigid and traction-free particles validate the accuracy of the formulation, and demonstrate its flexibility in both direct scattering simulations and inverse scattering applications.

Sea ice dynamics plays a crucial role in the Earth's climate system, making it an important component of weather and climate prediction models. Its numerical simulation remains challenging, however, as it exhibits complex mechanical behaviors due to a nonlinear ice rheology and interactions with external physical forcings. In particular, sea ice presents linear kinematic features (LKFs), i.e., narrow bands of intense deformation associated with processes such as lead opening or pressure-ridge formation. Accurately representing these features is necessary as they affect thermodynamics and ocean-atmosphere exchange. Yet their number, localization, and structure are highly sensitive to spatial resolution and to the chosen discretization of the sea ice velocity. In this work, we investigate the spatial discretization of Hibler's viscous plastic sea-ice model using a fully discontinuous Galerkin (DG) representation of all variables, including the velocity field. The ability of DG elements to represent discontinuities while having high order local polynomial approximation makes them well suited for resolving the complex ice deformation features, and insuring robustness towards mesh-induced numerical artefacts. We assess the fully DG method using an established sea ice dynamics benchmark and compare the obtained sea ice deformation with state-of-the-art discretizations. Additionally, we study the theoretical convergence of the modified Elastic-Viscous-Plastic (mEVP) formulation as a pseudo-time iterative solver for the viscous-plastic (VP) momentum equations. We prove the convergence of the underlying continuous pseudo-time dynamical system towards the viscous plastic formal limit, thereby providing a theoretical foundation for the use of mEVP as an iterative solver for Hibler's VP model.

We develop a linear, unconditionally energy-dissipative, mass-lumped finite element method for the highly nonlinear Landau–Lifshitz–Gilbert (LLG) equation on quasi-uniform triangular meshes. The method is built on a projection strategy for enforcing the nonconvex pointwise constraint $|\mathbf{m}| = 1$, whose simultaneous preservation with unconditional energy stability remains challenging for standard finite element discretizations. The key innovation is a unified hybrid finite element-finite difference framework that underlies both the design and the analysis of the proposed method. In the scheme construction, we exploit the weak formulation and nodal structure of mass-lumped finite element method, while incorporating suitable interpolation operators and a node-wise length-preserving mechanism inspired by finite difference discretizations. This combination yields a linear scheme that preserves the node-wise unit-length constraint and satisfies a discrete energy dissipation law. The same hybrid framework also plays a central role in the error analysis, where the weak formulation and quasi-uniform mesh structure of finite elements are combined with interpolation-based and nodewise finite difference method to control the strongly nonlinear damping term and to establish an optimal-order error estimate. More importantly, the proposed method provides a unified framework that systematically integrates the geometric flexibility of finite element method with the pointwise constraint-preserving property of finite difference method, and thus offers a general strategy for designing and analyzing structure-preserving discretizations of constrained dissipative systems. Numerical experiments, including a classical blow-up simulation, confirm the predicted accuracy, energy dissipation, and robustness of the method.

We present a reduced-order modeling framework based on proper orthogonal decomposition for dynamic linear viscoelasticity governed by generalized Maxwell models. The hereditary constitutive law is reformulated using internal variables, and the resulting problem is discretized in space by a symmetric interior penalty discontinuous Galerkin method and in time by either the backward Euler or Crank–Nicolson scheme. The reduced model is obtained by projecting the fully discrete system onto low-dimensional spaces constructed from displacement snapshots. For the full-order discontinuous Galerkin formulation, we show well-posedness and derive an a priori stability bound without using Gronwall's inequality. Consequently, the bound grows at most linearly, rather than exponentially, with the final time. We further derive spatial error estimates and an error decomposition that separates the effects of spatial discretization, temporal discretization, and proper orthogonal decomposition truncation. The reduced formulation replaces the large full-order systems with low-dimensional problems, enabling efficient long-time simulations. Numerical experiments verify the predicted convergence behavior, assess different inner products for constructing the reduced basis, and demonstrate accurate recovery of transient oscillations and viscoelastic relaxation with substantial computational savings, even after accounting for snapshot generation and basis construction.

While Physics-Informed Neural Networks (PINNs) have emerged as a transformative paradigm for solving complex differential equations, their reliance on backpropagation-based gradient descent and automatic differentiation (AD) imposes significant computational bottlenecks and severe non-convex optimization challenges. To overcome these fundamental limitations, we propose the Physics-Informed Stochastic Configuration Machine (PI-SCM), a novel backpropagation-free framework for both forward and inverse problems in differential equations. The core mathematical contribution lies in the analytical evaluation of local Jacobians for nonlinear differential operators, which facilitates a linearized representation of the physical loss and projects it into a unified, linearized algebraic subspace. This reformulation allows for the explicit determination of optimal network weights via a sequence of generalized linear least squares solvers, effectively bypassing the iterative traps of traditional nonlinear optimizers. We develop a progressive algorithmic suite comprising localized construction (PI-SC-I), sliding-window updating (PI-SC-II), and global updating (PI-SC-III), and rigorously establish their universal approximation properties. Extensive experiments demonstrate that PI-SCM achieves high-fidelity predictive accuracy and robust parameter identification while accelerating the training process by orders of magnitude compared to standard PINNs. Our work provides a highly efficient and scalable foundation for next-generation, real-time Scientific Machine Learning applications.

We establish high-probability bounds for mixed input derivatives of wide random neural networks whose activation derivatives satisfy a factorial growth bound. Our main result specializes these estimates to $\tanh$ networks with Xavier initialization. A direct deterministic analysis based on Euclidean operator norms of the weight matrices yields derivative bounds that generally grow exponentially with the depth. We show that this growth can be substantially improved for sufficiently wide Gaussian networks by isolating the term that is linear in the highest-order derivative and controlling the corresponding tangent directions by measurable finite nets. For scalar-output $\tanh$ networks with Gaussian weights and Xavier initialization, we prove that there exist constants $C,C_0,C_1>0$ such that, whenever the common hidden width satisfies $n \geq C\left(L^3n_0^2(1+\log n_0)+L^2\left(1+\log(L/\eta)\right)\right)$, then, with probability at least $1-\eta$, the estimate $\left|D^u\mathcal{R}_{\Phi^{(L)}}(x)\right| \leq C_0 |u|! (C_1L)^{|u|-1}\prod_{j\in u}\beta_j(\eta,n_0)$ holds simultaneously for every non-empty $u\subseteq[n_0]$ and every $x\in[0,1]^{n_0}$. Thus, the first-order derivative bound is independent of the depth, while a square-free mixed derivative of order $|u|$ grows at most polynomially as $L^{|u|-1}$, apart from the coordinate factors. As consequences, we obtain high-probability bounds for the Euclidean Lipschitz constant and for weighted Sobolev norms of the network realization. The latter connect the derivative estimates to quasi-Monte Carlo integration and indicate how such regularity can enter the analysis of QMC-based training.

The aim of this work is to provide a rigorous mathematical analysis of a well-known SIS epidemic model with a saturating contact rate. First, we establish the global asymptotic stability (GAS) of the model's equilibria by employing a suitable Lyapunov function in combination with the Poincaré–Bendixson theorem and the Bendixson–Dulac criterion. The resulting GAS results improve upon previous findings for this SIS model and may also be applicable to its extensions with more general saturating contact rates. Second, we construct families of first- and second-order nonstandard finite difference (NSFD) schemes that preserve the positivity and asymptotic stability properties of the SIS model for arbitrary step sizes. All these schemes are formulated within Mickens' framework; however, compared with their first-order counterparts, the second-order schemes employ a more elaborate construction that combines a weighted nonlocal approximation of the right-hand side with suitably renormalized denominator functions. The weighted nonlocal approximation guarantees dynamic consistency, whereas the denominator functions ensure second-order convergence. Finally, we conduct numerical experiments to validate the theoretical results and demonstrate the advantages of the proposed NSFD schemes. The numerical results are in good agreement with the theoretical results. The approach developed in this work is applicable not only to other epidemiological systems but also, more generally, to mathematical models arising in a wide range of real-world applications.

Hydrodynamic ram (HRAM) loading remains a persistent challenge in impact mechanics, less due to exotic physics than to the difficulty of achieving convincing quantitative agreement in reduced-order representations. Prevailing studies typically employ a single mesh and boundary treatment validated against one experimental dataset, leaving unresolved whether agreement reflects genuine fidelity or compensating discretization errors. This work presents a planar 2D, GPU-accelerated WCSPH solver employing an Adami-type ghost-particle boundary condition, extended for the first time to a moving, decelerating disk penetrating a confined liquid channel at ballistic velocity. The boundary condition is examined against a purpose-derived potential-flow added-mass solution and the analytical acoustic reflection coefficient at the steel-water interface; this exploration reveals that the stiff equation of state and cavitation cutoff impose identifiable delineation of the solver's operating regime on quantitative agreement. A three-point resolution convergence study reveals a non-monotone near-probe pressure peak, attributed to a bounded ghost-fluid density inconsistency. A comprehensive coefficient sweep, conservation diagnostics, and O(N) GPU throughput scaling enable exhaustive characterization impractical in 3D. Qualitative comparison against literature evidence at 900 and 600m/s clarifies load bearing versus approximate elements. Notably, energy diagnostics reveal a modest total energy growth over the simulated window, with the fluid absorbing disproportionately more energy than the projectile relinquishes - an instructive signature of prescribed, non-back-reacting projectile kinematics that helps map the formulation's operating envelope. The resulting solver, while confined to planar geometry, delivers verified, convergence-checked performance well suited for rapid parametric design exploration.