Papers

Most of qualitative properties of finite element (FE) approximations for various PDE models, in particular the validity of the discrete maximum principles (DMPs), are usually guaranteed by certain geometric properties (e.g. acuteness/nonobtuseness) of FE meshes used. In this work, we revisit the linear FE discretization of the Poisson equation and discuss, in the context of DMPs, the representation of the stiffness matrix entries in terms of the edge lengths and cotangents of the dihedral angles of the tetrahedral meshes, thus avoiding the common use of face areas and element volumes. Using this type of entry representation, certain structural properties of meshes, and simple trigonometric relations, as proposed in the two-dimensional case by Strang and Fix, we show that DMPs still hold (as the stiffness matrices remain $M$-matrices) for many tetrahedral configurations with obtuse dihedral angles. Such a relaxation of the standard acuteness/nonobtuseness conditions can be vital in practical calculations, as these angle properties are difficult to preserve in global and local mesh refinements. Several representative examples are presented.

A central goal in solving a polynomial eigenvalue problem is to compute eigenpairs that are exact for a nearby polynomial, with each coefficient perturbation bounded by a modest multiple of unit roundoff times that coefficient's norm. A backward stable eigensolver for a linearization does not automatically provide this coefficientwise relative guarantee. Scaling can improve accuracy, but a single scaling may not work across a widely spread spectrum. We construct a family of scaled companion linearizations that covers the nonzero complex plane. Each linearization is paired with a region and an eigenvector recovery rule for which the backward-error amplification factor is bounded at exact pencil eigenvectors, independently of coefficient norm disparities. Tropical roots guide the construction, and their multiplicities reduce the number of required solves. These results provide a basis for coefficientwise relative backward stability across the spectrum.

Deep Gaussian processes (DGPs) use iterative sampling to create non-stationary random fields. This is notoriously expensive because the covariance at each non-initial layer depends on the realised field from the layer below, so it changes with every sample. Moreover, standard fast sampling techniques like circulant embedding or matrix compression fail since the underlying kernel usually lacks the necessary regularity and structure. We develop and analyse a wavelet compression method for convolution Matérn DGPs. At each layer the covariance matrix is compressed using wavelets and a Lanczos algorithm computes the random sample. For specific parameter ranges involving the smoothness of the Matérn kernel as well as the smoothness of the wavelets, we show that the accumulative compression and sampling errors through all layers decay exponentially with increasing truncation parameter. We obtain a practical DGP sampling algorithm by choosing a truncation that gives log-linear covariance storage and matrix-vector costs, enabling DGP sampling at nearly linear computational cost.

We develop novel methods for the fast construction of spectral approximations to smooth linear and weakly singular Volterra integral operators and for the fast evaluation of smooth nonlinear Volterra integrals. The proposed methods employ the Duffy mapping, which transforms the original indefinite integral into a definite integral and maps the original causal triangular domain to a square, thereby avoiding any artificial extension of the kernel function. For smooth linear and weakly singular Volterra operators, the resulting approximations are almost-banded and lower-banded, respectively. The $N\times N$ principal finite section of these approximations can be constructed in $O(N)$ and $O(N^2)$ flops, respectively, via recurrence relations. For Volterra integrals with a smooth kernel that depends nonlinearly on a given function, recurrence relations lead to fast evaluation with a complexity of $O(N^2\log N)$ when the function is approximated by a degree-$N$ Chebyshev approximant. These spectral approximations and evaluation procedures give rise to fast and stable spectral methods for solving Volterra integral equations, as well as a broad class of problems involving Volterra integral operators. Extensive numerical experiments show that these spectral methods significantly outperform existing methods in terms of stability, speed, accuracy, and ease of implementation.

We propose and analyze a projection-based model reduction method for the rapid and reliable solution of time-harmonic acoustic wave propagation problems over a frequency band, with particular emphasis on seismic applications. The high-fidelity model is based on a hybridizable discontinuous Galerkin (HDG) discretization, formulated as a generalized saddle-point system coupling the volume and trace unknowns. The basis construction combines residual-based greedy sampling of the frequency parameter, enrichment with frequency derivatives of the HDG solution, and proper orthogonal decomposition (POD) compression. We establish the well posedness of both the high-fidelity and reduced problems and characterize their stability with inf-sup constants. We investigate the interaction between static condensation, inherent to the HDG discretization, and projection onto the reduced spaces, characterizing its implications for the reduced formulation and approximation errors. Since these operations do not commute in general, we derive a priori error estimates for both the trace and volume variables to quantify the gap between the two operations. We further identify a lifting condition on the volume reduced space that eliminates these gaps and yields a trace error estimate independent of the volume approximation error. Numerical experiments on the two-dimensional Marmousi model and a three-dimensional heterogeneous medium, with more than $10^{7}$ unknowns, illustrate the theoretical estimates and show that the reduced model achieves speedups of several orders of magnitude compared to the high-fidelity solver.

