Papers

We develop a factorization framework for magnetic induction tomography in the eddy-current regime. The induced current density is represented by a divergence-free current vector potential, and the external excitation by magnetic dipole densities on an observation sphere. This choice reveals a natural physical factorization of the near-field operator into a data operator, an interior boundary response operator, and an adjoint counterpart. By introducing the Riesz map on the flux space, we recast this factorization in a Hilbert-space setting, which yields a positive-semidefinite operator with a coercive middle factor. Douglas' range theorem then identifies the range of the square root of this dissipative operator with the range of the data operator. For inclusions with a regular boundary, this range identity yields a uniqueness theorem for the support of the conductivity. This characterization leads to a non-iterative spectral indicator based on a regularized Picard quotient computed from the Hermitian imaginary part of the near-field data matrix. The reconstruction stage requires no forward solves. Numerical experiments confirm that the indicator localizes and separates the inclusions accurately, degrades gracefully under increasing noise, and captures non-convex support geometry without any convexity prior.

We employ uncertainty quantification methods to investigate how parameter identifiability changes in the vicinity of a bifurcation point. We perform numerical experiments on the non-spatial Klausmeier vegetation model with random coefficients, which describes biomass-water interactions and exhibits a fold bifurcation. We partition the domain around the bifurcation into regions with distinct convergence behaviors. In each region, we follow a UQ workflow that includes sensitivity analysis, Bayesian inference, and Fisher information evaluation to assess parameter identifiability. The results show that proximity to the bifurcation point is decisive for parameter estimation. Reliable joint inference of both model parameters is not possible when the system exhibits bistability. Furthermore, we can reconstruct the model's bifurcation pattern by Fisher information heatmaps, presenting a practical tool for bifurcation detection. Lastly, we highlight the importance of transient data for successful parameter estimation.

In this paper, a lattice Boltzmann method with an adaptive relaxation parameter and a Lax–Friedrichs-type equilibrium is proposed for hyperbolic systems with source terms. The equilibrium distribution recovers the conservative variables and physical fluxes while incorporating dissipation determined by characteristic-speed bounds. To balance accuracy and robustness, the relaxation parameter is selected from a local smoothness indicator based on characteristic projections. In smooth regions, the parameter approaches the low-dissipation limit, retaining second-order accuracy; near discontinuities, it is automatically reduced to introduce localized dissipation and suppress nonphysical oscillations. The stabilization acts directly through the local collision step and preserves the standard collide-and-stream structure, without a posteriori recomputation or interface-based limiting. Maxwell iteration establishes second-order consistency in smooth regions, and a weighted $L^2$-stability estimate is proved for linear hyperbolic systems with periodic boundary conditions under a standard CFL condition. Numerical experiments for scalar advection, Euler, shallow-water, and reactive Euler equations demonstrate the expected accuracy for smooth solutions and robust resolution of challenging one- and two-dimensional discontinuous problems, including wet–dry fronts and reactive discontinuities. A large-scale simulation of a 2D cellular detonation further demonstrates the capability of the method to resolve long-time multidimensional shock–reaction interactions.

MLegS is an open-source Modern Fortran code package for solving linear and nonlinear time-dependent partial differential equations on radially unbounded domains. It builds on an earlier in-house prototype, recasting it as a modular solver with a unified interface for fields, transforms, and spectral operators. Within this framework, the numerical formulation maps $0\leq r<\infty$ onto a finite interval, where associated Legendre functions provide the radial basis, while Fourier series represent the periodic azimuthal and axial directions. The mapped radial basis represents radial infinity without imposing an artificial boundary condition at finite radius while incorporating the regularity required at $r=0$ according to azimuthal order. The implementation distributes the corresponding field layouts using the Message-Passing Interface (MPI), with strong- and weak-scaling measurements extending to 512 cores across four compute nodes. Beyond parallel performance, three examples assess the solver across diffusion, reaction-diffusion, and vortex dynamics: closed-form diffusion solutions establish spectral spatial convergence and the design orders of both semi-implicit integrators; a Fisher-Kolmogorov-Petrovsky-Piskunov front matches the linearized solution and propagates at the classical speed within its logarithmic correction; and vortex-pair simulations based on toroidal-poloidal decomposition yield a measured Crow growth rate $2.7\%$ above the theoretical prediction and demonstrate the reconnection of the vortex pair.

