Papers

We propose and analyze the random batch vortex blob method for the 2D Navier–Stokes equation in vorticity form on the whole plane. The vortex blob method is based on an interacting particle system of $N$ particles with computational complexity of $O(N^2)$, which is reduced to $O(N)$ by the random batch method \cite{JinLiLiu2020}. Our main result is a quantitative law level mean field error estimate whose dependence on the blob radius remains algebraic. We treat the two main error mechanisms separately. The random batch error is controlled through a locally coupled auxiliary partition, a symmetric law comparison, and Fisher-information dissipation. The mean field fluctuation is estimated by exploiting the oddness and divergence-free structure of the Biot–Savart kernel. For smooth, strictly positive initial vorticity, we prove on every finite time interval a normalized relative-entropy bound of order $ O\!\left(\varepsilon^{-4}\tau^2+N^{-1} \right), $ with constants independent of $N$, $\tau$, and $varepsilon$. Here $\tau$ is the batch refreshing interval and $\varepsilon$ is the blob radius. As a consequence, the fixed-particle marginals converge strongly in $L^1$ to tensor products of the regularized vorticity solution when $N\to\infty$ and $\varepsilon^{-2}\tau\to 0$.

We give a negative solution to the inverse generator problem on Hilbert spaces. More precisely, we construct a bounded operator $A$ with dense range on a Hilbert space $H$ that generates a bounded, strongly stable $C_0$-semigroup, while $A^{-1}$ does not generate a $C_0$-semigroup. We also construct an exponentially stable generator $A$ with $0 \in \rho(A)$ such that the inverse semigroup is unbounded and grows at least double logarithmically. For the latter generator, every Cayley transform satisfies the ordinary Kreiss resolvent condition but is neither strongly Kreiss bounded nor power bounded. Moreover, its powers satisfy a doubly logarithmic lower bound. Therefore, the Crank–Nicolson scheme is unstable in operator norm both for every fixed step size over long times and under mesh refinement at any fixed final time. Our counterexamples are deduced from a common finite-dimensional construction. For $\alpha\in(0,1)$, we use explicit bases of $\mathbb C^{2n}$ whose partial-sum projections are uniformly bounded and whose unconditionality constants are comparable to $n^\alpha$. The matrices underlying the counterexamples are then obtained as Schauder multipliers with respect to these bases, using a sequence of eigenvalues whose moduli decay doubly exponentially.

We consider the problem of minimising an energy over an unknown manifold that is given implicitly by a point cloud. For a known manifold one can define a gradient descent scheme analogously to the classical construction in Euclidean spaces. However, when the manifold is not known one has to simultaneously estimate the manifold whilst minimising the energy. We define a gradient descent scheme which remains in a neighbourhood of the manifold and, under suitable stability and sampling assumptions, converges to a neighbourhood of a local minimiser whose size vanishes as the time step and sampling errors vanish. As an example we show the application of the methodology to learning operator corrections in inverse problems.

We study Muon, a recently proposed matrix-aware optimization method, in the context of the Stiefel manifold. This manifold consists of matrices with orthonormal columns and is ubiquitous in machine learning and scientific computing. Existing extensions of Muon to this manifold rely on heuristic, approximate, or iterative updates with varying computational efficiency. We show that the corresponding Stiefel Muon update admits an exact closed-form solution and use this result to develop Skewon, a practical algorithm for orthogonality-constrained optimization with an efficient implementation. We further establish first-order convergence guarantees for Skewon in the smooth non-convex setting.

Conditional expectation operators (CEOs) and their associated conditional mean embeddings (CMEs) play a central role across applied mathematics and machine learning, appearing in nonparametric regression, Bayesian inverse problems, and Koopman operator theory. A fundamental question is when a CEO maps a function space on $\mathcal{Y}$ into a prescribed function space on $\mathcal{X}$, particularly a reproducing kernel Hilbert space (RKHS). We show that such mapping properties are characterized by the regularity of the Radon–Nikodym density of the conditional law, and establish a simple, verifiable sufficient condition under which the CEO is bounded and Hilbert–Schmidt. For RKHSs norm-equivalent to Sobolev spaces, this condition reduces to Sobolev regularity of the conditional density. The result yields a direct route to validate CME representations and error bounds for Galerkin-type and CME-based estimators. We verify the regularity condition in three settings: nonparametric regression, Bayesian inverse problems, and Koopman operator theory for stochastic dynamical systems. We show in each case that classical regularity results on the underlying probabilistic model imply the required mapping properties. The resulting framework offers a unified perspective on conditional expectation operators across probability, operator theory, kernel methods, and stochastic dynamics.

