Papers

This paper introduces a framework that integrates the dynamic stiffness matrix (DSM) with physics-informed neural networks (PINN). The DSM-PINN embeds physical constraints within the model and demonstrates robustness, particularly when addressing limited datasets across diverse investigations. In this approach, deep neural network outputs approximate the displacement fields of element nodes. Unlike the finite element method (FEM), the element shape functions are homogeneous solutions to the governing partial differential equation, forming the basis of the exact dynamic stiffness matrix, thereby avoiding high-order derivative terms. This matrix also serves as a frequency-domain spectral element, resulting in a strong-form PINN. The loss function is produced by connecting neural networks with dynamic stiffness matrices. We focus on utilising PINNs to resolve eigenvalue problems by employing the Wittrick-Williams algorithm, which overcomes the challenge of neural networks failing to converge to higher-order eigenvalues. Additionally, the frequency domain-PINN method is used to analyse structural dynamic responses under moving and impulsive loads, addressing the limitation of neural networks in handling complex numbers. Theoretical convergence stability of the suggested approach is also analysed even DSM is an indefinite matrix after implementing the boundary condition. The numerical results validate the practicality and efficacy of the recommended approach.

Physics-Informed Neural Networks (PINNs) have recently gained considerable attention as a mesh-free framework for solving partial differential equations. Nevertheless, their performance deteriorates when applied to strongly coupled multiphysics systems, such as Biot's consolidation model, due to severely ill-conditioned optimization landscapes. In this work, we propose a robust PINN-based solver, termed FS-ENGD-PINN, which synergistically integrates physics-based decoupling with geometry-aware optimization. Specifically, the Fixed-Stress (FS) splitting scheme is employed to decompose the coupled poroelastic system into contractive mechanics and flow subproblems, thereby significantly improving training stability and convergence. To further accelerate optimization, we adopt Energy Natural Gradient Descent (ENGD), which approximates the Newton direction in function space effectively mitigates stiffness-induced slow convergence. Moreover, to address volumetric locking arising in the nearly incompressible regime, we incorporate a three-field mixed formulation with an additional total pressure variable into the PINN framework. Extensive numerical experiments demonstrate that the proposed FS-ENGD-PINN consistently outperforms standard PINN formulations in terms of accuracy and robustness, providing a unified and reliable learning-based solver for poroelasticity across a wide range of material parameters.

We study a thermodynamically consistent non-isothermal phase-field system coupling two order parameters and the inverse temperature through a fully non-diagonal Onsager mobility. The model describes the interaction of mass diffusion, heat conduction, and local phase relaxation while conserving mass and internal energy and producing entropy. Global generalised dissipative weak solutions are constructed using a fully discrete approximation. The discrete scheme conserves mass and internal energy and satisfies a discrete entropy inequality. An availability estimate, together with a dimension-adapted barrier, yields strict positivity at fixed mesh and uniform estimates. Compactness then allow passage to the continuous system. Concentration of the singular part of the internal energy is represented by a non-negative defect measure. The entropy-availability structure yields finite dissipation on the infinite time interval and the existence of stationary $\omega$-limit states. Finally, a relative-entropy argument establishes weak-strong uniqueness whenever the weak and strong solutions remain in a bounded thermodynamic state range.

This paper deals with the computation of the sum of alternating series, whose general terms can be expressed by means of analytic functions. After rewriting the series as a weighted integral with the Abel weight, we employ the sinc (trapezoidal) rule and analyze the remainder term with respect to the number of quadrature points. We provide some numerical experiments to show the reliability of the derived error estimate and to test the algorithm for automatic summation with prescribed accuracy. All the Matlab codes used in the present paper can be found as open-source software at the authors' homepage.

We present a family of second-order, linear, unconditionally stable implicit-explicit (IMEX) methods for the Cahn-Hilliard-Navier-Stokes (CHNS) equations modeling matched-density two-phase flows. The proposed semi-discrete scheme combines extrapolation of the nonlinear terms with an auxiliary-variable formulation of the nonlinear free-energy term and a temporal-curvature regularization controlled by a parameter $\epsilon$. We establish a discrete energy estimate showing unconditional long-time stability of the method for $\theta\in(1/2,1]$ and $\epsilon\geq0$. The resulting scheme requires only linear solves at each time step. Numerical experiments demonstrate approximately second-order temporal convergence and examine mass conservation, energy dissipation, numerical robustness, and several representative interfacial-flow problems, including spinodal decomposition, droplet shape relaxation, two-phase lid-driven cavity flow, and Rayleigh-Taylor instability.

