Papers

We present an intrinsic finite element framework for the numerical approximation of scalar- and vector-valued partial differential equations on general manifolds described by atlases of charts. Weak formulations and their discretisation are expressed exclusively in the flat parametric space of the manifold, where the exact geometry enters only through the metric tensor that is evaluated at quadrature points. Conforming finite element spaces of arbitrary order for the whole de Rham complex are obtained on multi-chart atlases through purely topological inter-chart gluing. We develop a nodal change of basis construction for vector-valued spaces on multi-chart atlases, where nodal degrees of freedom couple charts with incompatible coordinate systems through transmission maps. Inter-chart continuity holds at the nodes with jumps of the order of the approximation error, and the resulting discrete field is exactly tangent to the manifold. Optimal convergence rates for this method are observed numerically for the surface Stokes problem; no a priori error analysis is provided. Numerical experiments on the cubed sphere manifold in two and three dimensions for the Hodge Laplacian problems demonstrate the intrinsic framework is free of the geometric consistency errors for a sufficient quadrature degree, and that intrinsic assembly is substantially cheaper than its extrinsic counterpart. Further considering the rotating shallow water equations that include orography as a perturbation of the geometry demonstrates that mass is conserved exactly, and energy is conserved up to the time discretisation error.

We extend the convergent reaction-diffusion master equation (CRDME) framework on unstructured grids to include two-body interaction potentials between particles. Such interactions are essential for modeling systems in which particles experience electrostatic repulsion, volume exclusion, or van der Waals attraction. Our approach discretizes the weak form of the forward equation for a particle-based stochastic reaction- drift-diffusion model to obtain a new CRDME by combining two complementary techniques: an edge- averaged finite element method with mass-lumping quadrature for spatial transport terms that include potential interactions, and a mix of finite volume and mass-lumping quadrature approximations for reaction terms. This yields transition rates for a spatial jump process of particles hopping between dual mesh voxels and undergoing reactions. The resulting spatially discrete CRDME retains key physical properties of the continuum system in closed domains. It preserves the Gibbs-Boltzmann equilibrium distribution in the absence of reactions, and satisfies detailed balance of both drift-diffusion and reactive fluxes for systems with reversible reactions. Numerical comparisons with Brownian dynamics simulations validate our method, and we demonstrate empirical second-order convergence as the mesh is refined.

We develop a Balancing Domain Decomposition by Constraints (BDDC) preconditioner for the interface systems arising from active-set semi-smooth Newton linearizations of elliptic optimal control problems with box-constrained controls, $L^1$-sparsity, and Moreau–Yosida regularized state constraints. The domain decomposition and BDDC construction are formulated directly at the infinite-dimensional level. Eliminating the interior variables yields a Schur complement for the state and adjoint traces assembled from local subdomain operators. We prove a global–local equivalence result, establish well-posedness of the local and partially assembled problems under verifiable conditions on the primal constraints. The BDDC operator admits a two-level additive Schwarz representation with independent local solves and a finite-dimensional coarse solve. We implement the preconditioner using conforming finite elements and test it for distributed and Neumann boundary control. For a fixed ratio of subdomain diameter $H$ to mesh size $h$, the GMRES iteration counts are nearly independent of mesh refinement and grow only moderately as this ratio increases. Enriching the coarse space substantially improves convergence and reduces sensitivity to the control regularization and sparsity weights. For the state-constrained problem, scaling the Moreau–Yosida penalty parameter proportionally to $h^2$ yields nearly constant GMRES iteration counts as the mesh is refined.

We consider using randomized least squares to construct quadrature rules exact on a subspace of functions. Using new conditions that we call relative admissibility and reference weight concentration, we establish both that the quadrature weights from such a procedure concentrate close to their asymptotic values with prescribed and arbitrarily large probability, and this in turn provides useful finite-sample probabilistic bounds on the stability of the resulting quadrature rules for very general classes of possibly complex-valued functions. Our analysis both significantly generalizes the existing analysis of randomized least squares quadrature construction, and provides new bounds on stability for these rules. These results specialize to existence results for positive quadrature rules, when such rules can be theoretically expected. Because our procedures are formally algorithmic, our analysis is a substantive advance toward computationally constructive generalized Tchakaloff theorems.

