Papers

Two-level overlapping Schwarz preconditioners with Generalized Dryja-Smith-Widlund (GDSW) and Reduced GDSW (RGDSW) coarse spaces for high-order finite element methods (FEM) and virtual element methods (VEM) are investigated and applied to the Poisson problem in 2D. The GDSW coarse space is a well-established robust choice for preconditioners in domain decomposition, while the RGDSW coarse space offers a competitive alternative by reducing the coarse problem dimension while maintaining robustness. On the theoretical side, a condition number bound for both preconditioners is proved for conforming $\mathbb{P}_k$ and $\mathbb{Q}_k$ finite element spaces with explicit dependence on the polynomial degree $k$. Through numerical experiments on structured triangular and quadrilateral meshes for FEM and on polygonal meshes for VEM, it is demonstrated that both GDSW and RGDSW maintain scalability properties across polynomial orders up to $k=9$ for selected degrees. The RGDSW approach achieves a reduction in the coarse problem dimension by a factor of approximately three compared to GDSW while maintaining competitive iteration counts.

In this work, we develop low-cost non-iterative Sylvester solvers for dynamical low-rank approximation of reaction-diffusion equations and incorporate them into a basis-update & Galerkin mass-lumped finite element method (\mbox{BUG-MLFEM}). By incorporating mass lumping via Gauss-Lobatto quadrature, the spatial discretization yields a matrix differential equation equipped with diagonal positive-definite mass matrices. We establish the weighted dynamical low-rank approximation (WDLRA) within the Galerkin framework of the MLFEM system, proving the theoretical equivalence between the BUG integrator and MLFEM on low-rank subspaces, which guarantees the discrete well-posedness of the implicit BUG steps. To overcome the severe computational bottleneck of standard $O(m^3+n^3)$ Sylvester solvers, we develop tailored, low-cost non-iterative solvers that exploit the inherent sparse-by-small algebraic structure. Under Dirichlet and Neumann boundary conditions (BCs), the large coefficient matrices are banded, enabling efficient diagonalization of the small $r \times r$ factor matrix and reducing the per-step computational complexity to $O((m+n)r^2)$ and storage to $O((m+n)r)$. For periodic BCs, we combine this diagonalization strategy with a rank-2 Sherman-Morrison-Woodbury update to handle wrap-around corner entries at negligible additional cost. Numerical experiments demonstrate first-order temporal and $(k+1)$th-order spatial convergence of the proposed BUG-MLFEM in the discrete $L^2$ norm, where $k$ denotes the polynomial degree of the finite element space, together with significant computational speedups over existing low-rank integrators.

This paper proposes a Projection Anderson Acceleration (PAA) algorithm for variational inequality (VI) problems, where the underlying function is expressed as a finite sum of Hölder continuous functions and satisfies the $\sigma$-star strongly monotone condition that includes the star strong monotonicity with $\sigma = 1$ as a special case. In the absence of quasi-monotonicity and Lipschitz continuity, we establish the evaluation complexity of the PAA algorithm for finding an approximate solution. Specifically, for the star strongly monotone VI problems, the PAA algorithm needs at most $O (\log (\epsilon^{-1}) \epsilon^{2 (\underline{\alpha} - 1)})$ evaluations of the function and projection for finding an $\epsilon$-approximate solution, where $\underline{\alpha} \in (0, 1]$ is the smallest exponent of Hölder continuity of the component functions. We illustrate our theoretical results by numerical examples of non-quasi-monotone complementarity problems and non-monotone semilinear elliptic equations with non-Lipschitz terms.

We propose a model for the one-dimensional isothermal gas flow through a valve in a pipe. The valve aims at keeping the outgoing flow in an assigned interval; moreover, its opening and closing depend on a further parameter. The model is realized by a system of two conservation laws with suitable conditions coupling the ingoing and the outgoing flow. We study the corresponding Riemann problem, where non-classical shocks show up. Although simple, the model allows to study two key oscillatory behaviors of the flow, which are due, respectively, to the chattering of the valve, especially seen in the numerical simulations, and to the stabilizing hysteresis induced by the valve operation. We characterize the states leading to chattering and show how they depend on the parameters of the model. Through several numerical simulations, we study the interplay between the two above phenomena.