Mechanistic chromatography models are valuable for process development, but gradient-based parameter estimation and optimization can be hindered by computational cost and the effort of deriving objective-function gradients. To address this concern, a fully differentiable Julia chromatography solver, ChromOps.jl, is presented that combines high-order spatial discretization with discrete forward sensitivity analysis (DFSA). Two high-order spatial discretizations, finite difference summation-by-parts (FD-SBP) and the discontinuous Galerkin spectral element method (DG-SEM), are compared on a 6-component ion-exchange chromatography problem with steric mass action kinetics. Both attain their theoretical convergence rates on manufactured problems, while FD-SBP shows favourable work-precision performance for both forward (primal) and dual-valued simulations and is conceptually simpler to implement. DFSA scales as tDFSA approx 1.4(1+Np)tfwd for up to Np = 24 parameters propagated in a single ForwardDiff.jl chunk, i.e. about 1.4 forward solves per additional parameter. Because the primal and dual-valued solves share the same code path, DFSA provides accurate gradients of user-defined objectives without any hand-derived chain-rule code, enabling accessible gradient-based parameter estimation and optimization.

Machine learning offers a promising avenue to accelerate physical simulations by replacing computationally expensive traditional Partial Differential Equation (PDE) solvers with fast, differentiable surrogate models. However, standard auto-regressive ML emulators often suffer from error accumulation over long horizons and struggle to capture the stochasticity of complex physical systems. In this paper, we propose Kastor, a comprehensive methodology to adapt a deterministic physics foundation model into a highly efficient and accurate generative surrogate. First, we introduce a two-stage inference scheme that combines a large-stride causal auto-regressive model with a non-causal temporal super-resolution network, significantly reducing error accumulation while minimizing computational cost. Second, we present Mean prediction regularization (MPR), a novel training objective that constrains the generative model to predict the deterministic distribution mean under null noise conditioning. This regularization dramatically improves the performance and stability of both Functional Generative Networks (FGN) and diffusion-based emulators. Finally, we demonstrate that incorporating spatial gradient matching improves the accuracy and physical fidelity of the simulations as measured by power spectrum density. Extensive evaluations on diverse simulation datasets of the benchmark The Well show that with these components, our model outperforms competing methods in forecasting accuracy, spectral consistency, and computational efficiency. Our model achieves a 42.9% average reduction in forecasting compared to our reference based on the Walrus finetuning methodology, and outperforms Walrus for 8 out of 10 datasets on variance-normalized RMSE (VRMSE).

This work develops a thermodynamically consistent phase-field model for tumor growth based on the energetic variational framework. The model couples the Cahn-Hilliard equations for tumor evolution and nutrient transport with the incompressible Navier-Stokes equations. A first-order time discretization scheme based on the Multiple Scalar Auxiliary Variables (MSAV) approach together with a pressure-correction strategy is proposed to efficiently handle the nonlinear and coupled structure of the system. The proposed scheme is rigorously proved to be unconditionally energy stable and mass conservative. Furthermore, optimal first-order temporal error estimates are established for the tumor phase-field variable, the nutrient concentration, and the fluid velocity. Finally, numerical experiments demonstrate the effectiveness and robustness of the proposed method and verify the theoretical convergence rates.

Reduced Basis (RB) methods provide low-dimensional approximations of parametrized partial differential equations with controllable accuracy. We discuss the construction of RB approximations for transient thermal simulation of an induction motor and compare their performance to Finite Element (FE) models and Lumped Parameter Thermal Networks (LPTNs) after calibration to measurements. Numerical results demonstrate that RB models can achieve FE accuracy, allow for a fully automatic construction, and offer computational efficiency comparable to traditional LPTNs.