Non-invasive pressure field estimation from velocity measurements is a longstanding engineering problem. We propose, analyze, and test a pressure-recovery method that computes a full stress field from velocity measurements, and leaves the pressure estimation as a cheap post-processing step. The method relies on a stress-velocity first order formulation of the Navier-Stokes equations, and we show that the formulation accounts for deviations from incompressibility in the measured velocity field by construction. In addition, we theoretically establish the convergence of the finite element (FE) approximation scheme, the stability of the stress recovery with respect to finite-resolution velocity measurements, and then validate this theory numerically. Our results show that the proposed estimator is robust in convective flow regimes and remains accurate at reduced spatial resolution, improving upon state-of-the-art pressure-recovery strategies.

Quantum state tomography is a fundamental technique for estimating the state of a quantum system from measured data and plays a crucial role in evaluating the performance of quantum devices. However, standard estimation methods become computationally prohibitive as the system size increases due to the exponential growth of the density matrix, describing a quantum state, with the number of qubits. We propose a low-rank tensor-network framework for mixed-state quantum state tomography based on a block tensor train (Block-TT) factorization. Specifically, the density matrix is represented as the contraction of a Block-TT with its Hermitian transpose, yielding a TT analogue of the Burer-Monteiro factorization. This parameterization guarantees Hermiticity and positive semidefiniteness by construction while compressing the number of optimization variables from exponential to linear in the number of qubits. Building on this representation, we develop single-site and two-site density matrix renormalization group (DMRG) algorithms for estimating quantum states from compressed measurements. The resulting methods operate directly on the compressed parameterization, support adaptive rank refinement, and exploit efficient tensor-network contractions for expectation-value evaluation. The framework is applicable to a broad class of low-rank quantum states, including pure states, nearly pure states, and ground states that admit accurate tensor-network approximations. Numerical experiments demonstrate accurate state reconstruction from limited measurements together with substantial reductions in memory requirements and computational cost compared with conventional low-rank tomography methods.

The fixed-charge network flow problem (FCNFP) couples continuous flow allocation with discrete arc-activation decisions, making it a canonical but computationally challenging model for a variety of network design and resource allocation problems. Exact mixed-integer linear programming formulations capture the fixed-charge structure faithfully, but often become difficult to solve on large networks. We propose a scalable continuous-optimization algorithm for large-scale single-commodity FCNFP based on an iteratively reweighted least-squares (IRLS) framework. The method replaces the discontinuous fixed-charge and linear arc cost objective with a smooth nonconvex Lasry–Lions surrogate and solves a sequence of weighted quadratic flow subproblems. Each subproblem is solved by a warm-started dual semismooth Newton method whose Newton systems have weighted graph-Laplacian structure, enabling the use of modern Laplacian solvers. To further improve the discovered arc supports of the challenging underlying combinatorial problem, we also develop an algorithmic variant that incorporates objective-driven perturbation restarts and an anchor-union restricted search that jointly leverages supports discovered by IRLS and by complementary FCNFP heuristics. Computational experiments on 410 benchmark, synthetic, and large-scale instances show that our method obtains the best objective quality among the evaluated scalable FCNFP algorithms, with a mean gap of $1.316\%$ to a time-limited MILP reference and a win-or-tie rate of $90.0\%$ among the non-MILP methods. The results indicate that combining smooth continuous optimization with support-level search is an effective strategy for producing high-quality feasible solutions to large-scale FCNFP.