This work proposes a biomechanical registration framework designed to combine functional optical-mapping data with structural high-resolution MRI of cardiac tissue. The two acquisitions are performed in different mechanical configurations: optical mapping is obtained after stretching the sample, whereas MRI is acquired in a fixed resting state. To estimate the deformation between these configurations, the tissue is modeled as a nearly incompressible transversely isotropic hyperelastic material, and the traction amplitude is identified from optical-map marker positions through an inverse problem. Two-and threedimensional simulations show physically admissible deformations and provide a first proof of concept for the proposed strategy.

In this work we introduce a numerical framework for the simulation of biological capsules based on meshes that adapt iteration-wise to the shape of the membrane. The discretization leverages a Hybrid High Order method (HHO) to solve a Stokes problem with a source term distributed on the interface between the internal and the external phases. The HHO discretization supports meshes composed of general polygons, permitting to use fitted meshes that can be easily generated by cutting the elements intersected by the interface. Moreover, polygonal meshes allow arbitrarily accurate interface discretizations without introducing bulk degrees of freedom. The presentation of the simulation algorithm is completed with a proof-of-concept numerical test on a well-established benchmark: the deformation of an elastic capsule immersed in a shear flow.

Immersed-boundary methods avoid body-fitted meshes, but the energy or entropy stability of the underlying flow discretisation does not automatically extend to the immersed treatment. This manuscript analyses the volume-penalised compressible Navier–Stokes equations for both energy and entropy stability, on the continuous and semi-discrete levels, and integrates the isolated source exactly in time. The source is formulated as a one-parameter family indexed by $\theta$, whose momentum component relaxes the fluid velocity towards the prescribed solid velocity; $\theta$ sets how much of the resulting kinetic-energy loss leaves the domain rather than being converted to internal energy. The two stability requirements impose opposing bounds on $\theta$, with dissipation of the relative total energy requiring $\theta\geq0$ and stability of the mathematical entropy requiring $\theta\leq1$. Where penalisation acts, both conditions are therefore satisfied if and only if $0\leq\theta\leq1$, and the symmetric value $\theta=1/2$ divides the kinetic-energy loss equally between the domain and the internal energy. We conclude that any entropy-stable split-form discontinuous Galerkin spectral element method remains semi-discretely entropy stable once the source is added at the collocation nodes, for stationary and prescribed moving bodies alike.

Rayleigh-Bénard convection (RBC) is a model for turbulent flow and challenging to simulate by direct numerical simulation (DNS). Adequate resolution of all relevant length scales requires very fine grids in space and stability restrictions can dictate small step sizes. In this paper, we investigate the use of parallel spectral deferred correction (SDC) as a high-order time-stepping method to reduce wall-clock times of direct numerical simulation of RBC. Performance results on JUSUF (CPUs) and JUPITER (GPUs) with two different codes show parallel SDC can offer some concurrency along the temporal axis plus larger maximum stable time steps. Taken together, this allows for a reduction of wall-clock time required for simulations over one free-fall time by having to compute fewer steps and being able to mitigate the higher cost per time step through temporal parallelization.

We consider multiple saddle-point linear systems with block-arrow structure. We analyze exact and inexact Schur complement based block diagonal preconditioner, by characterizing the eigenvalues of the preconditioned matrices in terms of the zeros of specific rational functions. Numerical results onto large tests coming from liquid crystal modeling and PDE-constrained optimization validate the obtained eigenvalue bounds

We study the recovery of a three-dimensional volume in single-particle cryo-electron microscopy (cryo-EM) from noisy tomographic projections acquired at unknown viewing rotations, drawn independently and uniformly from $\mathrm{SO}(3)$. We show that, in a spherical-shell representation satisfying an explicit shell-separation condition, the first three moments generically determine the underlying real volume up to its intrinsic global $\mathrm{O}(3)$ ambiguity. The proof is constructive and yields a frequency-marching procedure that successively recovers the volume. Consequently, consistent recovery in the high-noise regime has sample complexity $n=\omega(\sigma^6)$, where $n$ is the number of observations and $\sigma^2$ is the noise variance. Numerical experiments validate the constructive recovery scheme under finite sampling and additive noise, and indicate that direct third-moment optimization can succeed beyond the structural assumptions required by the constructive proof.