We present eMOOD, an entropy-guided evolution of the Multidimensional Optimal Order Detection (MOOD) framework for high-order finite-volume methods on unstructured meshes. A novel troubled-cell indicator based on numerical entropy production replaces the classical relaxed discrete maximum principle. Using entropy generation as a physically motivated measure of non-smoothness, the detector targets shocks, contact discontinuities and under-resolved flow features with a compact stencil. Compared with conventional MOOD detection, eMOOD reduces flagged cells by up to 75%, improving efficiency and reducing memory-access overhead. The detector is used within both the classical MOOD cascade of progressively reduced-order schemes and a novel MOOD-CWENO hybrid, which retains inexpensive linear reconstruction in smooth regions and activates nonlinear high-order reconstruction only where needed. Implemented in the open-source parallel solver UCNS3D, the methodology is assessed on the Shu-Osher shock-entropy-wave interaction, double Mach reflection, supersonic forward-facing step, shock-helium-bubble interaction and triple material point problems. These benchmarks encompass strong shock interactions, shock-induced instabilities, multi-material interfaces and vortex-dominated flows, testing robustness, shock capturing and accuracy on unstructured meshes. Results show that eMOOD preserves the accuracy and robustness of reference CWENO and MOOD formulations while achieving speed-ups of up to 20% without loss of solution quality. Its compact detection stencil, reduced memory footprint and improved localisation offer a more scalable and efficient framework for massively parallel high-order compressible-flow simulations. UK Ministry of Defence (copyright) Crown owned copyright 2026 AWE.

We approximate the solution operator $(a,f)\mapsto u$ of $-\operatorname{div}(a\nabla u)=f$ by operator networks whose only primitives are heat-semigroup functions of a fixed reference operator and pointwise products. The main result is an emulation theorem on strongly local Dirichlet spaces with a spectral gap and a Bakry–Émery bound $\mathrm{BE}(-K,\infty)$: an explicit scheme reaches $\dot H^1$-error $\varepsilon$ uniformly over coefficients with bounded carré du champ and bounded generator image and over $L^2$ data, with parameters depending only on the coefficient bounds, the accuracy and $K$; on the torus this gives a $C^2$-type coefficient class with constants independent of the dimension. The torus version refines this to Wiener classes of fractional regularity $A^{1+\theta}$, with an explicit rate attained by the scheme's product-defect block. The scheme unrolls the Laplacian-preconditioned Richardson iteration, controls each inexact step in the operator norm of $\dot H^1$ without Galerkin truncation, and synthesises the coupling $\nabla a\cdot\nabla u$ from the product defect of one heat layer. All parameters are explicit; the results are approximation theorems for continuum operators.

We study a spherical perfectly matched layer (PML) method for the time-harmonic Maxwell scattering problem with large wavenumber. Using vector spherical harmonics and a modal transformation for the coupled transverse magnetic components, we establish the sharp wavenumber-explicit stability estimate for the truncated PML problem and prove that its inf-sup constant is of order $O(\kappa^{-1})$. We then analyze the first-order Nédélec edge element method and its continuous interior penalty counterpart. Under the condition that $\kappa^3h^2$ is sufficiently small, the preasymptotic errors in the energy norm and the $\kappa$-scaled $\mathbf{L}^2$ norm are of orders $O(\kappa h+\kappa^3h^2)$ and $O(\kappa^2h^2+\kappa^3h^2)$, respectively. Numerical experiments confirm the theoretical results.