Graph-data workloads such as diffusion estimation, ranking, semi-supervised learning, and network optimization often solve many Laplacian or symmetric diagonally dominant M-matrix (SDDM) systems with the same coefficient matrix. Approximate Cholesky preconditioners eliminate vertices one at a time and store the resulting sparse approximate factorization, the factor, whose construction cost is amortized across these solves. But eliminating a vertex, the pivot, creates a dense Schur-complement clique among its $d$ active neighbors. We introduce CAST (Canonical Approximate Schur Tree), which replaces this clique with a weighted random spanning tree sampled directly from it. Every realization is connected and contains exactly d-1 edges, while reweighting each selected edge by the reciprocal of its tree-inclusion probability makes the update unbiased. The distribution is independent of the ordering of the pivot neighbors, and we prove that its leverage-score marginals minimize the largest normalized reweighted-edge contribution among unbiased inverse-marginal one-tree estimators. We also introduce CAST-$\rho$, which replaces each pivot neighbor with $\rho$ copies, each carrying a $1/\rho$ share of that neighbor's incident weight, samples a weighted random spanning tree on the expanded clique, and contracts the copies back to the original neighborhood. The resulting update remains unbiased and connected, can be sampled exactly in $O(\rho d)$ time, and satisfies a $1/\rho$ bound on the second moment of the normalized local Schur error. Increasing $\rho$ therefore reduces certified local sampling variability, but may increase construction cost and downstream fill. Empirically, we observe that CAST-1 is the faster default, whereas CAST-2 is preferable when its additional edge contributions remain inexpensive.

As an extension of the Point Estimate Method (PEM) to evaluate probabilistic moments of quantities of interest (QoI) in general $n$-dimensional spaces, the Quadratic Point Estimate Method (QPEM) has been recently developed. This new method is defined to fully represent up to fifth-order input moments in the Gaussian space, providing general analytical expressions for sample locations and weights, without requiring any numerical optimization. The QPEM can significantly improve the estimation accuracy of the output QoI moments, in relation to PEM-based methods whose numbers of sigma points grow linearly with the problem dimension, while at the same time having an affordable and competitive computational cost up to a considerable number of dimensions. The QPEM is further enhanced in this work by enabling copula integration into the framework, which enables effective modeling of the joint input probability density function by estimating marginals and the dependence structure of the involved random variables. The validity and efficient performance of the copula-based QPEM are showcased against numerous other sampling methods in various examples considering two practical scenarios: (i) when the joint dependence structure can be inferred from data, and (ii) when only marginal distributions and correlation matrices are known.

We prove sharp bounds for determining the spectral-type decomposition of Schrödinger operators in the spirit of Smale's program on the foundations of computation. For explicit one-dimensional self-adjoint Schrödinger operators $H=-{\mathrm d^2}/{\mathrm dx^2}+V$ on $L^2(\mathbb R)$, where $V\in C^\infty(\mathbb R;\mathbb R)$ is given by a finite description of all derivatives and derivative bounds, the pure point and absolutely continuous spectral sets cannot, in general, be recovered by any single limiting procedure. The singular continuous spectral set is strictly harder: it cannot, in general, be recovered by two nested limiting procedures. Analytic constructions of dichotomies realize the lower bounds: Gordon-type repetitions for pure point spectrum, high barriers for absolutely continuous spectrum, and an inverse spectral construction for singular continuous spectrum based on Riesz products, moment-killing perturbations, and a computational Gelfand–Levitan scheme. The finite-description framework also implies corresponding limitations on what can be certified in fixed formal systems (e.g., when used in computer-assisted proofs). Conversely, using wavelet-based certified computation, we prove matching upper bounds for broad classes of self-adjoint differential operators on $\mathbb R^d$ with coefficients of locally bounded variation and quantitative local variation control: two limits suffice for the pure point and absolutely continuous parts, and three for the singular continuous part. This provides a sharp hierarchy for spectral types.