Adaptive cross approximation (ACA) constructs low-rank CUR approximations from selected rows and columns of a matrix, making it attractive when individual entries are inexpensive to query but forming or repeatedly multiplying by the full matrix is not. Its practical performance, however, is not well explained by existing worst-case bounds, which can be exponentially large and often fail to predict when greedy largest-entry pivoting performs poorly. By viewing the matrix as a cross-gram matrix, we reinterpret the classical determinantal formula for the CUR residual through the lens of exterior algebra. In this representation, each residual entry is a weighted inner product between two $(k+1)$-blades divided by the normalized volume of the pivot block. The formula separates three effects: the weighted magnitudes of singular-vector rows, their angular geometry, and the conditioning of the selected pivots. From this we recover a classical $\sigma_{k+1}$-type estimate, characterize when one pivot annihilates or nearly annihilates additional rows and columns, obtain a quasi-optimality estimate for positive-semidefinite matrices, and bound the rank-one update through an explicit angle-based formula. We use this framework to explain the failure of greedy ACA on an asymmetric spiral-kernel example and to motivate a geometry-aware weighted-mass pivoting rule for radial-basis Galerkin stiffness matrices. In the reported experiments, weighted-mass pivoting often reduces the Frobenius residual relative to greedy pivoting and randomly pivoted Cholesky while retaining comparable runtime on parallel hardware.

We propose a variational mollification approach to circular deconvolution based on the reconstruction of a mollified target object rather than the exact solution itself. The method is formulated as a convex variational problem whose solution admits an explicit Fourier representation. We establish the consistency of the proposed reconstruction as the target resolution increases and derive convergence rates for deterministic data perturbations. For ordinary smooth kernels, the method achieves the classical order-optimal algebraic convergence rates under Besov-Nikolskii smoothness assumptions, while for supersmooth kernels it attains the corresponding order-optimal logarithmic rates. Numerical experiments on synthetic and wind-direction data illustrate the effectiveness of the proposed approach and confirm the theoretical predictions.

We study wetting on a doubly periodic rough surface in three dimensions. Although the liquid-vapor interface meets the solid at the local Young's angle, microscale roughness can cause the macroscopic apparent angle to differ substantially from this value. We formulate the directional contact angle hysteresis (CAH) interval in terms of apparent angles associated with pinned microscopic configurations. To approximate its receding and advancing endpoints, we evolve capillary mean curvature flow (CMCF) toward extremal stationary states. Computing these states is difficult because the pinning that produces hysteresis is generated at the scale of the roughness, whereas the apparent angle is only meaningful at the macroscopic scale, so a single uniform grid must resolve both. We therefore introduce a two-scale alternating (TSA) method based on a Schwarz decomposition: a Merriman-Bence-Osher (MBO) diffusion-generated scheme resolves the contact-line near region, while a linearized minimal-surface problem updates the far region. For an idealized reference iteration, we prove decay of an approximate interfacial energy. Numerical experiments on a representative doubly periodic surface show a strongly anisotropic CAH interval whose width varies by more than a factor of three with contact-line orientation and changes sharply near the diagonal directions. Because a stationary droplet must meet the solid at an apparent angle inside this interval, the computed anisotropy constrains which macroscopic wetted regions the surface can support; it is consistent with a square-like stationary droplet whose sides align with the diagonal directions.

Measures on function spaces arise throughout Bayesian inverse problems and generative modeling, often with low-dimensional structure relative to a tractable reference measure. We introduce the class $\mathcal{P}_\psi(\mu)$ of measures that differ from a reference measure $\mu$ only through a finite-dimensional map $\psi$ while preserving the reference conditionals on its fibers. Class members are determined by their $d$-dimensional pushforwards under $\psi$ and admit convenient block-triangular transport map representations. Draws are taken from this $d$-dimensional distribution and then completed to function space through sampling of the reference conditionals. For Gaussian references, these transport maps are finite rank perturbations of the identity. In contrast, optimal transport maps do not preserve this low-dimensional structure. For covariance perturbations that are trace class in the Cameron-Martin geometry of the reference, we show the optimal map is a trace class perturbation of the identity. Even finite rank perturbations yield corrections whose rank, governed by an invariant subspace, is typically infinite. We develop approximation theory for $\mathcal{P}_\psi(\mu)$ and error analysis for fitting within it, splitting the total error into an irreducible class error and a finite-dimensional marginal term set by the dimension of $\psi$ rather than the ambient discretization. We present numerical results including inference from low-dimensional nonlinear observation maps, deconvolution under a jump process prior, and state estimation for Navier-Stokes flows.