We summarize our advances in the algorithmic development and hardware utilization for risk-averse optimization problems in random materials. This includes risk-averse optimization using the entropic risk measure, as well as recently developed sampling techniques for random materials, that are interoperable with the optimization framework. Furthermore, we discuss recent progress in the efficient utilization of modern hybrid hardware architectures for these methods to solve three-dimensional partial differential equations.

We present recent developments in multilevel sampling methods under resource constraints. Over the past 15 years, multilevel methods have become widely used for uncertainty quantification. However, scaling them to high-dimensional problems and high-performance computing (HPC) environments remains challenging. In this work, we discuss two algorithms designed to address these issues: the budgeted Multilevel Monte Carlo (MLMC) method and the Multilevel Stochastic Gradient Descent (MLSGD) method. We demonstrate their effectiveness on HPC systems under consideration of the available computational resources for applications in forward uncertainty quantification (UQ) and optimal control (OC) under uncertainty.

We determine, up to absolute constants, the worst-case number of point evaluations required for a weighted $L_2$ Marcinkiewicz-Zygmund inequality for an $m$-dimensional complex function space. If $0<\varepsilon<1$ is the relative distortion, this number is $$\Theta\Big(\min\Big\{m^2,\frac{m}{\varepsilon^2}\Big\}\Big),$$ and exact discretization has the sharp worst-case value $m^2$. While the upper bounds follow from recent constructions, our contribution is the construction of function spaces that are hard to discretize and yield matching lower bounds. We use a trace-variance inequality for weighted subframes of unit-norm tight frames. One such instance is the complete-graph edge frame, which yields a construction in every dimension. Singer equiangular tight frames improve the constant when $m-1$ is a prime power, while maximal equiangular tight frames give the strongest bound possible using our method whenever they exist. We also derive consequences for the conditioning of weighted least-squares systems and for standard condition-number-based iteration estimates when these systems are solved by LSQR.

We develop a neural-network-accelerated multiscale hybridizable discontinuous Galerkin method for elliptic problems with heterogeneous coefficients. The method preserves the standard MsHDG local-to-global structure: fine-scale HDG problems on coarse blocks define discrete Dirichlet-to-Neumann operators, which are assembled through the standard MsHDG global skeleton equations. To reduce the cost of constructing these local operators, we train a neural network on coefficient fields defined on a reference block and use the predicted operators in place of repeated fine-scale local solves. The numerical experiments assess both the accuracy and online efficiency of the resulting NN-MsHDG method. For two-dimensional binary permeability fields with moderate contrast, the neural method reproduces the standard MsHDG solution with relative modeling errors of a few percent while reducing the total online computational cost by factors of approximately 5 to 16, depending on the coarse trace dimension. In the high-contrast regime, the chosen polynomial coarse trace spaces already yield substantial discretization errors in the standard MsHDG method, indicating the need for more effective coarse spaces, such as coefficient-adapted spectral trace spaces. In addition, the learned local operators introduce modeling errors that become severe as the trace space is enriched. These results demonstrate the potential of neural surrogates for accelerating multiscale HDG computations while also highlighting the need for improved operator representations and a better understanding of error amplification in high-contrast problems.