Active subspaces identify low-dimensional linear structure in high-dimensional parameter-to-output maps by estimating the dominant eigenspace of a gradient covariance operator. In practice this covariance is replaced by a Monte Carlo estimator built from a limited number of gradient evaluations. Classical analyses based on controlling the covariance error in operator norm lead to sample-complexity estimates that can be substantially more pessimistic than the sampling rules commonly used in computations. This paper studies the empirical active subspace method directly in the projection-error metric relevant for ridge approximation. We derive non-asymptotic quasi-optimality bounds governed by a regularized inverse Christoffel function associated with the gradient field. Under a bounded-gradient assumption, the resulting estimates already improve the sample-complexity estimates obtained from operator-norm covariance bounds. We then show that additional smoothness of the gradient map, expressed through membership in a reproducing kernel Hilbert space, yields sharper coherence estimates and motivates tractable importance sampling from kernel diagonal measures. Furthermore, the same smoothness assumption yields a priori decay bounds for the population active subspace tail energy, which can be combined with our finite-sample estimate to prescribe rank, regularization scale, and sample size, allowing to fully characterize the a priori sample complexity. The abstract assumptions are verified for lognormal Gaussian and affine uniform parametric elliptic PDEs using weighted summability of Hermite and Legendre series expansions.

The classical Melan equation modeling suspension bridges is considered. We first study the explicit expression and global properties of the analytical solution for the simplified “less stiff” model, based on which we derive an a priori estimate for the original classical Melan equation and establish, under an explicit load condition, that it has a unique solution with nonnegative integral under downward live loads, thereby showing the uniqueness of the corresponding deflection curve in the engineering setting. We also develop an efficient iterative approximation method by taking the solution of the simplified “less stiff” model as the first iterate, and prove its geometric convergence with explicit error estimates. The applicability and computational efficiency of the method are demonstrated through calculations for two actual bridges, which also quantify the influence of the nonlinear nonlocal term on the solution and clarify the relationship between the simplified and original models. Several engineering observations are verified and explained, and some related open problems are suggested.

In this paper, we present a family of high-order cut finite element methods based on the discontinuous Galerkin (DG) framework for scalar hyperbolic conservation laws on complex domains. Building on our previous work, we develop a multidimensional formulation that combines macro-element stabilization with flux limiting to obtain a scheme that preserves the maximum principle and remains robust with respect to arbitrary boundary cuts of the background mesh. The physical domain is embedded in a regular background mesh, which may produce arbitrarily small cut cells. To avoid the severe time step restrictions typically associated with such cells, ghost penalty stabilization terms are added on interior facets of macro-elements. The resulting method exhibits stability and accuracy properties similar to those of standard DG methods on fitted meshes. An $L^2$-stability result is derived for the semi-discrete scheme under both periodic and inflow-outflow boundary conditions. To enforce the maximum principle and suppress nonphysical oscillations, we adapt limiter techniques from standard DG methods to the CutFEM setting by defining limiting parameters on macro-elements. In particular, we present a macro-element-based parameterized flux limiter together with adaptations of the Zhang-Shu bound-preserving limiter and the Barth-Jespersen slope limiter. Numerical experiments in two and three spatial dimensions demonstrate optimal convergence orders, preservation of the maximum principle, and accurate shock capturing without spurious oscillations, even for challenging cut configurations involving very small-cut cell intersections.

In this paper, we prove pathwise uniform convergence at rates up to $1/2$ for tamed exponential Euler schemes for semilinear hyperbolic stochastic evolution equations with superlinearly growing nonlinearities and multiplicative noise. We take the term hyperbolic to mean that the leading operator generates a contractive $C_0$-semigroup but no parabolic smoothing occurs. Under local Lipschitz, polynomial growth, coercivity, and monotonicity conditions on the nonlinearities, we establish pathwise uniform strong error estimates of the form \begin{equation*} \Big(\mathbb{E}\max_{0\le j \le N} \|U(t_j)-U^j\|_X^p\Big)^{1/p} \lesssim \sqrt{k} \end{equation*} on a Hilbert space $X$ for $p\in [2,\infty)$. Here, $U$ is the mild solution and $U^j$ is the tamed exponential Euler approximation at time $t_j=jk$ with step size $k>0$. This extends previous convergence results for non-parabolic SPDEs from globally to locally Lipschitz nonlinearities, allowing both drift and diffusion to grow polynomially. In a stochastic Kato framework, we further establish local and global well-posedness as well as uniform a priori estimates for the mild solution and its approximation. Applications to nonlinear stochastic transport, Airy, wave-type, and dissipatively damped nonlinear Schrödinger equations are included, covering different nonlinearities-stopped and fractionally tamed schemes. For the Klein-Gordon equation with cubic velocity damping, this complements previous results obtained for additive noise.