Physics-informed neural networks (PINNs) often fail to accurately resolve partial differential equations (PDEs) with high-frequency or multi-scale solutions, as well as strongly nonlinear problems. Two factors underlie this difficulty: spectral bias, the tendency of neural networks to underfit high-frequency features; and representation-coefficient coupling, the entanglement of representation learning and coefficient fitting within a single nonconvex optimization objective. In this work, we propose the Fourier-enhanced alternating Levenberg–Marquardt PINN (FALM-PINN), an optimization framework that decouples representation learning from coefficient fitting. The upper-level problem learns a Fourier-enhanced basis that enriches the latent space with high-frequency components, while the lower-level problem resolves the coupling by fitting the projection coefficients on this basis, solving a nonlinear least-squares problem with the Levenberg–Marquardt algorithm. The framework applies to general nonlinear and coupled PDE systems, and reduces to a single-step convex optimization problem for linear PDEs. We prove global convergence of the alternating training scheme in both cases. Numerical examples on multiple challenging high-frequency and nonlinear PDEs show that FALM-PINN achieves relative $L^2$ errors up to two orders of magnitude lower than state-of-the-art baselines.

We study a fully implicit Rothe time discretization for a cutoff first-order formulation of the Westervelt equation. The key ingredients are the enthalpy variable and the primitive mobility variable, which turn each nonlinear time step into a uniformly monotone elliptic problem and avoid higher-order energy estimates and inverse inequalities. For every time step, the discrete problem reduces to a monotone elliptic equation, which yields well-posedness of the Rothe scheme by the Browder-Minty theorem. We derive a discrete energy inequality, establish compactness for the transformed variable by an Alt-Luckhaus type argument, and pass to the limit to obtain existence of weak solutions for the cutoff first-order system. We prove a weak-strong uniqueness principle for the cutoff problem and formulate a conditional a posteriori criterion under which the cutoff is inactive. For sufficiently regular solutions of the forced cutoff problem, we establish consistency estimates which identify the temporal residuals produced by the Rothe approximation in the natural discrete-in-time weak norms. These estimates provide a rigorous consistency basis for the observed temporal behavior. Finally, numerical experiments illustrate the stability of the scheme, its observed near first-order temporal behavior, and inactivity of the cutoff along the computed discrete trajectories in the tested regimes.

Deep neural networks have shown great empirical success in the solution of a wide variety of ill-posed inverse problems in imaging. Yet, very few works have studied their behavior in the limit that turns the discretized ill-conditioned problems into truly ill-posed ones, i.e., for an increasing resolution of the discretization. In this work, we review common approaches to neural operator learning in architectures that resemble a U-Net, one of the most common classical architectures for inverse imaging problems. We discuss advantages and drawbacks of the respective approaches, consider a 1D toy example for improved interpretability, and present extensive numerical experiments on how different types of neural operator U-Nets can improve a first (crude) limited angle CT-reconstruction. In particular, we study how well networks trained for a certain resolution of the discretization generalize to other resolutions. Our finding is that while U-shaped neural operator architectures are by design resolution-invariant, the classical U-Net architecture seems to be more robust with respect to resolution changes than expected.

We study the identification of effective coefficients for wave equations in heterogeneous media. Such equations arise in the modeling of spatio-temporal metamaterials, where the underlying material properties exhibit variations in both space and time. While homogenization provides effective models in the asymptotic regime of vanishing microscopic scales, determining macroscopic parameters from observations of wave propagation remains challenging. We introduce an optimization-based method that identifies a constant effective coefficient by minimizing a cost functional. The approach is designed for situations in which the underlying space-time-dependent coefficient is unknown, while the solution is available in space and time by measurements. We extend optimization-based coefficient identification techniques from elliptic multiscale problems to wave equations and prove that, provided a homogenized limit exists, the identified coefficient converges to the homogenized coefficient as the microscopic scale tends to zero. Furthermore, we establish convergence of the corresponding effective solution towards the heterogeneous solution, including strong convergence in $L^2$ and weak convergence of first-order space and time derivatives. Numerical experiments further demonstrate the performance of the method for non-periodic space-time heterogeneous media, including cases for which a homogenized limit is not known to exist.

We derive an explicit one-step kernel formula for the kinetic Fokker-Planck equation. The construction begins with a reference dynamics admitting an explicit transition kernel and uses a carefully chosen auxiliary weight to encode the drift perturbation. The resulting kernel formula provides an explicit approximation of the density evolution. In a weighted phase-space setting, we establish local weak consistency and first-order finite-time weak convergence under suitable regularity and moment assumptions. We also give a formal variational interpretation of the formula in terms of a regularized Wasserstein-type proximal operator and present analogous constructions for a broader class of parabolic equations. Numerical experiments illustrate the accuracy of the kernel approximation and its application to deterministic particle schemes for kinetic sampling.