We apply the adaptive Nesterov momentum (ANM) method [30] to electrical impedance tomography under the complete electrode model. The forward problem is formulated in the variational CEM setting, accounting for finite electrode size, contact impedance, insulating gaps, and the mean-free voltage gauge. The resulting nonlinear inverse problem is treated within a unified dual-to-primal framework using three classes of strongly convex structural penalties: an L2-type penalty, an L1-type penalty promoting sparse deviations from a calibrated homogeneous background, and a TV-type penalty favoring approximately piecewise-constant conductivities with sharp interfaces. The TV class is implemented using both smoothed TV and Huber-TV formulations. The data- misfit gradient is computed through CEM adjoint equations and stabilized by Sobolev smoothing. The method is evaluated on the publicly available KIT4 tank measurement data after calibration of the background conductivity and contact impedance from no-object measurements. The experiments include single, multiple, mixed-conductivity, and geometrically challenging phantom configurations. The L2-type penalty generally produces smooth but diffuse reconstructions, whereas the L1-type penalty yields cleaner backgrounds with occasional geometric distortion. The smoothed TV and Huber-TV penalties provide more spatially coherent localization and exhibit similar reconstruction behavior across most tested configurations. These results demonstrate the practical applicability of the adaptive Nesterov framework to measured CEM-EIT data.

We study a fast-forwarded quantum algorithm for solving weakly nonlinear dissipative ordinary differential equations. Our approach is a combination of the Carleman embedding technique and the linear combination of Hamiltonian simulation algorithm for linearized systems with fast-forwarded scaling. The complexity of our algorithm does not explicitly depend on the evolution time $T$, thus greatly improving the previous state-of-the-art $\widetilde{\mathcal{O}}(\sqrt{T})$ to $\mathcal{O}(1)$, and any remaining time dependence enters through the output norm and forcing parameters. We rigorously analyze the performance of this approach by convergence guarantees of the Carleman embedding for time-dependent coefficient matrices and detailed complexity estimates, and improve the realization of the Carleman-embedding-based algorithms by simplifying the post-selection step. In addition, we perform a numerical study on differential equations beyond the weakly nonlinear case, and identify possibility of achieving fast-forwarding scaling for systems with stronger nonlinearity or linear non-resonant effect.

We study a class of nonlocal evolutionary partial differential equations with weakly singular temporal kernel and spatially variable diffusion coefficient. The model is posed on $\Omega \subset \mathbb{R}$, and involves a left-sided Riemann–Liouville fractional derivative in space multiplied by a variable coefficient $a(x)$. The temporal derivative is approximated by an $L1$ type scheme, while the spatial operator is discretized by finite difference techniques, resulting in large scale all at once linear systems with a twolevel Toeplitz like structure. We develop and analyze a block lower triangular strategy that mimics the structure of the coefficient matrix while simplifying its components for computational efficiency. The analysis is carried out at the level of matrix sequences by means of generalized locally Toeplitz (GLT) theory. Within this framework, we characterize the asymptotic spectral distribution of the discretized operators and use the associated GLT symbol to guide the construction of the structured approximation. Numerical experiments using the GMRES solver demonstrate that the proposed preconditioning strategy significantly improves convergence rates, robustness, and scalability for large-scale problems. Open problems and possible extensions are briefly discussed at the end of the present work.

We propose a rigorous computational method that certifies the separator required by Temple's inequality and its Lehmann–Goerisch refinements, and with it unconditional ground-state bounds for $H=-\Delta+V$ on $R^3$ with attractive Coulomb singularities. For H$_2^+$ the method certifies $\lambda_2^{\rm eee}\ge-0.2035$ in place of the analytic $-\tfrac14$, and Temple's inequality with a global trial then gives an enclosure of width below $6.5\times10^{-51}$. For linear H$_3^{2+}$, whose central nucleus defeats the analytic separator, the certified widths are $1.3\times10^{-12}$ and $3.1\times10^{-13}$ at two nuclear separations. Such a separator, a lower bound on the second eigenvalue, is supplied by theory only for special configurations; elsewhere it is estimated and the bound is conditional. Here every hypothesis is verified in interval arithmetic at run time. Three constructions make this possible on an unbounded domain with a potential of unbounded $L^\infty$ norm: a two-sided Neumann/Dirichlet truncation with a computable remainder; a localized Hardy inequality and sharp form bound with explicit constants, restoring coercivity; and a projection stage that certifies lower bounds on the first two eigenvalues, the second serving as separator. On a box, a Lehmann–Goerisch step places $\lambda_1\in[-0.553291,-0.550740]$ for H$_2^+$. The coercivity constant, we show, must bound its own optimum from above. All drivers and certificates are released.