We study the localization of $s$ distinct planar point scatterers from sparsely sampled single-incidence, single-frequency acoustic far-field data over a limited receiver aperture. Under the first Born approximation, the coherent data form a two-dimensional off-grid exponential sum and admit a low-rank Hankel representation. A limited aperture generally admits multiple aperture-admissible Hankel pencils with different lifted dimensions, sampling multiplicities, and Vandermonde conditioning. Hence, the choice of the Hankel pencil affects both data completion and subsequent localization. In this work, we introduce a unified framework for structured completion and off-grid localization based on general aperture admissible Hankel pencils. We establish exact recovery in the noiseless case and stable recovery under bounded noise from a number of random aperture samples that is linear in $s$, up to factors determined by the pencil geometry, source conditioning, and logarithmic terms, together with an explicit localization bound in terms of the completion error. The analysis separates the contribution of the pencil geometry from the conditioning induced by the source configuration and provides verifiable sufficient conditions for stable completion and MUSIC resolution. These results yield quantitative criteria for selecting sample efficient Hankel pencils under a prescribed physical aperture. Numerical experiments on connected and disconnected apertures demonstrate the effects of pencil geometry, sampling, and noise.

Extreme learning machines (ELMs) are a class of machine learning algorithms that randomly intialize and subsequently fix the hidden layer parameters, and solve for the last layer coefficients using a least squares method. ELMs have been successfully applied to the solution of various partial differential equations, but the cost of solving large least squares problems limits their use in large scale problems or when high accuracy is desired. To address this limitation, domain decomposition methods (DDMs) have been developed for ELMs, reducing training times via parallel computation. Yet these methods, even those incorporating a coarse space, have struggled with scalability. This paper introduces a novel nonoverlapping DDM for ELMs that treats and enforces continuity and flux conditions equally. Additive Schwarz preconditioners and GenEO (Generalized Eigenvalues problems in the Overlaps) coarse spaces are constructed to accelerate the method. Numerical experiments show the new method outperforms previous DDMs for ELMs in both iteration count and accuracy, especially for the Stokes and wave equations. Furthermore, weak scalability has been observed in up to 4,096 subdomains when employing GenEO coarse spaces.

Predicting tumor dynamics in biological systems under physiologically relevant conditions using mathematical and computational models remains a challenging problem. Continuum models based on phase-field (diffuse-interface) formulations have proven to be an effective modeling strategy to govern tumor dynamics and interactions of multiple species. Within this framework, the present work investigates a tumor growth model based on the Cahn-Hilliard (CH) equation. The formulation involves a fourth-order differential operator that imposes higher continuity requirement on approximation spaces for a well-defined primal variational formulation. To address this challenge, we use isogeometric analysis (IGA), which inherently satisfies this requirement through spline-based basis functions and eliminates the need for mixed or auxiliary-variable approaches commonly used in standard finite element discretizations. Additionally, a locally adaptive IGA scheme with truncated hierarchical B-splines (THB-splines) is used to reduce computational cost while maintaining accuracy. The model is first evaluated on standard benchmark cases and then applied to an organ-scale, patient-specific geometric model of the breast reconstructed from magnetic resonance imaging (MRI) data. Our results show that the model reproduces known tumor morphologies, ranging from a spheroidal pattern to fingered growth. A series of numerical experiments further shows the diversity of tumor dynamics produced by different model parameter choices. Our findings demonstrate the predictive potential of the CH-based phase-field tumor growth model integrated with a locally adaptive IGA framework.