We study the recovery of a Hamiltonian from discrete trajectory observations using structure-preserving kernel methods. When exact-flow data are fitted through a symplectic one-step scheme, the fixed-step learning target is generally not the physical Hamiltonian itself, but a scheme-dependent inverse modified Hamiltonian. For the implicit midpoint and $p$-implicit symplectic Euler methods, we establish the existence of analytic inverse targets and show that their discrepancies from the physical Hamiltonian are second and first order in the observation step, respectively. Within an RKHS framework, we derive finite-sample reconstruction bounds for random trajectories, allowing dependence among observations from the same trajectory, and separate the error into estimation, approximation, and integrator-induced components. We further establish stability under noisy trajectory observations and develop a multistride Richardson correction that systematically raises the order of the integrator-induced error without requiring additional trajectories. Numerical experiments confirm the predicted scheme-dependent orders, the effect of Richardson correction, and the applicability of the method across several Hamiltonian systems.

Multi-continuum models provide an effective macroscopic description of transient flow in fractured and highly heterogeneous porous media by introducing several interacting state variables at each macroscale location. Although high-order multiscale (HOMS) methods reduce the cost of resolving microscopic heterogeneities while recovering local oscillatory features, repeated solution of the coupled homogenized system and repeated high-order reconstruction remain expensive in long-time and large-scale three-dimensional simulations. This work develops a proper orthogonal decomposition (POD) accelerated HOMS framework for transient multi-continuum problems. A vector-valued POD space is constructed from snapshots of the coupled homogenized solution, and the resulting block system is projected onto a low-dimensional space. Two recovery strategies are considered. The first reconstructs the reduced homogenized fields and performs full HOMS assembly at every output time. The second applies the spatial part of the multiscale recovery operator to retained POD modes offline and evaluates the online solution by direct linear combinations of corrected bases and time-derivative bases. An error estimate between the fully discrete homogenized solution and its POD-reduced approximation is established for the coupled system. A recovered POD-HOMS estimate is then obtained by combining the HOMS consistency error with POD truncation errors propagated through value, gradient, inter-continuum exchange, second-derivative, and time-derivative recovery terms. Numerical experiments for two- and three-dimensional heterogeneous configurations show that the reduced solutions retain the accuracy and local oscillatory behavior of the corresponding HOMS solution while substantially decreasing computational cost.

We develop and analyze a fully discrete semi-Lagrangian method for vector master equations associated with mean field games (MFGs) with individual noise on finite graphs. The method discretizes both the characteristic evolution in time and the population variable on the probability simplex. A central challenge is the degeneracy of the transport coefficients near the boundary of the open simplex, where no boundary condition is prescribed. To overcome this difficulty in the analysis, we introduce an interior localization based on exact density characteristics, which keeps numerical arrivals within compact regions of uniform coefficient control. Combined with a closure argument for the coupled numerical recursion, this localization yields local solvability of the updates and quantitative error bounds on compact subsets of the open simplex. Under additional $C^2$-regularity of the exact solution and suitable space–time mesh coupling, these bounds imply first-order accuracy. Furthermore, the density and value trajectories generated by the numerical master feedback converge to the corresponding MFG characteristics at the same rate. In the potential setting, the associated reconstructed density–flux pair is admissible and asymptotically optimal for the planning problem, with an action gap quadratic in the feedback error. Numerical experiments illustrate the convergence rates, feedback responses to population redistribution, and the effects of graph topology.

Weakly imposed boundary conditions in first-order system least-squares methods go back to Jespersen (1977), but two questions have remained open for decades: the right weights, especially for the Neumann residual, whose natural norm $H^{-1/2}$ suggests $h_F$, and norm and error equivalences, especially for indefinite equations. We answer both for general, possibly indefinite, second-order elliptic equations with non-homogeneous mixed boundary conditions. First, for the weights $h_F^{-1}$ on the Dirichlet and $1$ on the Neumann residual, we prove norm and error equivalences with mesh-independent constants for every discrete pair: the error in a mesh-dependent norm is equivalent to the square root of the functional plus a Dirichlet data term. The least-squares solution converges optimally in $H(\mathrm{div};\Omega)\times H^{1}(\Omega)$, and the functional is a reliable and efficient error estimator, with a recovered flux also for approximations from other methods, even ones violating the boundary conditions. The key idea is to apply the stability of the boundary value problem to pairs that are discrete only in their Dirichlet traces, whose $H^{1/2}$ norm an inverse inequality bounds by the weighted penalty. Second, these weights are the right ones: weaker powers, including $h_F$ on the Neumann boundary, destroy the norm equivalence, and stronger weights are admissible, since the penalty multiplies zero at a comparison function meeting the projected data. Third, for the fully discontinuous counterpart, a pure penalty method with both jumps weighted by $h_F^{-1}$, the error equivalence holds for every discrete pair, also for exact fluxes with normal traces not square integrable on interior faces.