We study multiple orthogonal polynomials with orthogonality defined by multiple positive discrete measures on the real line. Focusing on type~I and type~II multiple orthogonal polynomials and their step-line recurrence relations, we consider the reconstruction of the associated banded upper Hessenberg recurrence matrix from given nodes and weights of the measures. We propose efficient updating and downdating procedures that modify an existing recurrence matrix when nodes are added to or removed from the discrete measures. The updating strategy is based on solving an inverse eigenvalue problem, while the downdating procedure relies on an eigenvalue deflation technique inspired by a QR-type of algorithm. Numerical experiments confirm the stability and efficiency of the proposed approaches.

In this work, we investigate the inverse problem of recovering the conductivity coefficient in an elliptic equation from noisy measurements collected at finitely many deterministic scattered points in the domain $\Omega$, and corrupted by random noise. Inspired by the regularity analysis, we propose a numerical scheme based on the regularized least-squares formulation with a $W^{1,4}(\Omega)$ penalty, and discretize the regularized problem using the Galerkin finite element method with continuous piecewise linear elements. Under suitable assumptions on the problem data, we provide an error analysis of the regularized solution and its Galerkin approximation. We establish $L^2(\Omega)$ error bounds in a high-probability sense, which depend explicitly on the regularization parameter $\gamma$, the number $n$ of data points and the mesh size $h$. We also present numerical experiments to illustrate the theoretical findings.

In the current work, we study univariate polynomial weighted histopolation on $[-1,1]$, where the data are weighted integrals over a family of intervals. After choosing a polynomial basis, the weighted moment conditions lead to a histopolation matrix whose structure depends on the weight and on the geometry of the cells. We investigate its nonsingularity, which guarantees unisolvence, together with exact diagonality of its Gram matrix, which allows its singular values and spectral condition number to be determined explicitly. For families of intervals whose endpoints belong to a fixed grid, we characterize unisolvence in terms of the connectedness of the associated endpoint graph. In the unisolvent case, this graph is a tree, and the unique paths joining consecutive grid points provide an explicit expression for the inverse matrix. This identity gives explicit formulas for the singular values of both matrices, and shows that their condition numbers in the two-norm coincide and grow linearly with the matrix size. Moreover, it yields the limiting singular value distributions of the two matrix sequences. We also establish a general diagonalization criterion based on discrete weighted orthogonality. The criterion recovers the first kind Chebyshev construction and leads to a diagonal configuration for the constant weight based on discrete sine orthogonality. For interval families with a connected endpoint graph, the corresponding moment vectors define an inner product on the polynomial space and lead to a monic basis with a diagonal weighted Gram matrix. Finally, we derive reduction formulas for cell moments associated with generalized Jacobi weights and introduce an alternative basis for shifted Jacobi weights. Applied to the Chebyshev weight of the fourth kind, this basis, together with a correction of one nonconstant element, yields an exactly diagonal Gram matrix.

We propose and analyze a first-order numerical scheme for a class of scalar Itô stochastic differential equations with almost surely positive solutions. We construct a new explicit numerical scheme, such that for any choice of the time-step size, the numerical solution is guaranteed to remain almost surely positive. The main result of this article is the first-order strong convergence of the proposed positivity preserving scheme, which is illustrated with numerical experiments and proved rigorously. It is also shown how to adapt the scheme and the results to Stratonovich stochastic differential equations with positive solutions.

We construct and analyze numerical schemes for systems of stochastic differential equations, which preserve almost surely a given hypercube of arbitrary dimension. We propose a new general class of explicit schemes, such that for any choice of the time-step size the numerical solution takes values in the hypercube. We prove strong and weak convergence results for this general class of domain preserving numerical schemes, with strong order $1/2$ and weak order $1$ in general. We also construct a variant of the scheme which achieves strong order $1$ when the stochastic differential equation is driven by a one-dimensional Brownian motion. The convergence results are illustrated with numerical experiments.