Parametric PDEs can admit different boundary conditions with different accuracy and computational cost. We introduce a framework for learning when one reduced boundary condition can replace another: paired solutions train a neural network to estimate the resulting domain and boundary errors, and the simpler condition is used only when both predicted errors meet prescribed tolerances. We focus on singular limits in applications, in which a stiff Robin or nonlinear boundary law is replaced by its limiting Dirichlet form. We evaluate the method on a galvanic corrosion problem and other nonlinear stationary and evolution problems.

We present an iteration algorithm for vacuum and Einstein scalar-field equations in double-null gauge, which transform the non-linear PDE into systems of ODE. The numerical realization combines characteristic constraint solves, LGL spectral elements, pole-free spherical operators, Galerkin projection, and independent first-order residual and consistency checks.

We present an energy-stable discontinuous Galerkin (DG) spatial discretization of the intrinsic beam equations with Kelvin-Voigt damping. We reformulate the governing equations by splitting the total sectional resultants into elastic and viscous parts, thereby exposing a mixed system without mixed space-time derivatives. For homogeneous cantilever boundary data, the continuous system obeys an energy dissipation identity. With characteristic upwinding and alternating auxiliary traces the DG discretization is energy stable at the semidiscrete level for the stated homogeneous split boundary closure. The method is implemented in the Julia-based simulation framework Trixi.jl. A physically compatible manufactured solution with all state and viscous components active assesses convergence, while a rotating-beam test verifies the analytic constant-spin branch and the actuated energy balance.

We present a compact implicit high-resolution finite volume scheme for one-dimensional multicomponent transport problems with nonlinear Langmuir sorption. The method combines a second order compact implicit discretization with WENO reconstruction and a time limiter for the treatment of steep gradients and discontinuities. A special attention is given to the Courant number used in the limiter, where local characteristic Courant numbers based on the eigenvalues of the inverse retardation matrix are considered. The compact implicit schemes are compared with analogous explicit first and second order methods. Numerical experiments show second order convergence for smooth solutions and demonstrate that the implicit approach can be used with significantly larger Courant numbers. The method is also applied to a three-component displacement chromatography problem with sharp concentration fronts.

This paper develops a finite element method for wave-structure interaction problems arising in the hydroelastic modelling of submerged elastic plates. The approach is formulated for both two-dimensional settings and three-dimensional channels. For the fluid part of the problem, we use a symmetric discontinuous Galerkin scheme that captures the discontinuity in the potential across the plate. To impose the Sommerfeld radiation conditions, we take inspiration from the Dirichlet-to-Neumann map approach and use the analytic solution in the semi-infinite domains to construct appropriate non-local boundary conditions. For the plate, we use a continuous/discontinuous Galerkin method to resolve the 4$^{\rm th}$ order operator without using continuously differentiable finite elements. We prove well-posedness of the method, showing that the resulting sesquilinear form is bounded and satisfies the Gårding inequality, leading to a unique and stable solution. We implement our formulation as an open-source tool using the Gridap package in Julia. We validate our method against a hypersingular boundary integral method, finding excellent agreement. The flexibility of the methodology is demonstrated through simulations of a two-dimensional slanted plate and three-dimensional plate geometries including centred rectangular plates, offset rectangular plates, and annulus-shaped plates. The finite element method presented here can be readily applied to model wave energy conversion of piezoelectric bimorphs. Furthermore, the method could easily be extended to consider wave scattering from fixed or rigid structures or from variable bottom topography, and is not limited to hydroelastic applications.

This work considers to numerically solve three kinds of reaction-diffusion equations with free boundaries. First, the popular front-fixing method is used to transform the considered free boundary problems into fixed boundary problems. Then, by employing the finite difference method, numerical schemes with $\mathrm{M}$-matrices as their coefficient matrices are developed for these transformed fixed-boundary problems. Next, for the developed numerical schemes, we establish numerical theory involving positivity preservation, monotonicity preservation, and stability. It is noteworthy that, unlike some existing works that impose tight restrictions on the time step-size to discuss the numerical theory, the proposed numerical schemes require only a mild restriction. Finally, numerical examples are provided to test the developed numerical schemes and to confirm the theoretical results.