We establish a strong convergence rate for the implicit midpoint method applied to the periodic stochastic Korteweg-de Vries equation with additive noise, addressing a question that has remained open since the work of Debussche and Printems. To control nonlinear error growth, we obtain exponential integrability estimates for the exact and comparison solutions using the coercivity of a classical KdV invariant. We then derive local error bounds from estimates for the comparison solutions and control error accumulation through an $L^2$ error propagation identity that isolates a martingale increment. Under $H_0^{9/2}$ regularity assumptions on the noise and the deterministic initial data, we prove that, on every fixed finite time interval, the maximum $L^2(\mathbb{T})$ error at the time points is $O(h)$ in $L^p(\Omega)$ for every $1\le p<\infty$. Notably, the exponential integrability and the convergence rate hold for every damping coefficient $\lambda\ge0$, including the undamped case $\lambda=0$. Numerical experiments confirm the predicted first-order convergence.

We construct a structure-preserving effective linear approximation of the quantum Boltzmann collision operator for Bose and Fermi gases in the normal phase. For variable hard-sphere kernels, the resulting Fourier spectral evaluation costs $O(N_v^4\log N_v)$, with $N_v$ velocity nodes per direction and no angular quadrature factor. The continuous operator conserves mass, momentum and energy and satisfies a linearized $H$-theorem. We also use the operator at local quantum equilibrium as an approximate collision Jacobian in a modified Newton solver that retains the full nonlinear residual. Numerical tests show smaller errors than the quantum BGK model in homogeneous relaxation and fewer full collision evaluations for steady Couette and Fourier flows than the same Newton solver using a quantum BGK Jacobian.

Artificial dissipation improves the robustness of high-order methods for under-resolved flows, but can compromise accuracy and suppress physically relevant structures if its strength and spatial distribution are not appropriately controlled. We present a sensor-activated nonlinear artificial-dissipation framework, inspired by the residual correction method, for continuous and discontinuous summation-by-parts discretizations. In contrast to many existing spectral-element dissipation operators that either reduce the formal order of accuracy or target only the highest polynomial modes, the proposed framework permits design- and higher-order dissipation with explicit control over its magnitude, spatial distribution, activation, and convergence rate under mesh refinement. This modular construction also recovers several existing methods through appropriate choices of its components. The dissipation preserves discrete conservation and semidiscrete entropy stability whenever the base method satisfies these properties. It is dimensionally consistent, element local, and discretization agnostic, as the added term is constructed independently of the base spatial residual. For smooth solutions under suitable mesh assumptions, we establish parameter choices that retain the formal accuracy of the base discretization. Numerical experiments with summation-by-parts discretizations of the Burgers and Euler equations demonstrate substantially improved robustness for under-resolved flows and flows with shock waves while retaining high-order accuracy for smooth solutions, using the same parameter choices across several test problems. Overall, these results show that targeted dissipation can improve the robustness of high-order methods without excessive damping of resolved flow structures.

An efficient shallow Ritz method was developed in [8, 9] for one-dimensional linear diffusion and diffusion-reaction problems. In this paper, the method is extended to a nonlinear diffusion-reaction equation with a cubic reaction term: (1) a basic error estimate of the shallow Ritz discretization is established and (2) a reduced block Newton (rBN) method for solving the resulting non-convex optimization problem is developed. Again, the computational cost of the rBN method is $\mathcal{O}(n)$ for one-dimensional problems, where $n$ is the number of neurons. Combining a convex splitting scheme for time discretization with the efficient shallow Ritz method, we develop a time-stepping neural network method for solving the Allen-Cahn equation. The unconditional energy stability of the time discretization and an error estimate in the $L^2$ norm are established. Numerical experiments demonstrate that the neural network approach yields accurate solutions for singularly perturbed nonlinear diffusion-reaction equations and the Allen-Cahn equation. The method effectively resolves interior and boundary layers without oscillations and significantly outperforms the finite element method on fixed uniform meshes.