In this work, we present two novel strategies to impose high-order boundary conditions on fixed and moving curved domains, approximated with piecewise affine triangulations. Achieving high-order accuracy on curved domains requires tackling both the PDE discretization error and the geometrical error simultaneously. While the former can be reduced by employing high-order numerical methods such as finite volume and discontinuous Galerkin, the latter demands either a high-order parametrization of the physical domain or a consistent approximation of the boundary conditions. Minimization-based approaches like the Reconstruction for Off-site Data (ROD) method allow one to skip the construction of high-order curvilinear meshes by defining high-order consistent boundary conditions on a computational boundary that does not match the physical one. The ROD approach mitigates the second-order geometrical error by retrieving a modified polynomial in each boundary cell, which enforces the boundary conditions exactly on the physical boundary. However, the standard ROD method requires the inversion of a local linear system, whose cost grows with the polynomial degree and mesh refinement. Inspired by a recent one-dimensional analysis, we show that ROD-type approaches can be recast as simple polynomial corrections, applicable without any linear system inversion. This greatly simplifies the development of minimization-based boundary treatments and reduces the associated computational cost. To prove the wide applicability of our strategy, we develop it within a Runge-Kutta discontinuous Galerkin framework and an ADER arbitrary-Lagrangian-Eulerian finite volume framework for compressible flows on fixed and moving domains. Several numerical experiments with Dirichlet and slip-wall boundary conditions are presented, with convergence analysis up to fifth order in both 2D and 3D.

We introduce Potential Matching Optimal Transport (PMOT), a potential-flow framework for general $p$-cost optimal transport with $c_p(x,y)=\|x-y\|^p$. PMOT parameterizes the CNF velocity field with a scalar potential in the generalized Benamou–Brenier form for the chosen exponent $p$. It trains the potential gradient with a self-induced matching loss along straight bridges determined by the model's own endpoints, while allowing flexible terminal distribution matching. Our main result establishes zero-loss exactness: under the stated regularity, exact terminal matching, and uniqueness assumptions, any zero-loss solution satisfies the generalized Benamou–Brenier optimality system and recovers the corresponding $p$-optimal transport map and dynamics. On synthetic benchmarks, PMOT learns $p$-specific maps that agree with the corresponding $p$-matched OT references. It also remains competitive as a likelihood-based density model on high-dimensional tabular data, and an MMD-based color transformation experiment demonstrates flexible sample-based terminal matching.

The higher-order Craig-Bampton (HCB) method improves reduced-model accuracy through successive residual-mode enrichments. This study proposes an a priori estimator for the eigenvalue error remaining in the first-order HCB model (HCB-1) without using the full-order eigensolution. The estimator is defined as the relative difference between the HCB-1 and HCB-2 eigenvalues. Under symmetric Ritz approximation assumptions, ordered positive eigenvalues, nested augmented spaces, and preservation of retained eigenvalues by the System Equivalent Reduction Expansion Process (SEREP), the estimator is shown to be non-negative and to provide a lower bound on the relative HCB-1 eigenvalue error. A projection-based expression and an exact error identity are also derived, clarifying the conditions under which the estimator is accurate and showing that the lower-bound property alone does not guarantee a prescribed tolerance. Numerical results for a plate, an elbow pipe, and a reactor pressure vessel show estimation-to-error ratios of about 0.85-0.99, while a cost assessment reports the additional expense associated with the HCB-2 reference solution.

This paper presents a novel higher-order multi-scale (HOMS) computational framework for efficient, high-accuracy, and low-cost simulation of nonlinear hygro-thermo-mechanical (H-T-M) coupling problems in heterogeneous structures. The inherent nonlinearity in the investigated model stems primarily from temperature- or moisture-dependent material properties, and this model also accounts for temperature-dependent internal heat source and moisture sink terms induced by exothermic, moisture-consuming chemical reactions (e.g., hydration). The main contributions of this work are as follows. First, a high-accuracy multi-scale asymptotic model incorporating higher-order correction terms is constructed for nonlinear H-T-M coupling problems in heterogeneous structures with highly spatial inhomogeneity, using the multi-scale asymptotic approach together with Taylor series expansions. Second, rigorous error estimates in both point-wise and integral senses are derived for the multi-scale asymptotic solutions, which theoretically demonstrate the necessity and superiority of the proposed HOMS method. Third, an efficient two-stage numerical algorithm with off-line and on-line stages is developed, based on finite difference and finite element methods, and its convergence is also proved rigorously. Finally, two- and three-dimensional numerical experiments are performed to assess the computational performance of the proposed HOMS approach, showing excellent numerical accuracy and robustness with low computational overhead.