We introduce a computational method that extends the range of current global ocean climate models to non-hydrostatic scales, and establishes a convection resolving ocean model. The main obstruction towards this goal at global scale -the calculation of the non-hydrostatic pressure- is resolved by computing the pressure locally in a manner that capitalizes on particular conditions of ocean dynamics, and the structure of ocean models. Our claims are substantiated 1) by a theoretical analysis, which shows that the proposed method is compatible with ocean physics, 2) by numerical experiments, which evidence by comparison that it represents faithfully and to high accuracy non-hydrostatic dynamics, and 3) by a performance analysis that establishes that our non-hydrostatic method exceeds the cost of a hydrostatic model by a fixed factor of about $1.2$ in operation count and $1.3$ in runtime at any resolution. We show how to use the additional information contained in non-hydrostatic scales can be used to calculate the diffusivity free from numerical contributions and the mixing efficiency, transforming both from prescribed values by the modeller to calculated values from the dynamics.

A method is developed for structural analysis in which parts of an analysis domain are modeled by reduced-order models (ROMs) and the remaining part is modeled by a high-dimensional model (HDM). The multiple ROMs and HDM are combined based on an iterative domain decomposition method (DDM) with an arbitrary number of non-overlapping subdomains. A diagonal preconditioner is derived for the iterative DDM. The method is implemented based on a framework of the existing DDM-based analysis code. The ROM is implemented on GPU. Snapshot generation methods for the ROM are also proposed. Static structural analyses of nuclear power plant models are conducted using the combined multiple ROMs and HDM with multiple-subdomains model. The analyses are performed in parallel by the hybrid MPI-CUDA implementation. The convergence property of the iterative DDM is improved by the preconditioner.

Computational models of tumor growth form a contemporary area of study. Diffuse interface models have been proposed as an important tool for cancer modeling. In this work, we consider preconditioners for numerical simulations of two tumor growth models based on diffuse interface models. For numerical efficiency, we propose block preconditioners that rely on using effective Schur-complement approximations. We show that these methods combined with state-of-the-art adaptive finite element discretizations lead to robust simulations in two and three dimensions. In extensive numerical experiments, we show that the proposed methods show robust convergence behavior.

The Levenberg-Marquardt (LM) algorithm is the most widely used method for solving nonlinear least-squares problems, as it combines the robustness of steepest descent with the fast local convergence of the Gauss-Newton method. However, its computational cost can become prohibitive for large-scale problems because each iteration requires solving a large damped linear system, and conventional step acceptance strategies may require repeated solves as the damping parameter is adjusted. Despite this computational challenge, many large-scale least-squares problems exhibit effective low-dimensional structure, with only a small number of parameter-space directions strongly informed by the data. We propose an adaptive hybrid subspace Levenberg-Marquardt (HSLM) algorithm that constructs a low-dimensional subspace from complementary sources of gradient, memory, Krylov-subspace, and randomized curvature information and computes a spectrally damped LM step within this subspace. A distinguishing feature of the method is a deterministic adequacy monitor that quantifies how much descent information is captured by the reduced space and adaptively enriches the subspace when necessary. Step acceptance is decoupled from damping adjustment: Armijo backtracking determines the accepted step length, while the ratio of actual to predicted reduction is used solely to update the damping parameter, thereby avoiding repeated damped-system solves during step acceptance. For the HSLM algorithm, we establish global convergence to stationarity and prove local linear and superlinear convergence. Numerical experiments on neural-network training problems show that HSLM achieves convergence behavior comparable to classical and Krylov subspace LM (KSLM) while substantially reducing per-iteration computational cost, with increasing advantages observed as the parameter dimension grows.

The Stochastic Simulation Algorithm (SSA), widely considered an exact algorithm for stochastic chemical reaction networks, suffers from high computational cost. In this work, we propose a data-driven effective model that operates on a user-defined coarse time step independent of the underlying microscopic reaction-event scale. This is accomplished by directly approximating the finite-time transition kernel of the continuous-time Markov chain induced by SSA, using a generative machine learning model trained on short bursts of SSA simulation data. The trained model constructs a stochastic propagator that recursively generates statistically consistent trajectories at the constant coarse time step, with significantly reduced computational cost. In this paper, we employ conditional normalizing flow as the stochastic propagator. A comprehensive set of numerical examples is presented to demonstrate the accuracy and efficiency of the proposed method.