In this paper, we propose a new one-point iterative method for finding simple roots of nonlinear equations by modifying Halley's method. Its correction function is positive and bounded on $\mathbb{R}$ and matches the local expansion of Halley's correction function near a simple root. The method avoids the singularity and sign reversal of Halley's correction function while retaining at least third-order local convergence. We also establish sufficient conditions for global convergence. Numerical experiments show that the proposed method achieves high convergence success rates and favorable iteration counts over a wide range of initial guesses. Further experiments near the singularity of Halley's correction function demonstrate robust convergence behavior of the proposed method. An application to the van der Waals equation illustrates the effectiveness of the proposed method.

Short-recurrence Krylov methods that are equivalent in exact arithmetic often diverge in floating-point arithmetic. To illustrate this, we provide an exactly representable two-cycle for steepest descent with a recursively updated residual. A complementary convergence theorem gives a sufficient condition under which the stored residual decreases geometrically. We then show that a second positive definite family yields different outcomes for a Hestenes–Stiefel Conjugate Gradient (CG) implementation and a direct Lanczos–Galerkin approach. CG finds the exact solution after four updates, whereas the two-step projected system becomes inconsistent. Under stated perturbation and transfer hypotheses, a common spectral enclosure gives comparable convergence bounds. For a dense, left-to-right evaluation order, we prove a sufficient precision bound for a prescribed backward error within $n$ updates for an $n\times n$ matrix, together with a computable stopping test and specified exponent-range assumptions. We use polynomial estimates and numerical experiments to examine the trade-off between working precision and the number of CG iterations needed to achieve a prescribed backward error. We also compare results on inexact matrix-vector products and preconditioning with the frameworks of Paige and Greenbaum. For block Lanczos, we analyze Householder orthogonalization and singular-value truncation as the block size changes. Under componentwise and normwise error bounds, the computed coefficients satisfy a controlled local recurrence and an exact block Lanczos relation for a nearby symmetric problem in a larger space. A block form of Paige's identity bounds the overlap with a Ritz vector along its residual coordinate direction. An inter-block recurrence describes the evolution of overlap, and a Gram-matrix argument gives a count of additional nearby Ritz values.

We present GridapGeosciences$.$jl, a new parallel distributed-memory Julia package for the numerical approximation of partial differential equations on general manifolds. Our abstract approach to generating discrete manifold domains relies on mesh refinement in general, and is implemented within the Gridap finite element library, whose high-level interface enables seamless definition of both intrinsic and extrinsic variational formulations. GridapGeosciences$.$jl features the cubed sphere manifold, which is used in this study to showcase the functionality of GridapGeosciences$.$jl in different applications. The first application involves scalar transport where we develop intrinsic formulations of common stabilisation terms, which are shown to be equivalent to well-known extrinsic stabilisation formalisations. Such stabilisation terms are used to develop an intrinsic formulation of the two-dimensional thermal shallow water equations that includes upwinded numerical fluxes. The final application is a linearisation of the Boussinesq equations where we demonstrate the capacity of GridapGeosciences$.$jl to capture both surface and vertical dynamics on a three-dimensional manifold.

Connectors are ubiquitous in technical systems and frequently combine mechanical contact with electric current transfer and heat generation. To capture these interacting processes, an electro-thermo-mechanical contact formulation based on the third-medium concept is proposed. Mechanical contact, electric current flow, and heat conduction are represented within a unified finite element framework without explicit contact-surface tracking. Deformation-dependent switching functions govern the onset of electrical and thermal transport during interface closure. Electrical conduction is coupled to the transient thermal problem through Joule heating, while the temperature dependence of the electrical conductivity provides the corresponding feedback on the current flow. The coupled displacement, electric potential, and temperature fields are solved monolithically and combined with an independently interpolated deformation-gradient-like field for robust third-medium regularization under severe compression. Numerical examples demonstrate contact-induced current transfer, subsequent Joule heating and transient temperature evolution, as well as localized current paths during progressive closure of rough interfaces.