Second-order flows represent an emerging computational framework for approaching stationary points of variational energies, characterized by damped hyperbolic partial differential equations. In this paper, we develop a scalar auxiliary variable (SAV) method for this class of flows, referred to as SAV-E. The resulting scheme is linear and decoupled, and we prove unconditional modified-energy stability. We further establish subsequential convergence to stationary points for the Ginzburg-Landau and conserved Swift-Hohenberg models under suitable assumptions. Numerical experiments support the analysis and demonstrate robustness and efficiency compared with representative baselines and the corresponding gradient flows.

Among the various space interpolation operators used in high-order semi-Lagrangian schemes, the case of finite element interpolation presents a somewhat incomplete stability theory, partly due to the impossibility of applying in a standard form the framework of Von Neumann analysis, even in the constant coefficient case. In this work, the stability of the P_2 semi-Lagrangian scheme is stated in the case of the constant coefficient equation on a uniform mesh, using a suitable extension of the Von Neumann technique to treat block circulant matrices.

Non-periodic functions on closed intervals are extended into sufficiently smooth periodic functions through a continuous Hermite interpolation framework with necessary and sufficient endpoint conditions, using smooth blend-to-zero operators and short extrapolation recurrences. A matrix-free method is introduced for generating these smooth periodic extensions on the fly directly from function values at equidistant points. No precomputed continuation operator or auxiliary basis construction is required. The resulting extensions can then be used for FFT-based differentiation of non-periodic functions. Numerical tests are provided to show the effectiveness of the method, including for the highly accurate solution of a 1D linear advection equation.

This chapter focuses on finite volume methods. The finite volume method is a discretization method that is well suited for the numerical simulation of various types (for instance, elliptic, parabolic, or hyperbolic) of conservation laws; it has been extensively used in several engineering fields, such as fluid mechanics, heat and mass transfer, or petroleum engineering. Some of the important features of the finite volume method are similar to those of the finite element method: it may be used on arbitrary geometries, using structured or unstructured meshes, and it leads to robust schemes. The finite volume method is locally conservative because it is based on a “balance” approach: a local balance is written on each discretization cell that is often called “control volume;” by the divergence formula, an integral formulation of the fluxes over the boundary of the control volume is then obtained. The fluxes on the boundary are discretized with respect to the discrete unknowns.

This work analyzes a $C^1$-conforming virtual element method for the optimal control of generalized Oseen equations, formulated using the stream-function on simply connected polygonal domains. The method is based on a conforming discretization in $H^2$ of arbitrary order $k \ge 2$. The discrete scheme considers an approximation of the state and co-state equations which involves a specifically designed operator for the convective term and a variational discretization of the control variable. We derive quasi-optimal a priori error estimates in the $H^1$ and $H^2$-norms for the stream-function, and in the $L^2$-norm for the control, under standard assumptions on the family of polygonal meshes. The proposed method recovers the fluid velocity, vorticity, and pressure using suitable projections. Finally, we report several numerical experiments to support the theoretical findings using different polygonal meshes including the cavity problem. We highlight that this analysis also covers OCPs governed by Brinkman and Stokes equation.

The relationship between classical methods for solving inverse problems and neural network methods is examined. An answer is provided to the question of what inverse problems can be solved using neural networks. The answer is that neural network algorithms can solve only conditionally well-posed problems (well-posed according to Tikhonov). The reason is that the design of standard neural network algorithms violates the well-known Bakushinsky veto from the theory of ill-posed problems. It requires that any method theoretically guaranteeing a stable solution to an ill-posed problem must depend explicitly on the error level of data. The article highlights the distinction between the concept of regularization (RA) in the theory of ill-posed problems and regularization as understood in neural network training, where Bakushinsky's veto is violated. Due to this violation, any changes to the neural network architecture, regularization during training, or training dataset augmentation cannot alter the answer specified above. It is theoretically demonstrated that an ill-posed inverse problem can be solved only by using a parametric set of neural networks that explicitly depend on the known error levels (as parameters) of data. Accordingly, as an alternative to standard neural network algorithms, \neural network based regularization algorithms (NNBRAs) are proposed to solve ill-posed inverse problems. One approach to their construction involves using a regularized dataset. Such a dataset can be generated from perturbed inverse-problem data using a specific base RA explicitly depended on the error level of the data. Test solving the simplest finite-dimensional ill-posed inverse problem demonstrate that the accuracy of the NNBRA approximate solutions is comparable to that of the base RA. At the same time, errors of conventional neural network methods may be unacceptable.