This paper summarises our previous work on a fully $3$D Summation-by-Parts scheme, derived for a class of linear wave equations on hyperboloidal slices on a fixed Minkowski background. The scheme is derived in spherical polar coordinates, and allows having grid points at the origin and on the $z$-axis, despite coordinate singularities, and at infinity, by introducing compactification followed by rescaling, and is proved to be stable. Reducing it to the standard Cauchy problem, or to finite spacelike slices with an outer boundary, will follow a similarly. Second-order accurate finite-difference methods are used to implement this scheme numerically, but higher-order finite-difference or spectral methods could also be used. Kreiss-Oliger dissipation operators are generalized to curvilinear coordinates and are defined everywhere in the domain, including at the boundary points, such that they satisfy the dissipative property in the energy norms. We also propose new norm convergence tests that include all the grid points at all resolutions and produce more accurate results. Promising results are obtained, giving hope for application to fully nonlinear systems, like the Einstein Field Equations, and extracting the resulting gravitational waves free of systematic errors or gauge ambiguities.

Gradient-flow models are characterized by an intrinsic energy-dissipation structure, and faithfully preserving this structure at the discrete level is important for stable and reliable long-time simulation. To this end, we develop a Conditioning-Aware Representation Enhancement with Scalar Auxiliary Variable (CARE-SAV) framework, which constructs a compact spatial approximation space from flexible candidate features and evolves the gradient-flow dynamics directly within this space. The resulting fully discrete scheme preserves the discrete energy-dissipation law while providing a flexible alternative to conventional prescribed spatial discretizations. Rigorous analysis establishes the approximation capability, solvability, stability and convergence of the proposed method. Numerical experiments on representative gradient-flow problems demonstrate its accuracy, robustness and computational efficiency. We believe that CARE-SAV could provide a simple, flexible, and computationally efficient paradigm for structure-preserving discretization of gradient-flow problems.

The OBABO and BAOAB schemes and the other standard Strang splittings of kinetic (underdamped) Langevin dynamics are widely used Markov chain Monte Carlo algorithms. Under a suitable friction scaling, the underlying diffusion relaxes on a ballistic time scale, suggesting that these discretizations, suitably tuned, sample targets with condition number $\kappa$ in $O(\sqrt{\kappa})$ iterations. We prove that no fixed choice of step size and friction, based only on the curvature bounds and the dimension, achieves this acceleration: total variation mixing time lower bounds for OBABO show that ballistic cold-start mixing fails uniformly over the smooth strongly convex class, and the lower bounds extend, with the same orders, to BAOAB and the other four Strang splittings. The proof transfers non-acceleration from optimization to sampling. Eliminating velocity gives an exact noisy heavy-ball recursion, and by the non-acceleration theorem of Goujaud, Taylor and Dieuleveut, for every tuning either some Gaussian target has a mode with relaxation time at least of order $\kappa$, or an attracting cycle exists on a smooth potential; dilating such a potential as $U_R(x)=R^2U(x/R)$ preserves its curvature bounds and produces metastability for a number of steps exponential in $R^2$, from an initial state at Wasserstein distance $O(R)$ from equilibrium. Using contraction estimates of Leimkuhler, Paulin and Whalley and a Wasserstein-to-total-variation regularization estimate, we prove a complementary upper bound of $O(\kappa)$ steps, up to logarithmic factors, for a fixed-parameter OBABO tuning. Hence, among fixed-parameter OBABO tunings, the optimal condition-number dependence of cold-start total variation mixing over this class is linear, up to logarithmic factors. A direct Gaussian calculation also rules out fixed-parameter acceleration for the left-endpoint exponential integrator.