This paper concerns the Helmholtz scattering problem for open waveguide networks, extending the analysis from [arXiv:2302.04353, arXiv:2310.05816, arXiv:2401.04674] to plane-wave incident data. We first reduce the problem to a Fredholm second-kind integral equation on an infinite interface. This produces an integral representation of the solution to the original scattering problem, which is then used to obtain detailed asymptotic estimates of the solution at infinity. To leading order and away from the waveguides, the scattered field is itself a plane wave. Its amplitude and phase shift exhibit a jump across a critical polar angle depending on the direction of the incident plane-wave. Our estimates also imply that the scattered field is outgoing in a weak sense, allowing us to establish a uniqueness result. We conclude by calculating the waveguide-mode contributions to the scattered field in a perturbative regime.

We present a method for enforcing exactly Dirichlet, Neumann, and Robin conditions on general triangular domains with arbitrarily curved boundaries. The algorithmic constructions developed here complements those in a previous work for general quadrilateral domains. Our method combines (i) the exact mapping between a general curved triangular domain and the standard triangle, (ii) TFC constrained expressions, and (iii) Barnhill transfinite interpolations, to provide explicit closed-form trial ansatze that mathematically satisfy the imposed boundary conditions exactly. Corner compatibility induced by the boundary conditions is pivotal, and enforcing exactly the imposed boundary conditions hinges on the exact enforcement of the induced compatibility constraints. The corner compatibility associated with the acute vertices of the standard triangle and the coupling between gradient components when the hypotenuse corresponds to a Neumann/Robin boundary pose particular challenges. Two strategies are critical to resolving these issues, by reformulating conditional corner constraints into equivalent unconditional forms and by introducing a separate TFC-style auxiliary edge construction together with an associated linear solve on the hypotenuse. We analyze the corner compatibility constraints and present formulations for the general form of functions that exactly satisfy the imposed Dirichlet, Neumann, or Robin conditions on curved triangular domains. The proposed method has been implemented using extreme learning machine type randomized neural networks for solving PDE problems. We present ample numerical experiments with linear/nonlinear PDE problems on a variety of domains with complex boundary geometries. Numerical results show that the current method has enforced the Dirichlet, Neumann, and Robin conditions on curved domain boundaries to the machine precision.

Fourier-based spectral methods provide compact global approximations, permit simple differentiation and integration, and achieve spectral convergence for smooth periodic functions. Their use in finite-horizon optimal control problems (OCPs) is nevertheless challenged by the generally nonperiodic state and control trajectories: endpoint mismatch induces artificial discontinuities and can destroy spectral convergence. Leveraging Eckhoff-type endpoint-correction techniques, we develop a Fourier-based spectral framework for OCPs with possibly nonlinear dynamics and constraints. The correction improves convergence without enlarging the set of decision variables, yielding an efficient finite-dimensional numerical method. Furthermore, we formulate selected dynamics and obstacle-avoidance constraints through Sobolev residual bounds, which incorporate information about the residuals over the entire time horizon and improve constraint satisfaction. Numerical experiments on benchmark problems demonstrate competitive computational performance, lower residual errors, and improved solution accuracy.

Fast Fourier Transform (FFT)-based preconditioners provide an efficient approach for approximating inverse operators in large-scale linear inverse problems when the underlying operator structure can be exploited in the Fourier domain. Optical Fourier transforms (OFTs), the analog of digital Fourier transforms, exploit light propagation through optical systems and offer the potential for greater throughput and energy efficiency than conventional digital FFT implementations. However, practical optical systems are inherently imperfect and introduce a variety of errors into the computed output. Thus, OFTs can only provide approximations to their digital counterpart, and the implications of these inexact transforms for iterative solvers remains largely unexplored. This work investigates the feasibility of OFT-based preconditioning of iterative Krylov subspace methods for the solution of discrete inverse problems. We digitally simulate the optical Fourier transform using the FFT and incorporate sources of error expected in practical optical systems into the simulation. Numerical experiments suggest that given sufficient optical system quality, iterative methods utilizing OFT-based preconditioners can maintain reconstruction quality and convergence comparable to those with fully digital FFT-based preconditioners.