We establish nearly sharp $W^{1,\infty}$ stability estimates for the Legendre spectral Ritz projection associated with Laplacian on a square domain. Its norm grows proportionally to the square root of the polynomial degree, up to a cubic logarithmic factor. The upper bound uses an integral representation of the discrete gradient projection in terms of one-dimensional semigroup kernels and a rank-two Duhamel comparison between the mass-matrix semigroup and its infinite-dimensional counterpart. The lower bound follows from the endpoint concentration of the Legendre kernel. As an application, we derive gradient maximum-norm error estimates for the Poisson and heat equations, with spatial factors optimal up to logarithmic terms.

We present a CR-CR method for strain gradient elasticity and prove its robustness with respect to the Lame parameter $\lambda$ and the length scale $\iota$. The method introduces the strain as an independent symmetric tensor and discretizes both the strain and the displacement using CR elements. The commuting properties of the CR interpolation operator with the relevant differential operators leads the robustness with respect to $\lambda$. A stabilization term weakly imposes the relation between the strain and the symmetric gradient of the displacement. The method also applies directly to the biharmonic equation and singularly perturbed fourth-order problems. For these problems, we prove that the CR-CR method is equivalent to an RT-CR method proposed by Huang and Tang. This result extends the classical equivalence between the RT-P0 mixed method and the CR method for the Poisson equation to fourth-order problems. The equivalence shows that the stabilization term arises naturally from the RT-CR mixed formulation.

Velocity space discretization is a primary source of computational cost in deterministic methods for multiscale gas flows. A fixed velocity grid often retains redundant kinetic degrees of freedom (DOFs), particularly under near equilibrium conditions. To address this issue, this work develops a nonequilibrium-capturing formulation based on the unified gas kinetic scheme (UGKS), in which the number of kinetic DOFs adapts locally to the degree of nonequilibrium. A discrete orthogonal basis is constructed for numerical implementation, and the velocity distribution function is decomposed into an analytical local equilibrium and a basis expanded nonequilibrium component. The retained, locally adapted DOFs are evolved directly, substantially lowering computational overhead while preserving the discrete conservation of mass, momentum, and energy, as well as accurately capturing nonequilibrium states. The formulation retains the core multiscale coupling of molecular transport and collision inherent to UGKS, providing explicit and implicit schemes for unsteady and steady flows, respectively. Numerical tests across a wide range of flow speeds and rarefaction regimes (from continuum to highly rarefied) demonstrate close agreement with benchmark solutions. For unsteady flows, the explicit scheme achieves speedups of 12.64 to 58.50 over the standard UGKS. In steady state tests, only about $2\%$ of the original kinetic DOFs are required in continuum regions, accelerating each implicit iteration by a factor of 33.87. Even in the most rarefied cases, the scheme captures nonequilibrium behavior using only about $27\%$ of the original DOFs.

We develop and analyze a three-dimensional pure source transfer domain decomposition method (PSTDDM) for the Helmholtz equation in free space. Extending the two-dimensional PSTDDM proposed by [Du and Wu, J. Sci. Comput., 83(2020),83:68], the proposed method transfers sources successively along the $x_1$, $x_2$, and $x_3$ directions, yielding layer-wise, bar-wise, and block-wise PSTDDMs. The main challenge in the convergence analysis is deriving wavenumber-explicit truncation-error estimates. To address this challenge, we establish inf–sup conditions for the uniaxial perfectly matched layer (UPML) problem under relaxed restrictions on the PML parameters and reformulate the source-transfer operators. We then introduce a positive-definite reconstruction step that produces a globally conforming approximation. In addition, a new quantitative lower bound is established for the imaginary part of the complex PML distance, more easily yielding exponential decay estimates in three dimensions. These lead to wavenumber-explicit convergence estimates for all three PSTDDMs. Numerical experiments demonstrate the effectiveness of the proposed methods both as direct solvers and as preconditioners for GMRES.