Many scientific problems require accurate modeling of complex physical phenomena, such as fluid dynamics or climate modeling. These phenomena often result in high-dimensional and thus computationally expensive computational models that can limit their application to real-time and multi-query problems. Model-order reduction (MOR) is an approach that seeks to approximate full-order models (FOMs) using reduced-order models (ROMs), trading a minor reduction in accuracy for a major reduction in computational cost. In this work, we propose a non-intrusive MOR method using generative machine learning by means of denoising diffusion probabilistic models (DDPMs) for generating solutions of the dynamical systems under different instances of their parameters. Unlike conventional DDPMs, which often operate purely in the spatial domain, we aim to generate spatio-temporal solutions to improve quality and temporal coherence. In addition, we embed the DDPM in a latent space obtained by sequentially applying proper orthogonal decomposition and an autoencoder to reduce the data dimensionality and the computational cost of the DDPM. We test our approach on a parametrized 2D fluid flow around an obstacle. The numerical experiments demonstrate that our latent DDPM can (i) produce accurate and temporally coherent solutions, (ii) achieve strong generalization capabilities to scenarios involving unseen parameter values, and (iii) extrapolate in time beyond the training horizon.

Stochastic Gradient Langevin Dynamics (SGLD) reduces the cost of Langevin-based sampling by replacing full-dataset drift evaluations with mini-batch approximations, but the resulting subsampling error may offset this computational saving. We study this trade-off for stochastic differential equations with finite-sum drifts and compare the computational cost of SGLD with that of the Euler-Maruyama (EM) method. For a prescribed mean-square accuracy $\varepsilon^2$, we derive complexity estimates that explicitly track the dependence on the dataset size $m$, mini-batch size $s$, and accuracy parameter $\varepsilon$. The resulting comparison reveals distinct parameter regimes in which either method is preferable. In particular, EM can have lower leading-order cost only in a small-data, aggressive-subsampling regime, whereas SGLD is favoured over most of the remaining parameter space. In the practically relevant regime $s \ll m$, the transition between the two methods occurs at the scale $m \asymp \varepsilon^{-1}$. We complement the theoretical analysis with numerical experiments based on a Gaussian Bayesian inference model, which examine the predicted cost regimes together with the underlying discretisation error and variance estimates.

In this work we propose a new existence theorem for numerical-flux functions for supercritical fluids that are pressure-equilibrium-preserving (PEP). In particular, we characterize the existence of consistent numerical-flux functions that satisfy an algebraic PEP property. Our theory links the existence of PEP schemes to geometric properties of the equation of state describing the thermodynamics of the fluid. When these geometric properties fail for a pair of states on the same isobar, no PEP schemes of the considered form can exist on a domain containing those states. In our analysis the equation of state itself can be fully general only needing to satisfy some fundamental thermodynamic principles. Recently, PEP compatibility conditions for general equations of state have been derived [1] that rely on the existence of certain thermodynamic derivatives which are not guaranteed to be defined under fundamental thermodynamic principles. The geometric conditions in our theory do not depend on the existence of these derivatives and recover the recent compatibility conditions in the case that these derivatives are defined. Finally, using numerical experiments we demonstrate for two supercritical fluids that our existence conditions are restrictive and thus that no PEP schemes of the form we consider exist on the domain we specify for these fluids. Using our geometric perspective, we also shed light on mechanisms by which PEP schemes can develop numerical issues, which we also demonstrate using numerical experiments.