Standard semidefinite programming (SDP) solvers working in double precision can become unreliable on difficult problems, even when their output appears plausible. We study this issue for SDP problems arising from high-order moment-SOS hierarchies. Three polynomial optimization problems are analysed by comparing standard double-precision computations with reference solutions obtained using the extended-precision solver \mbox{Loraine\,64x4}. We observe two distinct patterns of numerical failure. In false convergence, the computed bounds form a smooth sequence but approach the limit at the wrong rate. In stagnation, the sequence develops a plateau or oscillations and no longer reflects the convergence of the underlying hierarchy. We introduce an empirical effective condition number based on the positive spectral components of the primal and dual solution matrices. It serves as a warning indicator rather than an error bound, and its growth accompanies numerical deterioration in our examples. We also examine the imbalance between primal and dual spectral conditioning and a Cholesky inverse residuum as potential diagnostic indicators. The results show that solver termination messages and apparently regular convergence curves may be insufficient to establish the reliability of high-order moment-SOS bounds.

We study the complexity of high-dimensional integration of periodic functions from Besov spaces $S^r_{p,q}B(\mathbb{T}^d)$ with mixed smoothness in the quantum model of computation. We construct new quantum algorithms based on the multivariate Faber-Schauder decomposition of such functions and analyze their convergence rates. We show a gain of one convergence order compared to the well-known deterministic rates of convergence. We additionally obtain lower bounds showing the optimality of our results up to logarithms.

We study periodic homogenization of fully-nonlinear second-order Hamilton–Jacobi–Bellman and Isaacs equations subject to a homogeneous Dirichlet boundary condition and the numerical approximation of the effective Hamiltonian, where the diffusion matrix is assumed to satisfy the Cordes condition. In the first part, assuming the Cordes parameter is sufficiently large, we show existence and uniqueness of $W^{2,p}$ solutions to the true cell problem, prove homogenization under minimal regularity assumptions, and establish a quantitative homogenization result. In the second part, we propose and rigorously analyze $W^{2,p}$-conforming and mixed finite element methods for the approximation of the corrector, which are then used in the approximation of the effective Hamiltonian. The key tool used throughout this work is Campanato's concept of near operators.

Anomalous diffusion in heterogeneous media is ubiquitous across physical and biological systems, where spatially varying microstructures give rise to distinct transport regimes. Variable-order nonlocal operators provide a natural modeling framework for describing heterogeneity and transitions between diffusion regimes. However, the position-dependent long-range interactions inherent in these operators impose prohibitive memory costs and make numerical methods developed for constant-order operators computationally impractical. Consequently, accurate and computationally efficient numerical methods for variable-order nonlocal operators remain scarce. In this paper, we propose accurate finite difference discretizations for the variable-order fractional Laplacian $(-\Delta)^{\alpha(\mathbf{x})/2}$ and develop a fast Fourier transform (FFT)-accelerated scheme for their efficient implementation. The computational and memory costs are $O\!\big(MN\log(2N)\big)$ and $O(MN)$, respectively, where $N$ is the total number of spatial points and $M \ll N$. The method is also naturally parallelizable. We establish its consistency and accuracy through a rigorous error analysis that quantifies the local interplay between solution regularity and the spatially varying order $\alpha(\mathbf{x})$. To the best of our knowledge, this is the first work to quantify how these two factors jointly determine convergence rates, with numerical experiments demonstrating this dependence. Extensive numerical experiments confirm the accuracy and efficiency of the proposed method. Applications to anomalous diffusion in composite materials further demonstrate its ability to capture the effects of material heterogeneity.

In many practical applications of Bayesian inference, it is difficult or undesirable to compute the gradient of the negative log-posterior, motivating the study of derivative-free samplers, also called zeroth-order samplers. The aim of this paper is to generalize and rigorously analyze a multiscale derivative-free sampler recently proposed in the literature by Pavliotis, Stuart and Vaes (2022). In the discrete-time setting, the multiscale sampler we study bears a strong resemblance to other recently proposed methods, such as the so-called zeroth-order perturbed stochastic gradient descent algorithm (Jin et al., 2018) and the zeroth-order Langevin Monte Carlo algorithm by Roy et al. (2022). We present error estimates with respect to all the parameters of the problem, both in the continuous-time and discrete-time settings. The error estimates are in three different metrics: strong convergence, convergence of ergodic averages, and uniform-in-time weak convergence.