We propose a potential projection (PP) method for the nonlinear Schrödinger equation (NLSE) with low-regularity and singular potentials by replacing the potential with its Fourier projection (truncated Fourier series). This gives rise to an NLSE with projected potential (NLSEwPP). We rigorously establish its well-posedness under low-regularity potentials, and its convergence to the original NLSE. The NLSEwPP is then directly discretized by the standard time-splitting Fourier pseudospectral (TS-FP) method with the mesh size chosen according to the potential projection parameter, resulting in the fully discrete TS-PP-FP method. Numerically, TS-PP-FP consistently achieves optimal temporal and spatial convergence rates with respect to the regularity of the exact solution under low-regularity and even singular potentials. Compared with the time-splitting Fourier spectral method, TS-PP-FP attains the same convergence rates at significantly lower computational cost and preserves the mass conservation and time symmetry (that are lost in the Fourier spectral discretization of the time-splitting method); compared with the time-splitting FP method, it is substantially more accurate for low-regularity potentials with the same computational cost.

In this paper, we construct normalized periodic kernels from univariate functions satisfying suitable moment conditions and derive sufficient conditions for high-order convolution approximation. Two computable quasi-interpolation schemes are obtained by discretizing the convolution integral using either a full tensor grid or random samples. In each case, the analysis separates the continuous convolution error from the error introduced by discretization. The contribution of kernel energy to the error exposes the dimensional dependence of tensor-product kernels and motivates our weight-adapted ANOVA constructions. Under suitable summability conditions on the weights, we establish convolution error and kernel energy estimates that yield dimension-uniform error bounds for the resulting quasi-interpolants. Numerical experiments corroborate the theoretical convergence rates and illustrate the effect of coordinate weighting on dimensional dependence.

In this work, we develop a novel hybrid Schwarz method, termed as edge multiscale space based hybrid Schwarz (EMs-HS), for solving the Helmholtz problem with large wavenumbers. The problem is discretized using $H^1$-conforming nodal finite element methods on meshes of size $h$ decreasing faster than $k^{-3/2}$ such that the discretization error can be bounded by the data uniformly with respect to $k$. EMs-HS consists of a one-level Schwarz preconditioner (RAS-imp) and a coarse solver in a multiplicative way. The RAS-imp preconditioner solves local problems on overlapping subdomains with impedance boundary conditions in parallel, and combines the local solutions using partition of unity. On the PDE level, this coarse space corresponds to an edge multiscale space proposed in {[S. Fu, G. Li, R. Craster, and S. Guenneau, Multiscale Model. Simul., 19(4):1684-1709, 2021]}. The approximation properties of this coarse space to the discrete solutions from the primal and adjoint problems are derived by a novel conforming Galerkin finite element method based upon a nonstandard variational formulation. This work represents the first theoretical analysis of an edge multiscale space at the discrete level. We prove that the EMs-HS preconditioner leads to a convergent fixed-point iteration uniformly for large wavenumbers, by rigorously analyzing the approximation properties of the coarse space to the global Helmholtz harmonic extension part and to the solution of the adjoint problem. We provide extensive numerical results on the sharpness of the theoretical findings and also demonstrate the method on challenging heterogeneous models.

A sequence of eigenvalue problems requires solving $A_i\vec{x}_i= \lambda_i\vec{x}_i$ where $i \in \{1, \dots, N\}$, $N \gg 1$ and (usually) $\lambda_i$ is the largest or smallest eigenvalue of $A_i$. Often, the matrices $A_i$ come from a discretization of a smooth map depending on a parameter $p \rightarrow A(p)$, that is $A_i = A(p_i)$, so the matrices $A_i$ change gradually along the sequence. Under this latter assumption, we introduce an efficient algorithm for solving sequences of eigenvalue problems. Our algorithm combines the ideas of reduced order model and row subsampling, as done in [E. Jones and Y. Nakatsukasa, arXiv, 2025] for linear systems. In addition, in our approach, the same set of rows is guaranteed to work well independently of the value of the parameter $p$ or the matrix $A(p)$ as long as the snapshot matrix is unchanged and given that the row subset was chosen well at the beginning. Furthermore, we implement an update strategy of the snapshot matrix. Along these results, we present theoretical analysis of the snapshot matrix for eigenvalue problems, showing an exponential decay of the singular values, and we give bounds on the error on eigenvectors for unseen value of the parameter. This analysis generalizes that in [M. Guido, D. Kressner, and P. Ricci, J. Sci. Comput., 2024] and shows that, regardless of the distribution of the snapshot points, the convergence is always (asymptotically) exponential in the number of snapshot points if the function $A(p)$ admits an analytical extension. We further generalize our algorithm for computing largest or smallest singular values and relative singular vectors and to generalized eigenvalue problems. We conclude with numerical results showing the efficiency of our approach.