We bound the difference between two stochastic Volterra processes with identical Lipschitz coefficients but different kernels. For non-convolution kernels, we establish estimates in $C^0([0,T];L^p(\Omega))$, $p\geq 2$, and for convolution kernels in $L^p(\Omega;L^q(0,T))$, $q \in [1,p]$, and $C^\beta([0,T];L^p(\Omega))$, $L^p(\Omega;C^\beta([0,T]))$, where the range of the Hölder exponent $\beta \in (0,1]$ is the maximal permitted by the regularity of the processes. For the fractional kernel, we then construct Markovian approximations whose error we show to decay as $e^{-a\sqrt{N}}$ in the aforementioned norms, using an $N$-node quadrature based on sinc methods. Numerical experiments for the fractional Brownian motion verify our findings.

The derivative-free loss method (DFLM) is a mesh-free neural network approach for solving partial differential equations by simulating stochastic walkers rather than computing derivatives directly. Although DFLM has previously been applied to the Navier–Stokes equations, extending it to turbulent flows reveals two limitations of its simplest implementation. The walkers do not account for the swirling, directionally stretching motion of the flow, and randomly sampling the walkers used to evaluate the target introduces noise into the learned solution. We address both limitations with an analytic approximation of the walkers' local motion that captures directional stretching and eliminates this sampling noise, while remaining computationally efficient. We study the resulting method as a non-intrusive multiscale solver capable of learning macroscopic flow behavior without resolving every scale on a fine grid or explicitly coupling coarse- and fine-scale models. For two turbulent two-dimensional flows, the proposed method reproduces the macroscopic energy spectrum of fully resolved reference simulations more accurately than the standard method, particularly where directional stretching is strongest, confirming DFLM's efficacy as a multiscale solver for turbulent fluid systems.

We establish the sharp CFL condition for the piecewise constant sparse grid discontinuous Galerkin (DG) method with forward Euler time stepping, applied to transport equations with constant coefficients on periodic domains in arbitrary dimensions. For the transport velocity $\boldsymbol c=(c_1,\ldots,c_d)$ and a uniform mesh of size $h$, we prove that the scheme is $L^2$ stable if and only if $\Delta t \leq {h}/{\max_{1\leq \ell\leq d}|c_\ell|}$, whereas the corresponding full grid upwind scheme is well-known to require $\Delta t\leq h/\sum_{\ell=1}^d |c_\ell|$. The sparse grid discretization therefore enlarges the admissible time step by a factor of ${(\sum_{\ell=1}^d |c_\ell|)}/{(\max_{1\leq \ell\leq d}|c_\ell|)}$, which lies between $1$ and $d$ and reaches $d$ for isotropic transport. The proof of sufficiency relies on projection leakage identities for the multilevel Haar decomposition, which allow the mixed directional terms in the energy estimate to be absorbed by the energy discarded by the sparse grid projection. The proof of sharpness follows from alternating modes in one dimension at the finest level. As a by-product, we obtain explicit formulas for the $L^2$ operator norm and the spectral radius of the amplification operator. For spaces over general downward closed index sets, we derive an explicit sufficient CFL condition and a geometric criterion for its sharpness. Numerical experiments in two and four dimensions confirm the theoretical results.

We develop the convergence analysis of the Interface Control Domain Decomposition (ICDD) method, an overlapping domain decomposition method based on an optimal control framework with Dirichlet interface control functions and interface observation. We consider elliptic problems with possible discontinuous coefficients approximated by $hp-$FEM in each subdomain. When the discretizations are conforming on the overlap between 2D domains, we provide theoretical estimates of the number of GMRES iterations needed to solve the non-symmetric interface Schur complement system associated with ICDD. Our results are obtained by combining novel spectral estimates for the Schur complement matrix of ICDD and classical GMRES convergence theory. We prove that the convergence rate behaves as $\mathcal{O}(\delta^{-1} p^{3/2}\log p)$, where $\delta$ denotes the overlap width and $p$ the local polynomial degree, while remaining independent of the mesh size $h$. Numerical experiments verify the theoretical predictions and show the effectiveness of ICDD in the presence of large coefficient jumps in the computational domain. Since in the conforming case, the considered ICDD formulation coincides with the Substructured Restricted Additive Schwarz (SRAS) method, the analysis also provides convergence estimates for SRAS in two dimensions and, through its known equivalence, for the Restricted Additive Schwarz (RAS) method.