We present and analyze a spatially and temporally adaptive high-order method for the one-dimensional Richards equation describing vertical water transport in unsaturated soil, with an emphasis on resolving sharp wetting fronts while satisfying a fully discrete water mass balance. The spatial discretization uses an arbitrary-order local discontinuous Galerkin spectral-element formulation with collocated Legendre-Gauss-Lobatto quadrature, and the time integration is based on a high-order Rosenbrock-Wanner method. We focus primarily on the mixed formulation of the Richards equation, which advances water content as the conserved variable while evaluating the nonlinear diffusion operator in terms of pressure head; the resulting fully discrete scheme satisfies the water mass balance in exact arithmetic, independent of spatial resolution and time step size, including under spatial and temporal adaptivity. In contrast, the pressure-head formulation advances pressure head through a nonlinear capacity function and satisfies the corresponding conservation identity only at the semi-discrete level. Numerical experiments examine convergence through spatial and temporal refinement studies, with mass conservation assessed under spatial and temporal adaptivity using the solver flux output method (SFOM). The SFOM results confirm machine-precision mass conservation for the mixed formulation, whereas the pressure-head formulation exhibits mass balance errors above roundoff. Moreover, in the temporal refinement studies, the mixed formulation yields more accurate solutions than the pressure-head formulation for a prescribed time step size. The combination of high-order accuracy, spatial and temporal adaptivity, and fully discrete mass conservation makes the proposed mixed-form discontinuous spectral-element method a promising approach for simulating unsaturated flow in next-generation land models.

We study the boundary corrector $\theta_\epsilon$ in the two-scale expansion for periodic elliptic problems with Dirichlet conditions. By results of Gérard-Varet–Masmoudi and Prange, on smooth uniformly convex domains $\theta_\epsilon$ converges to the solution $\theta^*$ of a homogenized problem with boundary data given by the boundary layer tail $d(n)$, which is determined by a half-space problem for each normal $n$. In the lifted form of Gérard-Varet and Masmoudi, $d(n)$ is given explicitly by the cell correctors and the Dirichlet-to-Neumann operator $N_n$ of this half-space problem. Our main result is that $N_n$ satisfies a nonlinear equation on the period cell alone, so that $d(n)$ and the boundary density that represents it can be computed on the cell. This equation is a Riccati equation in which the derivative in depth is replaced by a commutator with the derivative along $n$. This cell problem can therefore be solved by a Riccati doubling algorithm, and we can compute $d(n)$ at irrational normals without rational approximation. Dropping the commutator term in the cell equation gives an operator square root approximation for $N_n$ that captures the leading term of $d(n)$ for small contrast, for the truncated problem. Finally, we approximate $\theta_\epsilon$ at finite $\epsilon$ by the sum of two solutions of constant coefficient problems, one with oscillatory boundary data and the other $\theta^*$, and test it on the disk.

Recovering piecewise-smooth functions with jump discontinuities from noisy, indirect, and incomplete data is a central task in computational inverse problems. Bayesian neural networks (BNNs) with heavy-tailed weight priors have emerged as effective edge-preserving priors for such problems. However, posterior inference for these models remains computationally challenging, and existing approaches settle for point estimates obtained by directly optimizing heavy-tailed log-densities. We propose GSBL-BNNs, which model the neural network parameters hierarchically using sparsity-promoting priors from generalized sparse Bayesian learning (GSBL). The hierarchical formulation recovers several popular marginal priors, including Laplace and Student's $t$ (including Cauchy). A key modeling ingredient, novel among edge-preserving BNNs and native to the hierarchical formulation, couples each neuron's bias to its weights through a shared hyperparameter, so that the induced jump discontinuities are free to occur anywhere in the domain while the conditionally Gaussian structure underlying efficient inference is preserved. Furthermore, we develop a block-coordinate descent (BCD) algorithm for maximum a posteriori (MAP) estimation. This BCD approach alternates between standard network training with iteratively reweighted $\ell^2$-regularization and exact, closed-form hyperparameter updates for gamma and inverse gamma hyperpriors. Building on the BCD algorithm, we perform approximate Bayesian inference and uncertainty quantification (UQ) using deep ensembles. Computational experiments—including signal denoising, deblurring, and inverting a two-dimensional hyperbolic conservation law— show that GSBL-BNNs deliver accurate, edge-preserving reconstructions with built-in UQ. Compared with Gaussian- and Cauchy-BNNs, they yield more accurate predictive means and avoid overconfident uncertainty estimates.