Quasi-static electroporoelasticity couples Maxwell's equations with Biot's poroelasticity through electrokinetic coupling between the electric field and the pressure gradient. By introducing the total pressure, the electroporoelasticity equations are reformulated as a five-field system to address poroelastic locking in the nearly incompressible regime. For the resulting five-field system, a monolithic weak formulation is derived, together with a continuous stability estimate. A second-order backward differentiation formula (BDF2) time discretization and a mixed finite- element spatial discretization are then introduced, yieling to a fully discrete monolithic scheme. Building on this scheme, we develop an iterative decoupling method that alternates between an electromagnetic subproblem and a poroelastic subproblem, and prove its geometric convergence to the monolithic solution with an explicit mesh-independent contraction factor under the physical coupling condition. An algebraically equivalent reduced form is also presented, in which the electromagnetic block is solved only once per time step, while the poroelastic block is solved iteratively with electric-field correction updates obtained from the pressure-gradient feedback. Numerical experiments verify the theoretical predictions and demonstrate the locking-free performance.

We develop and analyze a family of divergence-free kernel interpolation methods for tangential vector fields on the unit sphere. Starting from scalar radial kernels in Euclidean space, we construct tangent-valued, surface divergence-free matrix kernels without repeatedly applying surface differential operators. The construction separates the geometric enforcement of the divergence-free constraint from the choice of scalar generator and admits lower-order variants with reduced regularity requirements compared with classical potential-based methods. Using vector spherical harmonics, we derive explicit kernel representations and characterize their Fourier multipliers. We also introduce an inverse Laplace-Beltrami construction that preserves the multipliers of the underlying scalar zonal kernel. For interpolation at scattered nodes, we establish a lower bound for the smallest eigenvalue of the interpolation matrix and derive pointwise and Sobolev error estimates, including superconvergence for targets smoother than the native space. Numerical experiments corroborate the theoretical convergence and stability results, while additional examples on nonspherical surfaces illustrate the applicability of the kernel formula beyond the sphere.

Recently, novel integrability techniques, and most notably Birkhoff coordinates followed by an explicit formula, have been introduced to solve the Benjamin–Ono (BO) equation, leading to major theoretical and computational advances. We propose to perturb this explicit BO formula via a novel splitting method for numerically solving a class of PDEs, which includes the Intermediate Long Wave (ILW) equation, as well as other quasilinear equations that are not necessarily integrable. Using the Birkhoff coordinates of BO we establish first-order convergence in the $H^s$-norm ($s \ge 0$) while requiring only one additional derivative on the initial data $u_0 \in H^{s+1}_{0}$, significantly improving upon the regularity requirements of typical splitting methods for nonlinear dispersive equations. Furthermore, we prove that in the deep-water limit, the scheme for the ILW equation converges to the BO solution. Computational advantages of these schemes are shown in simulations: unlike classical splitting methods, our approach does not require a restrictive quadratic time step condition to nearly preserve energy over long time scales, thus rendering efficient and accurate long-time simulations feasible. As an application, we numerically explore the soliton resolution conjecture for both the ILW and a KdV–BO equation.

We develop and analyze a second-order, fully implicit, monolithic time-discretization for the coupled time-dependent Stokes and quasi-static Biot system. The Biot subsystem is formulated in a three-field total-pressure formulation, in which the total pressure is introduced as an additional unknown together with the solid displacement and pore pressure. This reformulation improves robustness in nearly incompressible regimes and suppresses volumetric locking. All variables are discretized in time using the BDF2 scheme, and the coupled problem is solved monolithically without operator splitting. Discrete energy stability is established using the BDF2 $G$-stability identity. The resulting estimates exhibit robustness with respect to the Lamé parameter, Biot–Willis coefficient, and storage coefficient. A consistency–stability argument is then used to derive second-order convergence in time together with optimal-order convergence in space. Numerical experiments confirm the theoretical results and demonstrate second-order temporal accuracy in the corresponding energy norms.