The paper presents a higher-order global multi-trace integral equation for time-harmonic electromagnetic scattering by composite objects. The higher-order multi-trace formulation is preconditioned with a Calderón multiplicative preconditioner using higher-order quasi-Helmholtz projectors. The higher-order quasi-Helmholtz projectors separate the solenoidal and non-solenoidal components of the basis functions. Separate access to the Helmholtz components circumvents the explicit inversion of an ill-conditioned mixed Gram matrix. This enables the application of Calderón multiplicative preconditioning without refining the mesh and resorting to dual basis functions. Furthermore, it enables low-frequency stabilization, as the solenoidal and non-solenoidal components can be rescaled individually. The higher-order quasi-Helmholtz projectors are computed iteratively, enabling iterative solvers to efficiently solve the preconditioned matrix system, yielding accurate solutions in a few dozen iterations. Numerical experiments are conducted for composite dielectric bodies, confirming the effectiveness of the proposed preconditioner for the global multi-trace formulation discretized with higher-order basis functions for dense meshes and at very low frequencies

For high-quality images with sharp edges, a popular choice for edge-preserving regularization is using a general(ized) $\ell_q$-norm of the gradient of the image. This can be implemented efficiently using the $\ell_2$-norm and a sequence of weighted gradients, with weights derived from the current solution estimate. We can solve the resulting sequence of regularized least squares problems using hybrid Krylov subspace methods, which efficiently compute the regularization parameter using the problem projected on the Krylov subspace. However, each update of the regularization operator requires a new Krylov subspace. The majorization-minimization generalized Krylov subspace method (MM-GKS) addresses this problem by using a single, generalized, Krylov subspace (GKS). Unfortunately, for large-scale problems, if convergence is not fast, MM-GKS has overwhelming memory requirements and computational costs. We propose a variant of MM-GKS that alternately compresses and expands the search space while maintaining strict monotonic convergence. We show that our method provably converges to the minimum of the selected functional, even if the search space dimension is kept very small. This substantially improves on previous theoretical results for MM-GKS, where the convergence proof relies on the basis for the solution space (eventually) spanning the full space. We show that our method can solve large-scale problems efficiently both in terms of memory requirements and computational complexity. We further generalize our proposed method to handle streaming problems, where the data is either not all available simultaneously or needs to be treated as such because of the extreme memory requirements. We use numerical examples from image deblurring, dynamic photoacoustic tomography, and streaming X-ray computed tomography (CT) to illustrate the effectiveness of our proposed methods.

We establish two complementary lower-bound mechanisms for individual eigenvalues of the Dirichlet Laplacian. First, energy orthogonality yields a frequency-dependent cap on the Fourier density of a finite spectral projection. Combining this cap with the standard $L^2$ Bessel estimate and a radial-capacity bathtub principle gives, on every open set of finite positive measure in $\mathbb R^n$ with $n\geq2$, \[ \lambda_k\geq c_n(2\pi)^2\omega_n^{-2/n}|\Omega|^{-2/n}k^{2/n}, \qquad \frac{n}{n+2}<c_n<1. \] The constant $c_n$ is characterized by a scalar equation, with $c_2=0.5383068077\ldots$. This estimate preserves Weyl scaling and strictly improves the individual consequence of the Li-Yau sum inequality, although it does not improve the sharp leading coefficient in that sum inequality. Second, we retain part of the spectral deficit discarded when an eigenvalue sum is bounded by its largest term. A lower bound for the counting function, integrated through the exact first Riesz-mean identity, leads to a strictly monotone scalar equation. Its unique positive root is no weaker than the volume-only bound, and we give necessary and sufficient criteria for strict improvement over both that baseline and any independent lower bound. Quantitative estimates of Jiang-Lin provide an explicit implementation on bounded Lipschitz domains. The final comparisons and numerical example distinguish improvements within this framework from stronger estimates available under additional geometric or spectral assumptions.