We construct computable a posteriori error estimators for parametric families of nonlinear elliptic problems by measuring the residual of a variational formulation in the dual norm of a closed subspace of the test space. A single abstract theorem shows that the resulting estimator is bounded above and below by the error in the trial norm, i.e., it is a variationally correct loss function. Two instances are covered. The first considers an $L_2$-like test space, yielding a nonlinear first-order system least-squares loss function that certifies arbitrary, possibly meshless, approximations. The second breaks the test space, yielding a nonlinear discontinuous Petrov-Galerkin loss function computed by element-local problems. All assumptions are verified for quasilinear diffusion, including monotone fluxes and semilinear problems, and a nonlocal Kirchhoff model. Numerical studies train neural networks approximating the parameter-to-solution maps with the latter loss and measure the ratio of loss to squared error.

We propose a method for solving partial differential equations on complex domains by learning smooth, invertible mappings from a fixed reference domain using Neural Spline Flows (NSFs). Monotone rational-quadratic spline coupling layers ensure a strictly positive Jacobian determinant by construction, eliminating the need for Jacobian-sign penalties. A diffeomorphic pullback transfers the weak form of the PDE to the reference domain, enabling finite element analysis on reusable structured meshes. To generalize across shapes, a shared NSF is conditioned on low-dimensional latent codes learned through an auto-decoder framework; a subsequent encoder predicts these codes for unseen geometries in a single forward pass. For a given learned map, the pullback preserves the mathematical structure of the PDE and is exactly equivalent to solving on the domain represented by that map. Because the learned domain generally approximates the target geometry, geometric mapping error remains distinct from finite element discretization and quadrature errors. Numerical benchmarks show that sufficiently accurate learned mappings yield reliable finite element solutions. By avoiding repeated generation of unstructured meshes, the method addresses a bottleneck in simulation workflows. The fixed reference mesh also creates opportunities for efficient assembly, memory-aligned data access, and precomputed quadrature. The resulting framework assigns geometric complexity to the neural network while retaining conventional finite element discretization and convergence properties on the learned domain, with overall accuracy also governed by geometric fidelity.

Because its binary optimization landscape contains exponentially many local optima, MAX-CUT is highly sensitive to initialization. We study Particle Polynomial Relaxation (PPR), a moment-based nonconvex continuous reformulation possessing an advantageous optimization landscape. We warm-start PPR with Combinatorial Multigrid (CMG), a spectral graph-coarsening method that produces structure-aware initial solutions in near-linear time without training data. On t2g20 ($n=400$), benchmark instance with known exact optimum, the combined CMG+PPR pipeline reaches $97.5\%$ of optimum on average and $98.4\%$ on the best run. An extended variant with discrete local-search refinement reaches the published optimum on all sixty Biq Mac instances evaluated by a recent neural branch-and-bound method. Using a sparse CSR representation, the same approach scales to G-set instances with up to $20{,}000$ vertices reducing memory requirements by up to $3333\times$ while maintaining a consistent improvement over the CMG warm start. Gains are largest on sparse weighted graphs; on dense graphs we complement the pipeline with spectral sparsification.

In this paper, we investigate the exponential stability of coupled slow–fast stochastic systems with a small delay. Under an averaged Lyapunov condition, we establish practical $p$th-moment exponential stability for $0<p<2$ using a cutoff-function approach to handle singular Lyapunov derivatives near the origin. Exponential stability in the $p$th moment and almost surely is further obtained through two distinct mechanisms: homogeneous dissipation of the fast dynamics or vanishing sensitivity of the slow coefficients to the fast state near the origin. We extend these results to locally Lipschitz systems, establishing exponential stability with arbitrarily high probability for bounded sets of initial data, thereby broadening the applicability of stability analysis for stochastic systems with delayed feedback. Numerical examples illustrate our results.