The Material Point Method (MPM) has drawn great attention in the numerical modelling of large deformation, geotechnical problems. The popularity of MPM is mainly because its formulation shares significant similarities with the Finite Element Method. In MPM, the iteration points can move independently from the mesh, allowing for the resolution of large deformation problems. However, because of this, the original MPM formulation suffers from the well-known cell-crossing noise and volumetric-locking instabilities, resulting in a strongly oscillated stress field. A novel implicit locking-free B-spline MPM that controls stress oscillations to a negligible level is proposed in this paper. A novel, but very straightforward B-spline shape function implementation procedure, avoids the need for a complex material point searching algorithm, providing seamless transformation from the original MPM to this robust B-spline MPM, aiming at modelling large-strain geotechnical problems. The newly proposed volumetric locking mitigation strategy is also very easy to implement, which facilitates the reproducibility of this research. The proposed method is validated against three numerical studies: granular column collapse experiment, slope failure and footing with large penetration. The proposed numerical method agrees well with experiments reported in the literature and previous numerical studies. Also, these numerical examples show that the proposed method provides a more prominent stress field than other available methodologies.

In this work, we consider a structured tensor approximation problem, where only a limited number of lateral slices are observed. The proposed algorithm , called Basis and Manifold prior Tensor Approximation (BMTA), exploits both global and local structures of the evolution of a global tensor. Specifically, BMTA integrates two signal models: (i) a quasi-basis model that captures smooth global variations along a physical trajectory, and (ii) a manifold-guided interpolation model that characterizes local relationships among tensor slices. A low-rank Tucker reconstruction framework is incorporated to efficiently capture the priors, resulting in coefficients for basis function estimation and a tensor optimization. In addition, we provide a theoretical analysis which establishes a non-asymptotic reconstruction error bound that characterizes the effects of sampling complexity, optimization convergence, and model mismatch. Numerical experiments are performed on both synthetic and real-world datasets, including quantum chemistry and spatiotemporal sensing applications.

Supervised training of finite-element (FE) surrogate models requires reference solutions, and each reference solution is obtained by solving the system that the surrogate is intended to replace. The assembled discrete potential energy provides a training signal that requires no reference solution. This note records, with proofs, the identities that make this signal exact for linear elastostatics: the difference between the energy of a prediction and the energy of the reference solution equals one half of the squared stiffness-norm error, and the gradient of the energy equals the stiffness-weighted error. Label-free discrete-energy minimisation and supervised regression in the stiffness norm therefore have the same unique minimiser and identical gradients at every point. Around this central result, the note states a conditioning lemma that bounds the displacement error by the energy gap, a modewise contraction identity that explains why the Euclidean displacement error is an unsuitable primary metric, the Chebyshev bound that governs conjugate-gradient post-processing of surrogate predictions, and a conditional latent-separation proposition for joint-embedding predictive architecture (JEPA) pretraining on a shared stiffness operator, with an explicit numerical counterexample that delimits its scope. Every claim with numeric content is implemented as an executable falsification check; the checks were executed twice, on synthetic test problems and on a probe set of 16 instances from the validation split of a pre-registered experimental run, and every inequality holds, with the measured tightness reported. A closing section explains why the construction does not extend to elastodynamics through direct minimisation of the action functional, and which time-discrete formulation restores exactness.

This paper introduces a contour integral method (CIM) for efficiently computing outputs of parametric linear systems in control form over specified time intervals and to a user-prescribed accuracy. The CIM approximates the inverse Laplace transform via a quadrature rule applied along a modified integration contour. For parametric systems, we show how CIM integrates effectively with projection-based model order reduction (MOR) where a greedy algorithm builds the projection spaces following an error estimate we derive for this setting. We additionally demonstrate that the developed projection framework naturally enforces Hermite interpolation conditions. This combination substantially lowers the cost of evaluating the input-output relations across the parameter domain, for a wide range of input functions, and for initial conditions well captured by a low-dimensional subspace.. We demonstrate the accuracy and efficiency of the approach on benchmark non-parametric and parametric control systems, comparing against state-of-the-art projection-based MOR methods.