Wavelets and frames that provide strong guarantees for signal representations are a cornerstone in signal processing theory, but have long been restricted to non-causal settings. Recent work extended wavelet analysis to time-causal systems that operate on streaming signals, with exact reconstruction and time-recursive implementations. However, the frame guarantees that make such representations reliable under noise and quantization are uncharacterized. We derive closed-form frame bounds and conditioning on band-limited domains, together with peak frequencies, constant-Q bandwidth tiling, and spectral decay rates for time-causal bandpass wavelets. Our work provides design rules for the scale ratio and channel count, which we validate numerically against discrete implementations of the causal wavelets in terms of recursive filters.

We derive error bounds for a boundary-stopped positive cubature scheme for degenerate parabolic Hamilton–Jacobi–Bellman equations with Dirichlet data on bounded domains. The scheme balances antipodal branches at their first boundary contacts and uses a control-dependent effective time. Its stopped-barrier remainder is \(O(\tau\Delta t^{\gamma/2})\), and the total mass of interpolated branches is at most \(2\tau/\Delta t\), where \(\tau\) is the effective time. After division by the effective time, these estimates give a consistency bound that is uniform as \(\tau\to0\). Under a common strict boundary barrier and the stated mesh conditions, interior consistency, coefficient shaking, switching, and localized comparison yield \[ \|(u-u_h)^+\|_\infty\le C(\Delta t^{1/4}+h\Delta t^{-1/2}), \qquad \|(u_h-u)^+\|_\infty\le C(\Delta t^{1/10}+h^{1/2}\Delta t^{-1/4}). \] Thus \(\Delta t\asymp h^{10/7}\) gives an \(O(h^{1/7})\) maximum-norm bound. The analysis applies to any fixed centrally symmetric positive degree-two cubature and permits rank-deficient diffusion. Numerical examples examine parallel heat loss to a cold wall and control-dependent, rank-deficient diffusion.

We reconstruct a spatially varying material coefficient in a scalar elliptic equation from finitely many boundary excitations, each producing a full Dirichlet trace. To mitigate the ill-posedness and the cost of repeated PDE solves, we restrict the coefficient to a low-dimensional family, specified analytically or learned from samples, and solve the inverse problem in its latent coordinates. We consider a \(C^1\) parametrization with \(m\) latent coordinates and full-rank derivative at a reference point. If the continuous linearized Neumann-to-Dirichlet map is injective on the corresponding tangent space, at most \(m\) excitations suffice for local injectivity and Lipschitz stability. Convergence of the coefficient sensitivities then transfers this stability to conforming finite element discretizations. For sufficiently fine meshes, the stability constant and neighborhood can be chosen independently of the mesh size. Under a local residual-comparison condition, uniform accuracy of the surrogate forward map yields coefficient-error bounds separating representation error, data noise, finite element error, and surrogate error. Derivative accuracy additionally preserves the surrogate's own local stability. All stability statements are local to a reference coefficient. Two-dimensional numerical experiments combine analytic and learned representations of inclusions and crack-like coefficients with neural forward surrogates. They illustrate latent-space reconstruction, reduced online cost, and further improvement from optional FEM-based refinement.

We design and analyse an arbitrary-order numerical scheme for the Reissner–Mindlin plate problem on general polygonal meshes. The scheme is derived from the Hodge–Laplacian associated with a discrete Bernstein–Gelfand–Gelfand (BGG) twisted complex and exploits a discrete $H_2$-based construction for the transverse displacement. We establish a discrete Korn inequality for the underlying Discrete de Rham method (DDR) spaces and prove the convergence of the method. At the lowest order, the analysis yields an error estimate that is uniform with respect to the plate thickness, showing that the scheme is locking-free. Numerical experiments on several families of polygonal meshes support the theoretical results and illustrate the benefits of the enhanced continuity of the transverse displacement discretisation.