We develop fully discrete multi-entropy-stable high-order compact Runge–Kutta discontinuous Galerkin (cRKDG) methods for hyperbolic conservation laws. Our main result is a weak cell entropy inequality for multi-dimensional systems after one full high-order Runge–Kutta (RK) time step. The proof exploits the cell-local intermediate stages of cRKDG: their perturbations are controlled by the initial in-cell variation. The analysis does not require a strong stability preserving decomposition or stagewise entropy correction. For a sufficiently dissipative Lax–Friedrichs flux and explicit RK methods with nonnegative output weights, the inequality holds under a suitable time-step restriction for strongly convex entropy pairs. For a prescribed, uniformly conditioned family of such entropies, common dissipation and time-step conditions ensure weak cell entropy stability for all entropy pairs. A single cell-average-preserving scaling at the end of each step then enforces all corresponding strong cell entropy inequalities simultaneously. We also prove a Lax–Wendroff-type theorem: under the stated assumptions, any strongly convergent sequence of numerical solutions converges to a weak solution satisfying all prescribed entropy inequalities. Numerical experiments in one and two dimensions demonstrate high-order accuracy, the benefit of multiple entropy constraints for nonconvex scalar problems, and robust performance on shock problems.

We consider the problem of estimating the expectation of a function (observable) with respect to a high-dimensional probability distribution, under the assumption that a low-dimensional collective variable map captures all the directions of slow mixing for the distribution. More precisely, under the assumption that the fibers of the collective variable map satisfy a log-Sobolev inequality, we establish non-asymptotic relative-error guarantees for this problem using the eigenvector method for umbrella sampling (EMUS). Our analysis applies to both log-Lipschitz observables and indicator of convex sets, and gives complexity bounds that are polynomial in regularity parameters and inverse relative accuracy, and exponential in the dimension of the collective variable space. As a relative error guarantee, our bounds apply to rare event estimation. We consider Gaussian biasing windows defined on a low-dimensional collective variable space and analyze a modified sparse matrix version of EMUS. The bounds quantify the effect of the map through mixing of the conditional distribution and fluctuation of the observable along the fibers.

This paper extends the non-unique games framework introduced in Bandeira et al. (2020) to address symmetric structures. Precisely, we consider the problem \[ \begin{aligned} \underset{g_1,\ldots,g_n\in \mathcal{G}}{\operatorname{min}}\,\, \frac{1}{|S|}\sum_{i,j=1}^n \sum_{s\in S} f_{ij}(g_i s g_j^{-1}), \end{aligned} \] where $\mathcal{G}$ is a compact group, $S\subset \mathcal{G}$ is a finite subgroup of symmetry elements, and $f_{ij}:\mathcal{G}\rightarrow \mathbb{R}$ is a loss function. A notable example is common lines based orientation estimation of cryo-EM images for molecules with point group symmetry $S$, where $\mathcal{G}$ is the rotation group $\operatorname{SO}(3)$. With a semidefinite programming relaxation, we turn the problem into a convex program that is solved with the alternating direction method of multipliers. Unlike existing approaches which are tailored to specific symmetry groups, ours provides a unified computational framework for addressing all point group symmetries. We demonstrate the applicability of the proposed algorithm on the cryo-EM orientation estimation problem on both simulated and experimental datasets. In particular, our algorithm is effective in generating high-quality ab initio models from only a few class averages computed from experimental images.

We introduce and analyze a two-level overlapping additive Schwarz preconditioner for discontinuous Galerkin discretizations of virtual element spaces on general polygonal meshes in two dimensions, for a second-order elliptic problem whose coefficient is constant on each subdomain and may have large jumps across the interface. The discretization is a weighted symmetric interior penalty formulation: the averages carry the coefficient and the penalty is scaled by its harmonic mean across each edge, which keeps the method consistent and makes the constants of the analysis independent of the coefficient. The coarse space is obtained from a virtual, continuous partition of unity with a single degree of freedom per subdomain vertex, and is therefore applicable to the irregular subdomains produced by a graph partitioner; the discrete harmonic extension used to build it may be replaced by a cheaper projection onto polynomials of low degree. We prove that the condition number of the preconditioned system is bounded by $C(1+\log(H/h))(1+H/\delta)$, where $H$, $h$ and $\delta$ denote the subdomain diameter, the mesh size and the overlap, and where $C$ is independent of the number of subdomains and of the contrast of the coefficient, provided the latter is quasi-monotone. This is the same bound known for the conforming case and the analysis allows for hanging nodes. We also show, by an explicit example, that the quasi-monotonicity assumption cannot be dispensed with. Numerical experiments with structured triangular, hexagonal and Voronoi meshes, and with subdomains obtained both by uniform partitioning and by METIS, confirm the theory.