We construct exact transparent boundary conditions (TBCs) for a time-continuous, spatially discrete Schrödinger equation that models a one-dimensional quantum lattice. Using a recently developed exact solution for the discrete system, we derive the Dirichlet-to-Neumann maps analytically via Laplace transforms. This yields a convolution-type boundary condition governed by Bessel functions. We rigorously demonstrate the consistency of this discrete formulation with its continuous counterpart in the continuum limit. Additionally, we present an efficient time-discretization scheme based on the trapezoidal rule for practical implementation. Numerical experiments using a Crank-Nicolson solver verify that our proposed TBCs eliminate spurious backscattering entirely and preserve reflectionless propagation of a Gaussian wave packet exiting the computational domain

We develop an infinite-dimensional Bayesian framework for recovering spatially varying log-permeability from pressure observations in single-phase Darcy flow with a moving boundary, motivated by resin transfer moulding. The main contribution is an observation-wise representer formulation linking sensitivities of the parameter-to-observable map to the moving-boundary state and adjoint systems. This yields a reduced Levenberg–Marquardt method for maximum a posteriori (MAP) estimation on the Cameron–Martin space of a Gaussian prior and an explicit finite-rank covariance for the linearised MAP (LMAP) approximation. Once the representers are computed, the LM update and covariance assembly reduce to linear algebra in observation space, independently of the discretisation dimension. In 1D, the explicit moving-boundary solution is used to establish posterior well-posedness and existence of MAP estimators, and to derive closed-form Fréchet derivatives and representers. In 2D, we introduce a very weak mixed formulation coupling pressure, interface velocity and domain evolution. A formal one-sided directional shape linearisation yields the corresponding linearised state and adjoint systems, with observation-wise representers computed through independent, parallel adjoint problems. Numerical experiments show that the 1D LMAP approximation closely matches a reference MCMC posterior while requiring \(\mathcal{O}(10^1)\) rather than \(\mathcal{O}(10^5\text{--}10^6)\) forward simulations. In 2D, the method produces reconstructions and uncertainty estimates comparable to ensemble Kalman inversion with large ensembles, at substantially lower computational cost. The workflow is demonstrated across several geometries and experimental configurations without reformulating the inversion methodology.

We study exact simulation of diffusions via rejection sampling on path space using unbiased estimators of the density ratio obtained from Girsanov's theorem. When applied to the underdamped Langevin diffusion, it yields an algorithm for sampling from a strongly log-concave and log-smooth distribution with condition number $\kappa$, in dimension $d$, to accuracy $\varepsilon$ in Rényi divergence, in $\widetilde O(\kappa^{2/3} d^{1/3}\,\mathrm{polylog}(1/\varepsilon))$ queries. Under a third derivative bound, the dimension dependence improves to $d^{1/5}$. This improves substantially over the prior state-of-the-art complexity of $\widetilde O(\kappa d^{1/2}\,\mathrm{polylog}(1/\varepsilon))$ for the Metropolis-adjusted Langevin algorithm, and over the $d^{1/4}$ dimension dependence of Metropolized Hamiltonian Monte Carlo under the same third derivative bound. We also present applications to the mirror Langevin diffusion, and for obtaining Fisher information bounds in the non-log-concave case.

We consider the computation of internal solutions for a time domain plasma wave equation with an unknown potential $q$ from boundary response data. The internal solutions are computed by transforming known background snapshots using the Cholesky decomposition of the data-driven Gramian, or mass matrix. It was recently shown that in one dimension these data generated internal solutions converge in $L^2$ at order $\sqrt{\tau}$ for well chosen initial waves. Here we study the internal solution reconstruction in two dimensions, where a multiple input/multiple output (MIMO) setup and a block Gramian are needed, with the number of boundary sources increasing as the time sampling is refined. We show that the general error bound carries over to this setting: the distance between the data generated solutions and the best causal approximation from background snapshots is controlled by the best approximation mass matrix mismatch. Numerical refinement studies on a square domain measure the convergence of the data generated solutions alongside the best causal approximation from the background snapshots, with the relative errors appearing to go to zero at rate $\sqrt{\tau}$, and the absolute errors going to zero in $L^1$. We also show that the solution reconstructions remain accurate for high contrast composite media. Finally, since the mass matrix assembled from noisy response data can fail to be positive definite, regularization is needed; comparing a diagonal shift with an eigenvalue floor, we find the floor more robust, with reconstructions that remain more accurate than the unperturbed background field with noise that is up to ten percent of the root mean square data amplitude.