The purpose of this paper is the identification of high-fidelity digital twin data models from numerical code outputs by non-intrusive techniques (i.e., not requiring Galerkin projection of the governing equations onto the reduced modes basis). In this paper the author defines the concept of the digital twin data model (DTM) as a model of reduced complexity that has the main feature of mirroring the original process behavior. The significant advantage of a DTM is to reproduce the dynamics with high accuracy and reduced costs in CPU time and hardware for settings difficult to explore because of the complexity of the dynamics over time. This paper introduces a new framework for creating efficient digital twin data models by combining two state-of-the-art tools: randomized dynamic mode decomposition and deep learning artificial intelligence. It is shown that the outputs are consistent with the original source data with the advantage of reduced complexity. The DTMs are investigated in the numerical simulation of three shock wave phenomena with increasing complexity. The author performs a thorough assessment of the performance of the new digital twin data models in terms of numerical accuracy and computational efficiency.

This paper develops an embedded-hybridized discontinuous Galerkin (E-HDG) method combined with a first-order implicit-explicit scalar auxiliary variable (IMEX-SAV) time discretization for variational initial-state data assimilation governed by the unsteady incompressible Navier-Stokes equations. We adopt an optimize-then-discretize strategy. The spatial discretization uses discontinuous piecewise polynomials of degrees $k$ and $k-1$ for the element velocity and pressure, respectively, a continuous degree-$k$ velocity trace, and a discontinuous degree-$k$ pressure trace. The resulting state and adjoint velocities, as well as the reconstructed initial velocity, are globally divergence-free. The forward IMEX-SAV state scheme is unconditionally energy stable. Under suitable regularity and local-trajectory assumptions whose bounds introduce no explicit negative powers of the viscosity, the mesh-dependent condition $\Delta t\lesssim h^2$, and a contractivity condition on the Tikhonov parameter, we establish local existence and uniqueness of the fully discrete OTD optimality system and Reynolds-semi-robust $L^2$ error estimates of order $O(h^k+\Delta t)$ for the state, adjoint and reconstructed initial velocities; the constants contain no explicit negative powers of the viscosity. Numerical experiments confirm convergence on smooth tests and demonstrate machine-precision discrete incompressibility, effective nonlinear-solver behavior, and stable performance in small-viscosity regimes.

Multi-state HLL-type approximate Riemann solutions of ideal magnetohydrodynamics (MHD) typically assume that the medium within the Riemann fan is incompressible, and thus the slow magnetoacoustic mode cannot be included in their space-time wave configurations for the approximated states. We propose a new strategy to design multi-state HLL-type approximate solutions, allowing for the medium within the Riemann fan to be compressible. In particular, for the complete seven-wave configuration of the MHD Riemann problem, we first estimate the slow magnetoacoustic speeds and a longitudinal flow speed between the slow modes, and then follow the conservation laws and Rankine-Hugoniot jump relations across the fast, Alfvén, and slow waves, to calculate all the intermediate states within the Riemann fan. Moreover, we discuss the solutions when certain wave modes degenerate, ensuring well-posedness and smooth transitions between complete and degenerate wave configurations. Numerical simulations using a Finite Volume (FV) solver show that the proposed approximate Riemann solution is less diffusive than the classic HLLD scheme, particularly for slow mode waves. For example, in a 1D test case with a strong longitudinal magnetic field, the new scheme needs one order of magnitude fewer grid cells compared to the classic HLLD scheme to resolve all wave modes.

We study nonconvex methods for matrix completion, the problem of recovering a low-rank matrix from a subset of its entries. Convex methods achieve sample complexity linear in the matrix dimension and the rank, up to logarithmic factors, whereas global guarantees for commonly used nonconvex methods require a higher polynomial dependence on the rank. We close this gap by analyzing Riemannian gradient descent (RGD) and Riemannian Gauss–Newton (RGN) methods. For an $n\times n$ matrix of rank $r$ with incoherence parameter $\mu$ and condition number $\kappa$, the two methods achieve exact recovery with high probability from $O(\mu nr\log n\log(n\kappa))$ and $O(\mu nr\log n\log(2\mu r\kappa))$ observations, respectively. The methods use a multiscale residual initialization, while the analysis simultaneously controls the spectral error and incoherence. The resulting RGD iterates converge linearly, whereas RGN eventually converges Q-quadratically.