This paper introduces and analyzes a novel algorithm—Truncated Fourier Filtering (TFF)—for the fast, high-order accurate evaluation of standard integrals and convolutions involving piecewise-smooth (possibly discontinuous) integrands over general $m$-dimensional domains ($m \ge 1$) employing an $m$-dimensional Cartesian grid. For an $N$-point discretization, the method runs at a computational cost of $\mathcal{O}(N)$ operations for standard integrals and $\mathcal{O}(N \log N)$ operations for convolutions, following, in either case, a one-time $\mathcal{O}(N \log N)$ precomputation step (not required in dimension $m = 1$). The core idea underlying TFF is to approximate the characteristic function of the integration domain by a truncated Fourier expansion of it over a suitably extended periodic domain, and to evaluate the resulting integrals via trapezoidal quadrature on a Cartesian grid with an appropriately chosen discretization size. Despite its conceptual simplicity, TFF attains high-order accuracy even for complex, possibly non-smooth or even non-Lipschitz geometries. A complete theoretical analysis is provided that establishes the superalgebraic convergence (i.e., convergence faster than any negative power of $N$) of the overall approach.

We construct and analyze structure-preserving numerical schemes for a two-stage local–nonlocal dispersal system. On general bounded connected smooth habitats, we prove fixed-$\delta$ spatial consistency. We introduce a semi-implicit scheme that is uniquely solvable for every $\Delta t>0$, preserves nonnegativity, and exactly retains the sign of the semidiscrete spectral threshold. A negative threshold yields geometric extinction of the numerical solution, while a positive threshold makes the zero state of the numerical scheme linearly unstable. On rectangular habitats, we construct a reflected Cartesian discretization that is asymptotically compatible with the Neumann local limit, with convergence estimates uniform with respect to the ratio between mesh size and interaction scale. Numerical experiments illustrate the structural, threshold, and local-limit behavior.

Simulating large-scale free-surface water by coupling a localized 3D fluid solver to a cheaper 2D surface model has long faced a mismatch in wave dynamics: efficient 2D wave models used in graphics are typically either linear or non-dispersive. These models are fast, simple, and accurate for calm, small-amplitude seas, but coupling them with strongly nonlinear 3D solvers produces visible reflections and artifacts at the 2D–3D interface. We address this problem by introducing a nonlinear and dispersive 2D wave model based on the canonical Zakharov formulation. Its Hamiltonian structure, in which the surface elevation and surface potential form a canonical pair ($\eta$, $\psi$) governed by the wave energy, enables a canonically consistent two-way coupling scheme, allowing information to pass smoothly across the 2D–3D interface. Our 2D solver reduces mean wave-height error by 1.7–5$\times$ over SWE, BEM, and Airy baselines while running more than $10^3\times$ faster than BEM; it achieves greater nonlinear accuracy and coupling fidelity than SWE and Airy, with minor losses in speed and stability. Coupling it with a 3D Navier–Stokes solver yields a full system that suppresses visible seam artifacts across a range of experiments, including dispersion-matching and Kelvin-wake tests, and runs over 4$\times$ faster than a pure GPU NB-FLIP simulation on the same domain.

Waveform inversion seeks to estimate the wave speed of a heterogeneous, inaccessible medium, from time-resolved measurements of the waves at user controlled sensors. We consider this inverse problem for acoustic waves and an active array of source/receiver sensors that emit probing signals and measure the generated pressure waves. The forward map, from the wave speed to the measurements, is nonlinear and oscillatory. The oscillations cause cycle skipping, the main impediment to using the standard, nonlinear least-squares data fitting formulation, known as full waveform inversion (FWI). A recently introduced alternative waveform inversion approach computes from the measurements an algebraic surrogate of the wave operator, a reduced order model (ROM) matrix, which is then used to estimate the wave speed. The mapping from the measurements to the ROM is nonlinear, but well understood. It is computed efficiently, in a non-iterative manner. The nonlinear mapping from the ROM to the wave speed is less understood, and its approximation involves time-consuming optimization. Our goal in this paper is to use a neural network to map the ROM matrix to a nearby one, that has a simpler and explicit dependence on the wave speed. This simplifies and reduces the computational cost of the ROM-based waveform inversion. We introduce the methodology, called ROMNet, and test it with numerical simulations, using two training data sets: The first set consists of random media with variations of the wave speed modeled by a superposition of Gaussians with random amplitudes and standard deviations. The second is the publicly available GeoFWI dataset introduced for benchmarking FWI using deep learning. We compare the performance of ROMNet with the direct ROM-based inversion and with two representative deep learning approaches to FWI: “Fourier-DeepONet" and “InversionNet".