Fourier-Motzkin elimination computes an inequality description of the projection of a polyhedron onto a subset of its coordinates by eliminating one variable at a time. It is used in several areas of optimisation and computer science, and it is a standard way of obtaining the entropic constraints of a causal structure, where the marginalisation over the latent variables produces such a projection. Its limitation is the growth of the intermediate systems of inequalities, which can be doubly exponential in the number of eliminated variables even though the projection itself grows only as a single exponential. In practice the computational overload of the method therefore depends on two choices: how the redundant inequalities are removed after each step, and the order in which the variables are eliminated. We consider both. We first show, by an explicit example, that Imbert's redundancy test cannot be interleaved with redundancy removal by linear programming. We show that the two methods, however, can be combined soundly if the derivation records used by Imbert's test are re-initialised after every step at which linear programming is used. We then propose a rule for choosing the elimination order of the variables that gives a significant computational advantage, however, at the cost of increased resource usage. We demonstrate this advantage on some random polytopes, where the rule reduces the running time by factors of between 6 and 25 compared with the same elimination under a fixed order. For entropic descriptions of causal structures, with more than 250 inequalities and more than 100 variables to eliminate, our rule keeps the number of inequalities handled at each step one to two orders of magnitude lower than a fixed order.

Layer potentials provide a classical approach to boundary value problems for Laplace's equation on Lipschitz domains. Kenig's 1994 spectral-radius conjecture for the double-layer operator would ensure operator-norm convergence of the associated Neumann series on mean-zero $L^2$ densities when the boundary is connected. We disprove this conjecture by constructing a bounded simply connected planar Lipschitz domain whose double-layer operator on arclength $L^2$ has essential spectral radius strictly greater than $1/2$. More precisely, for every $t>1/2$ sufficiently close to $1/2$, we obtain such a domain with $\pm i t$ in its Fredholm essential spectrum. The construction starts from smooth graphs whose shapes repeat under translation. In the limit of separated scales, refinement makes solutions of adjoint resolvent equations grow with fixed forcing. The graph slopes remain uniformly bounded. A computer-assisted certificate proves this growth through an inequality for Hermitian $2\times2$ matrices. Its strict margin at $- i/2$ persists at nearby spectral parameters. Normalisation and a Floquet transform then give compactly supported densities with small residuals on the full graphs. We insert rescaled segments of successive graphs into one bounded boundary, where these densities form a weakly null sequence of approximate eigenvectors. The same spectral conclusion holds on a single periodic Lipschitz graph. The certificate combines continuous estimates, exact rational arithmetic and rigorous interval enclosures.

We propose a new fractional threshold-dynamics algorithm for subdiffusive interface motion and apply it to multiphase grain-growth simulations. The method combines the L1 discretization of the Caputo time derivative with a Helmholtz-resolvent thresholding step. At each time level, the L1 convolution produces a convex history average of the previously thresholded states; this history-dependent field is then propagated by a Helmholtz diffusion solve and projected back to pure phases by thresholding. For $\alpha=1$, the method reduces to a Helmholtz-resolvent analogue of the Merriman-Bence-Osher scheme, whereas for $0<\alpha<1$ it introduces a subdiffusive temporal memory. We establish a variational characterization, maximum and comparison principles, and stability for the continuous Helmholtz diffusion step. We also derive a memory-tail estimate and a finite-dimensional pinning criterion under an explicit resolvent-contraction hypothesis. The equal-tension multiphase scheme is implemented on periodic grids. The numerical experiments investigate the joint influence of the L1 history and the smoothing length. The fixed-step runs show faster coarsening for smaller fractional orders in the tested regime, while a comparison at fixed smoothing length examines the influence of the L1 history on successive updates.