Particle swarm optimization (PSO) is a widely used algorithm featured in many state-of-the-art optimization tool-kits. However, rigorous performance guarantees are still lacking. The standard PSO dynamics do not admit a natural mean-field description, which would provide an avenue for theoretical analysis. By modifying the PSO formulation, one can recover the consensus-based optimization (CBO) algorithm with memory, which admits a mean-field limit and facilitates rigorous convergence analysis. These theoretical guarantees rely heavily on the fact that for CBO, the drift and noise strengths can be chosen independently, whereas they are coupled for PSO. We analyze how the PSO parameter coupling affects existing convergence guarantees for CBO and its variant with memory effect. We show, by an explicit construction, that the coupling still leaves a non-empty set of admissible parameters for these convergence guarantees to hold. However, the admissible parameter ranges shrink in the limits used to recover PSO. The resulting convergence guarantees from CBO therefore do not directly extend to the classical PSO model. We provide numerical simulations illustrating the parameter tradeoffs shown in the theoretical analysis.

Numerical simulations of ocean flows typically require fine computational meshes to resolve the Munk scale, leading to high computational costs. Filtering-based large eddy simulation (LES) provides a way to relax the mesh size requirement by modeling the effects of the unresolved scales. For the implementation of this strategy, we propose a three-step algorithm called Evolve-Filter-Relax (EFR) that requires (i) the solution of a QGE problem, (ii) a nonlinear Helmholtz filter for the potential vorticity field leveraging an indicator function, and (iii) a final relaxation step. We show that the EFR algorithm can be interpreted as a splitting scheme for a perturbed QGE problem with additional dissipation and provide a practical choice for the relaxation parameter. For comparison, we also investigate a nonlinear Bardina regularization of the QGE. Numerical results on a classical benchmark show that both the EFR approach and the nonlinear Bardina regularization significantly improve the accuracy and stability of coarse mesh simulations with no LES model. Additionally, the EFR method with a deconvolution-based indicator function delivered the best balance between accuracy, stability, and computational efficiency in a test case involving a more realistic geometry (Mediterranean Sea).

This paper develops a general framework for the well-posedness, regularity, and strong approximation of the stochastic reaction–diffusion equation (SRDE) with superlinear drift and diffusion coefficients. We first extend the well-posedness results in W. Liu and M. Röckner, J. Funct. Anal., 2902–2922, 2010 and W. Liu, J. Differential Equations, 572–592, 2013 to the case of superlinear diffusion in the Gelfand triple \(V \hookrightarrow H \hookrightarrow V^*\), with \(V\) equipped with the norm \(\|\cdot\|_V\), and derive a moment estimate by establishing a new Itô formula for \(\|X\|_V^p\) with general \(p \ge 2\). We then apply this abstract result to the SRDE, establish higher spatial regularity \(\dot H^{1+\gamma}\) for any \(\gamma \in [0,1]\) whenever the initial datum lies in the same Sobolev space, and obtain temporal Hölder regularity. Finally, we construct a family of tamed finite element methods (tamed-FEMs) for the SRDE under general assumptions on the tamed functions, derive their long-time unconditional stability, and establish optimal strong convergence rates. To our knowledge, this is the first strong approximation result for SPDEs with superlinear diffusion coefficients.

A novel semi-implicit second-order finite volume scheme integrating the scalar auxiliary variable (SAV) approach is proposed for solving the pseudo-time Eikonal equation in wall-distance computation. Unconditional energy stability under zero boundary conditions is rigorously proved, eliminating the dependence of the time step on the grid scale and enabling large-time-step computation. The scheme is a priori accuracy-preserving, and its discretization matrix forms an M-matrix, thereby guaranteeing strict non-negativity of the numerical solution inherently. The framework extends readily to any non-conservative scalar equation and, being independent of the specific finite-volume reconstruction, is compatible with schemes of arbitrary order of accuracy. A vanishing artificial viscosity is introduced to smooth the solution without compromising formal accuracy, and upwinding is incorporated through directional weighting in the weighted least-squares (WLS) reconstruction. Numerical experiments confirm that the scheme achieves the designed accuracy and permits stable computations with uniformly large time steps. For complex configurations such as a three-element airfoil and the three-dimensional ONERA M6 wing, accurate results are obtained on high-aspect-ratio grids, with relative errors in the computed wall distance below 3% relative to the search-based reference, except near geometric singularities.