Papers

The linearized buckling load of a stiffened panel depends on the stress stiffness its shell element assembles. We read it from exported operators against three truncations of one second variation. The classic pass of ANSYS SHELL181 carries a rotation-rotation block pairing the drilling freedom with the bending rotations and its perturbation pass does not; removing the block recovers the perturbation load factor to 0.02%. SHELL281 carries block and couplings in both passes. Abaqus S4 matches the critical mode of the complete second variation to 1.0000 on the translations and its load factor to 1.7%, against 17% and 34% for the other two forms. On an optimized panel under cantilevered shear a 20-node continuum lies 3% to 6% above that form, S4 and SHELL281, 11% and 23% below both SHELL181 passes and 25% above Abaqus S8R, at the finest meshes. On a conventionally stiffened panel the SHELL181 passes stand 1.0% and 3.5% above the complete form, 9% and 21% at half the rib pitch; under a shear flow, on cylinders, open beams and under axial compression the three forms coincide and no pass parts by more than 0.3%.

This 30 years old book is posted for free download by the generosity of Oxford University Press that granted permission to allow access to future generations of young scientists. A printed copy can be purchased from the publisher. The authors received for this book "The 1999 INFORMS Computing Society Prize for Research Excellence in the Interface Between Operations Research and Computer Science" from the Institute for Operations Research and the Management Sciences (INFORMS). For a book review by Janos D. Pinter go to: Journal of Global Optimization, 16: 107-108, 2000.

We use the results of Lorist and Schwenninger, A solution to Crouzeix's conjecture, arXiv:2608.03841v2, https://arxiv.org/abs/2608.03841, to update the 2019 paper of Crouzeix and Greenbaum, Spectral sets: Numerical range and beyond, SIAM J.~Matrix Anal.~Appl., 40(3):1087-1101. For all regions considered in the 2019 paper, we are able to show that the region is a $K$-spectral set with a bound on $K$ that is smaller than that established in the 2019 paper. In particular, we show that various annular regions are $2$-spectral sets and that a more general convex region with a circular hole or cutout is a $4$-spectral set.

A fully discrete local discontinuous Galerkin (LDG) method coupled with an implicit-explicit (IMEX) Euler time discretization is presented and analyzed for a class of high dimensional nonlinear stochastic convection-diffusion equations driven by multiplicative $\mathcal Q$-Wiener noise. The model allows nonlinear leading coefficients, nonlinear convection terms, dissipative source terms, and gradient-dependent noise. The diffusion operator is treated implicitly through the LDG formulation, while the nonlinear convection, lower-order drift, and stochastic terms are evaluated explicitly. The main contribution is a high-moment stability and error analysis for the fully discrete scheme. A central difficulty is that the nonlinear terms lead to pathwise growth factors that cannot be controlled uniformly on the full sample space. To provide the stability and error estimate, we introduce recursively defined nested subsets adapted to the numerical solution. Under the stated stochastic parabolicity and refinement conditions, we prove that these subsets have probabilities converging to one. On these subsets, the numerical solution satisfies high-moment stability, and the fully discrete error converges with order arbitrarily close to $r+1$ in space and $1/2$ in time. We also derive a pathwise error estimate by combining the high-moment error bound with a discrete Kolmogorov argument. Numerical experiments for stochastic Burgers' and Allen-Cahn equations confirm the theoretical rates and demonstrate the robustness of the proposed method for nonlinear stochastic models.

We propose a filtered time stepping finite element scheme for curve shortening flow of open and closed curves in arbitrary codimension that is second-order accurate in time. Open curves are assumed to evolve inside a given domain $\Omega \subset \mathbb R^n$, $n\geq2$, and meet the external boundary $\partial\Omega$ orthogonally. We prove optimal error bounds for the $L^2$– and $H^1$–norms. In practice only a single linear system needs to be solved at each time step. Numerical experiments confirm the accuracy and practicality of the introduced method, including an asymptotic equidistribution property.

We prove fully discrete convergence of the classical Barrett–Garcke–Nürnberg (BGN) scheme for curve-shortening flow of smooth embedded closed planar curves. The main obstruction is that the mass form controls only normal motion, whereas the tangential motion is determined implicitly by the curvature equation and governs the parametrization. The usual length-decay estimate therefore does not control perturbations of the full position update. We separate temporal and spatial errors through the time-semidiscrete BGN solution. A scalar normal resolvent and exact curvature and length identities yield uniform regularity and first-order time convergence. For the spatial analysis, an adapted normal–tangential norm gives a near-contractive estimate for the linearized update. Gauss–Lobatto cancellations and an exact covariance identity for the assembled nodal normals produce an $H^1$ one-step defect of order $h^{k+1}$, while an exact difference identity controls the nonlinear remainder. For every fixed $k\ge1$, including the original piecewise linear method, we obtain the matched $W^{1,\infty}$ error bound $C(\tau+h^k)$ on periodic quasi-uniform meshes with $h\le c\tau^2$. All discrete steps are uniquely solvable, and the numerical curves remain regular and embedded. To the best of our knowledge, this is the first fully discrete convergence result for the classical BGN curve-shortening scheme without additional stabilization.

The transmission eigenvalue problem plays an increasingly important role in inverse scattering theory. Although significant progress has been made in developing efficient numerical methods for the problem over the last two decades, its numerical treatment for curved domains in $\mathbb{R}^d$ ($d=2,3$) remains challenging. In this paper, we introduce and analyze a geometrically flexible and $H^2$-conforming isogeometric method for solving a fourth-order, quadratic and non-self-adjoint eigenvalue problem arising from the Helmholtz transmission eigenvalue problem. Using the spectral approximation theory for compact non-self-adjoint operators, we derive optimal error estimates for the discrete eigenvalues and eigenfunctions. Numerical results for the two- and three-dimensional benchmark problems, including those defined in curved domains, are presented to verify our theoretical results, and to demonstrate the advantages of the method over existing numerical methods in terms of both accuracy and geometric flexibility.

The energy dissipation law and the maximum bound principle (MBP) are two important physical features of the well-known Allen–Cahn equation. In this paper, we develop and analyze novel second-order linear numerical schemes for a class of Allen–Cahn type gradient flows. Our scheme is based on the generalized scalar auxiliary variable (GSAV) approach and a novel second-order exponential time-differencing Runge–Kutta (ETDRK2) method. The resulting formulation overcomes a longstanding difficulty in combining these two techniques while retaining both the MBP and energy stability. We prove that the proposed scheme unconditionally preserves both the MBP and energy stability. In addition, rigorous error analysis is carried out for the proposed scheme, establishing second-order accuracy in both time and space without imposing any coupling condition between the time step $\tau$ and the spatial mesh size $h$. We also present some numerical experiments to demonstrate the efficiency of the proposed scheme and its preservation of the theoretical properties.

This paper proposes and analyzes a high-order finite volume element method on curved-edge quadrilateral meshes for elliptic equations. Unlike existing theories, which are primarily based on straight-edge meshes, this study is the first to establish an analysis of the stability and optimal convergence of the finite volume method on curved-edge meshes. By constructing a dual mesh based on Gaussian points, we overcome the accuracy degradation issues caused by geometric deformation and Jacobian non-uniformity of curved-edge meshes. We prove the coercivity of the discrete bilinear form under weak mesh regularity conditions, thereby obtaining an optimal error estimate in the energy norm. Furthermore, using orthogonality and the Aubin-Nitsche technique, we derive an optimal $L^2$ error estimate. Numerical experiments cover problems with constant and anisotropic coefficients, different dual partition strategies, complex curved boundary domains, and interfaces with large deformations. Numerical results indicate that this method consistently achieves the optimal convergence order in both the $H^1$ and $L^2$ norms on a variety of curved-edge meshes. Compared to straight-edge meshes, curved-edge meshes offer significant advantages in approximating complex curved boundaries and demonstrate better resistance to distortion in cases involving sudden changes in coefficients and large deformations at interfaces. This paper provides a unified theoretical framework for the finite volume element method on curved-edge meshes and verifies the efficiency and robustness of the proposed method.

Discovering the structure and kinetics of chemical reaction networks (CRNs) from time-series concentration data is a fundamental challenge in chemical kinetics, with existing approaches relying on prior mechanistic assumptions or suffering from high computational cost and noise sensitivity. In this work, we present integral chemical reaction neural networks (iCRNNs), a framework that combines the interpretable, physics-constrained, architecture of chemical reaction neural networks (CRNNs) with an integral collocation formulation of the governing dynamics. Rather than solving ODEs forward in time at each training step, we approximate the integral of the learned rate functions directly using numerical quadrature. This yields an entirely algebraic forward pass consisting only of matrix operations, eliminating the repeated adaptive ODE solves of standard CRNN and producing smoother, more predictable training. We provide a recovery error analysis characterising the structural sources of ill-conditioning (conservation laws, reaction reversibility, and shared reactant sets) that limit network identifiability. On two benchmark CRNs, iCRNN trains in roughly half the wall-clock time of the baseline CRNN method on a four-species mechanism and between three and four times faster on a five-species mechanism, while completing every training run reliably and attaining comparable or lower loss.

This paper studies the element-wise convergence behavior of subspace iteration. All results are established for both \(\mathbb{F}=\mathbb{R}\) and \(\mathbb{F}=\mathbb{C}\). For a diagonal matrix \(\Lambda=\mathrm{diag}(\lambda_1,\ldots,\lambda_n)\) with \(|\lambda_1|>\cdots>|\lambda_n|>0\), we analyze the modulus of each entry of the iterative matrix sequence generated by subspace iteration, both without and with the Rayleigh–Ritz procedure. Under mild assumptions on the initial matrix \(X\in\mathbb{F}^{n\times m}\), we first derive exact asymptotic expressions for \(|Q_k(i,j)|\) in the subspace iteration without the Rayleigh–Ritz procedure: entries with \(i\neq j\) decay as \((|\lambda_{\max\{i,j\}}|/|\lambda_{\min\{i,j\}}|)^k\), and the deviation of \(|Q_k(j,j)|\) from \(1\) decays as \(\max\{|\lambda_j|/|\lambda_{j-1}|,|\lambda_{j+1}|/|\lambda_j|\}^{2k}\), with explicit coefficients determined by the LU factorization of \(X\). For the Rayleigh–Ritz variant, we obtain element-wise bounds for the normalized Ritz vectors \(Z_k\). Specifically, for \(|Z_k(i,j)|\), off-diagonal entries with \(i\geq m+1\) decay as \((|\lambda_i|/|\lambda_j|)^k\), off-diagonal entries with \(i\leq m\) decay as \((|\lambda_{m+1}|^2/(|\lambda_i||\lambda_j|))^k\), and the deviation of \(|Z_k(j,j)|\) from \(1\) decays as \((|\lambda_{m+1}|/|\lambda_j|)^{2k}\). These results give an element-wise description of the convergence behavior of subspace iteration. After an orthogonal or unitary change of basis, the results apply to real symmetric or complex normal matrices.

In this article we claim that we are going to give a priori and a posteriori error estimates for a Crank Nicolson type scheme. The problem is discretized by the finite elements in space. The main result of this paper consists in establishing two types of error indicators, the first one linked to the time discretization and the second one to the space discretization.

Fractional partial differential equations couple a memory-dependent time derivative with a nonlocal fractional Laplacian, and their repeated solution under varying fractional orders, initial data, or diffusivity fields is computationally expensive. Neural operators offer a fast surrogate, but standard architectures must relearn the dominant linear fractional evolution from data. We propose the Fractional Fourier Neural Operator (FrFNO), a resolution-robust conditional operator. It injects the analytic Mittag–Leffler response of the frozen linear part as a once-precomputed, resolution-independent propagator table, and trains a spectral convolutional network only on the residual induced by variable coefficients and nonlinear advection, conditioned continuously on both fractional orders. The central theoretical result is that under zero-shot super-resolution the injected propagator fills the out-of-band modes that a standard Fourier neural operator sets to zero, replacing the out-of-band full-field tail by the smaller residual tail (a $K$-independent constant-factor reduction under weak perturbation). The residual and the full solution carry the same Sobolev order, so the advantage is amplitude reduction rather than a steeper tail. Mesh refinement drives the error to a residual-controlled floor. On a nonlinear two-dimensional space–time fractional Burgers problem FrFNO achieves the lowest relative $L^2$ error among six baselines (FNO, PINO, PDNO, CNO, DeepONet, U-Net) at the training resolution and under zero-shot super-resolution, yields the smallest spectral phase error, and remains the best at the integer-order limit and for long integration windows. The same construction extends to coupled fractional systems such as fractional Allen-Cahn and Navier–Stokes equations. Code, training scripts, and reference outputs are available at https://github.com/Derek2021Pang/FrFNO.

We consider variable-density Cahn–Hilliard–Navier–Stokes (CHNS) discretizations on finite element meshes that may change between accepted time levels through fixed-topology motion or topology-changing remeshing. When the discrete spaces vary in time, the phase, kinetic, and pressure histories entering a multistep scheme are measured in different discrete structures and cannot, in general, be transferred by a single operator. We develop decoupled backward Euler (BE) and second-order backward differentiation formula (BDF2) schemes by combining exact physical cross-mesh pairings with history representations compatible with the corresponding phase-energy, kinetic-energy, and pressure-gradient storages. The phase update also determines an Abels–Garcke–Grün-consistent mass flux used in the momentum transport. A scalar capillary-exchange equation separates the phase and fluid solves while retaining the discrete energy exchange. The resulting field subproblems are linear, the scalar equation has a unique positive solution, and the schemes satisfy modified energy balances without a time-step restriction under the stated admissibility assumptions. Numerical experiments confirm second-order temporal convergence under both mesh updates, phase-mass conservation, modified-energy decay in the unforced tests, and comparable Rayleigh–Taylor and rising-bubble dynamics.

A unified kernel-free boundary integral (KFBI) framework is proposed for the Allen-Cahn and Cahn Hilliard equations with homogeneous no-flux boundary conditions on two- and three-dimensional irregular domains. With a stabilized first-order implici-explicit (IMEX) discretization, the Allen-Cahn update reduces to a Neumann modified Helmholtz problem. An auxiliary-variable reformulation reduces the Cahn-Hilliard update to subproblems of the same type without evaluating the Laplacian of the nonlinear source. The Cahn-Hilliard update requires two sequential solves for real shifts, while complex-conjugate shifts allow the real-valued solution to be reconstructed from a single complex solve. All subproblems are handled by the same KFBI solver that indirectly evaluates potentials through equivalent interface problems. To correct discrete mass defects in the Cahn-Hilliard update, a geometry-weighted discrete $L\v{s}$ mass projection is constructed tailored to the KFBI treatment of irregular boundaries. Exterior overlap weights are redistributed to physical nodes, excluding auxiliary KFBI extension values from the mass calculation. Numerical experiments demonstrate second-order spatial convergence, first order convergence of paired temporal-error indicators, and preservation of the prescribed discrete mass. The reported source-free simulations exhibit energy decay, with cumulative projection-induced energy perturbations remaining small relative to the observed dissipation.

We derive a thermodynamically-consistent model of fluid flow through a poroelastic medium. Starting from elastic and fluid free-energy densities, an energy-dissipation rate, and a kinematic constraint, the force-balance equations are derived using variational principles, with the pressure–density constitutive relation emerging as a direct consequence of the variational structure; the same kinematic constraint also supplies the total-flux transport structure. In the ideal-gas limit, the model linearization recovers the classical linear Biot equations. For power-law fluid energies, it yields isentropic pressure–density relations. A key advantage of the variational formulation is that extensions to richer physics, such as thermal effects, chemical reactions, or multi-component fluids, can be incorporated systematically by augmenting the energy and dissipation functionals without redesigning the force-balance or transport closure. We support the model with an energy-compatible two-field discretization and study consolidation under a surface load with three lateral-boundary treatments and three fluid-compressibility exponents.

Discrete-time non-symmetric algebraic Riccati equations (DTNAREs) arise in game-theoretic computations of Nash equilibria. Solving such equations at large scale is often computationally prohibitive. This paper develops a numerical approach for large-scale DTNAREs whose solutions are low rank. A low-rank alternating direction implicit (ADI) method is introduced that recursively constructs a low-rank stabilizing solution without explicitly solving any projected DTNARE. Through the pole-placement property of the low-rank ADI iteration, the method ensures that the implicitly solved projected DTNARE always admits a stabilizing solution. An automatic shift-generation strategy is also developed for the ADI iterations. Once an initial shift is provided, the algorithm computes the low-rank solution without further user intervention. Numerical experiments on DTNAREs with dimensions between \(10^6\) and \(10^7\) demonstrate the accuracy and efficiency of the method. The results confirm that the proposed low-rank ADI algorithm is an effective solver for large-scale DTNAREs that would otherwise be computationally prohibitive.

In this work, we develop a second-order, linear, decoupled, and structure-preserving numerical scheme for the modified Landau–de Gennes model of smectic-A (SmA) liquid crystals. The main contributions are threefold. First, to the best of our knowledge, we propose the first integration of the generalized scalar auxiliary variable (GSAV) approach with a second-order exponential time-differencing Runge–Kutta (ETDRK2) discretization, leading to a second-order GSAV–ETD2 scheme. Second, we prove that the proposed scheme satisfies an unconditional energy-dissipation law, thereby closing the theoretical gap in the energy-stability analysis of second-order GSAV exponential integrators of this class. Third, by deriving a coercive discrete reformulation, we establish a fully discrete error estimate without imposing any coupling condition between $\tau$ and $h$, achieving the optimal convergence rate $\mathcal{O}(\tau^2+h^2)$. Numerical experiments are presented to verify our theoretical results and to simulate the self-assembly dynamics of the SmA phase.

Every smooth closed Riemannian surface of genus zero is isometric to the unit sphere carrying a conformal metric, the round metric multiplied by a positive factor. This paper solves the Poisson equation and the eigenvalue problem of the Laplace–Beltrami operator of such a metric with smooth spherical splines, taking the conformal factor as the only description of the geometry. Because the Dirichlet energy is conformally invariant in two dimensions, the stiffness matrix is the round-sphere matrix for every conformal metric, and the factor enters only the load vector, the mean-value constraint and the weighted mass matrix. A surface embedded in space with a conformal parametrization, a surface of revolution, and a surface given only as a triangle mesh are instances of the same problem that differ in how the factor is obtained; the first two are treated here, and the factor is checked against an area identity before any equation is discretized. The unknown is a spline on a fixed spherical triangulation in broken Bernstein–Bézier form, with $C^r$ smoothness imposed algebraically through edge functionals and the conforming space realized by a null-space matrix computed once and reused for every metric. Each step of the computation is stated with the formulas it needs. Experiments on the round metric, a prolate spheroid, a dumbbell of revolution, an Evans–Fung red-blood-cell profile and a prescribed factor with no embedding show the rates $d+1$ in $L^2$ and $d$ in energy for spline degree $d$, the rate $2d$ for eigenvalues with reproduction of the even-order eigenvalues of the round metric up to degree $d$ to roundoff, and, against parametric surface finite elements on the same triangulations, smaller errors at comparable numbers of unknowns for degrees four and six.

We introduce Bernstein constraint complexes, a representation of finite element differential complexes in which every triangle of a triangulation keeps its own Bernstein–Bézier coefficients and global continuity is imposed only through smoothness functionals attached to edges and vertices. Continuity of scalar splines, tangential continuity of spline vector fields, normal continuity of spline vector fields, and higher smoothness differ only in the edge functionals. No global basis is constructed. The structural result is a coefficient-level commuting identity: on every edge the smoothness functionals of a Bernstein derivative are explicit combinations of the smoothness functionals of its argument, through univariate Bernstein difference matrices. This gives a local proof that the kernels of the smoothness matrices form a subcomplex, for the de Rham traces and for componentwise $C^r$ spline profiles, and a computable row-space test for candidate smooth profiles. We realize the lowest-order Powell–Sabin exact sequence in this form and confirm its exactness by ranks. The constrained Galerkin problem is solved in the broken coefficients with exact smoothness equations; we prove equivalence with the conforming method, treat semidefinite element operators through a pseudoinverse interface system, and give a singular pencil whose finite spectrum is the constrained Maxwell spectrum. Under the covariant and contravariant surface Piola maps the reference functionals of the de Rham sequence transfer unchanged to conforming curved triangulations. Experiments compare the null-space method, the augmented Lagrangian iteration of Awanou, Lai and Wenston, and direct solution of the saddle system on the same problems, and cover planar source and eigenvalue problems, a $C^1$ biharmonic solve, the exact sphere, and an embedded hyperboloid.

Transient pressure diffusion in heterogeneous porous media becomes difficult to resolve efficiently when permeability is discontinuous and spans several orders of magnitude. We develop a domain-decomposed random-feature evolutional deep neural network (DD RF-EDNN) that separates spatial approximation from temporal evolution. Permeability-informed random features are compressed into an orthonormal pressure space, and a conservative finite-volume operator is projected onto this space so that only the reduced coordinates are advanced in time. This formulation retains the dissipative structure of the discrete flow problem without repeated neural-network optimization. The analysis quantifies dictionary and singular-value truncation errors, proves contractivity under an energy-consistent dissipativity correction, and derives a conditional grid-level estimate for homogeneous dynamics and time-independent data admitting a steady lifting. Experiments on low-permeability inclusions, high-conductivity channels and blocks, and three-dimensional grids demonstrate consistent accuracy across distinct coefficient structures. The method reaches a final-time relative $L^2$ error of $9.54\times10^{-4}$ in the principal benchmark and maintains errors on the order of $10^{-2}$ at a permeability contrast of $10^3$. Numerical diagnostics further show that the observed error balance and temporal convergence are consistent with the analysis. These results support DD RF-EDNN as a structure-preserving and interpretable reduced formulation for transient simulations on prescribed heterogeneous media.

GPU implementations of the Muon optimizer and of first-order semidefinite programming (SDP) solvers share one computational pattern: a matrix factorization is replaced by a fixed-depth polynomial filter applied after normalization. When a few dominant spectral components carry most of the input's scale, normalization pushes the remaining spectrum toward zero, where the filter is least accurate. We propose spectral deflation: estimate the dominant components, remove them, filter the normalized residual with the unchanged filter, and restore them. Deflation preserves the target matrix function, and we prove that it strictly reduces the finite-step error of the classical Newton–Schulz family. Implemented with batched randomized SVD for Muon and warm-started subspace tracking for ADMM, deflation consistently improves GPT-2 pretraining over the corresponding Muon baselines with both the Newton–Schulz and Polar Express mappings, also in wall-clock time, and lowers the KKT residuals of factorization-free ADMM on large-scale SDPs within a similar projection time.

A pervasive structural pattern in modern deep learning is the linear factorization block: a submodule of the form $W = BA$ in which two parameter matrices are multiplied directly, with no intervening nonlinearity. Such blocks appear in LoRA adapters, low-rank compressed layers, query-key products of self-attention, and share a common pathology: the factorization is non-unique, which can destabilize training and limit usable learning rates. Despite this, factorization blocks are typically optimized with standard Euclidean methods that ignore the underlying geometry. We introduce Stiefel-AdamW, a near drop-in replacement for AdamW for use wherever such blocks appear. By constraining one factor on the Stiefel manifold while leaving the other Euclidean, Stiefel-AdamW relaxes the full $\mathrm{GL}(\mathbb{R}^r)$ gauge symmetry to a compact orthogonal symmetry, ruling out factor blow-up while retaining the coordinate-wise diagonal preconditioning that gives AdamW its practical strength. Moment estimation is performed in the ambient Euclidean space, with geometry entering only through a tangent-space projection and a manifold retraction. The implementation overhead over AdamW is minimal, and we show that the resulting optimizer inherits both the stability benefits of Riemannian methods and standard convergence guarantees. We validate Stiefel-AdamW on LoRA-style fine-tuning of GPT2, ViT, and Mistral 7B and on full pretraining of GPT2 on OpenWebText, showing consistent improvements over strong baselines at essentially no additional cost over AdamW.

We study uniform approximation of $\exp(tz)$, $t>0$, on $(-\infty,0]$ by rational functions $P_m(z)/(q_m-z)^m$ with a prescribed repeated pole $q_m>0$. A Möbius transformation reduces the problem to polynomial approximation of $F_\lambda(x)=\exp\!\left(-\lambda\frac{1-x}{1+x}\right)$ on $[-1,1]$, where $\lambda=tq_m$. We derive a uniform two-saddle asymptotic formula for the Chebyshev coefficients, including explicit amplitude and phase and a relative remainder for each localized complex saddle contribution, covering fixed, sublinear, and linear pole scalings away from saddle coalescence. Because the two saddle contributions can cancel in a single coefficient, we pass to a growing block of neighboring coefficients and prove that the whole block cannot cancel. This transfers the coefficient asymptotics to approximation errors. For $q_m=(\alpha/t)m$, $0<\alpha<3\sqrt3/2$, the best uniform error has two-sided order $m^{-1/2}H_e(\alpha)^m$; at the optimal ratio $\alpha=1/\sqrt2$ this becomes $m^{-1/2}(\sqrt2-1)^m$. The normalized Chebyshev-weighted $L^2$ projection error has an explicit bounded oscillatory profile. These prefactor-resolved estimates yield a two-term precision-to-work law and quantify mismatch between pole-design and stopping degrees. Finally, the scalar error gives an exact worst-case matrix-action benchmark for self-adjoint negative semidefinite matrices and dimension-independent shift-and-invert Krylov bounds.

We introduce a new class of backward logarithmic orthogonal functions and generalized backward logarithmic orthogonal functions, constructed by applying a terminal-endpoint logarithmic mapping to generalized Laguerre polynomials. These functions are designed for backward spectral approximations of problems whose solutions exhibit weak singularities at the terminal endpoint. The proposed basis functions generate non-polynomial weighted approximation spaces with nodes naturally clustered near the singular endpoint, and therefore provide an effective framework for resolving algebraic and logarithmic endpoint singularities. We develop the basic approximation theory for these backward logarithmic orthogonal functions, including recurrence relations, derivative formulas, orthogonality, Sturm–Liouville characterization, mapped Laguerre–Gauss quadrature rules, weighted projection estimates, backward Lagrange interpolation estimates, inverse inequalities, and stability properties in weighted Sobolev-type spaces defined through a terminal logarithmic pseudo-derivative. A generalized version of the basis is also introduced by incorporating an algebraic scaling parameter, which improves the flexibility of the approximation space and allows singular factors to be represented more effectively. The error analysis and numerical results show that the proposed backward logarithmic basis is particularly suitable for weakly regular functions with terminal-endpoint singularities and can recover exponential or high-order convergence rates that are typically lost when usual polynomial approximations are applied directly.

We investigate the sensitivity of solutions to multi-dimensional scalar balance laws with respect to perturbations of the initial data analytically and numerically. Reliable first-order sensitivity information is essential in gradient-based optimization and inverse problems constrained by hyperbolic balance laws. In hyperbolic problems, such perturbations affect both the smooth components of the solution and the locations of shocks. Consequently, classical difference quotients of the solution operator generally fail to converge in $L^1$, even in one space dimension. Although generalized tangent-vector techniques have been developed for one-dimensional problems, extending these ideas to multiple space dimensions is considerably more challenging because it requires a geometric description of hypersurfaces of discontinuity. We represent perturbations of hypersurfaces of discontinuity by normal displacements and account for the induced variation of the normal direction. This representation yields evolution equations for the components of a generalized tangent vector. Based on this calculus, we develop a numerical method for computing first-order variations with respect to the initial data. Numerical experiments in two space dimensions confirm the expected first-order accuracy of the resulting approximation.

As a fundamental task across computational statistics, scientific computing and machine learning, sampling from high-dimensional probability distributions has received increasing attention in recent years. Numerous sampling algorithms have been proposed, among which underdamped Langevin Monte Carlo (ULMC) methods based on underdamped Langevin dynamics (ULD) have emerged as a class of efficient ones. In this work, we introduce a “universal" predictor-corrector formulation that bridges Euler-type, UBU-type and randomized schemes through different choices of method parameters. Notably, the “universal" integrator induces two novel classes of low-cost integrators, termed low-cost randomized integrators (LC-RIs) and low-cost UBU integrators (LC-UBUIs), as well as their exponential-free variants based on polynomial and rational approximations. The resulting new UBU-type and randomized schemes require only one gradient evaluation and two Gaussians per iteration, considerably reducing the number of gradient evaluations or Gaussians per iteration required by existing counterparts. Further, a general framework of long-time error analysis is developed for general discretization schemes in a probability metric. Under certain smoothness and non-log-concavity conditions, we rely on the unified framework to establish non-asymptotic $\mathcal{W}_1$-error bounds of both old and new schemes, revealing convergence rates of order $\mathcal{O}(d^{\frac{1}{2}}h)$ for Euler-type schemes, order $\mathcal{O}(d h^2)$ for UBU-type ones and order $\mathcal{O}(d^{\frac{1}{2}}h^{\frac{3}{2}})$ for randomized ones. In the strongly convex setting, the same non-asymptotic error bounds can be recovered in $\mathcal{W}_2$-distance. Numerical experiments corroborate the theoretical findings.

Physics-informed neural networks and finite element methods provide two different paradigms for the numerical approximation of partial differential equations: the former are commonly trained by minimizing pointwise strong residuals, whereas the latter are naturally built from weak variational formulations and the finite-dimensional systems obtained after discretization. In this work, we introduce a common framework based on the discretization of functional Gauss–Newton problems by finite families of linear measurements. We show that, through an appropriate duality pairing, the linear measurements can be represented by test functions. The resulting Gauss–Newton system is then precisely a Petrov–Galerkin discretization of the linearized functional problem. This perspective recovers pointwise collocation and natural-gradient constructions as particular cases, while making the choice of test functions an explicit algorithmic design choice. We specialize this framework to elliptic problems, where it naturally leads to weak residual formulations and to a hybrid finite element–neural construction acting on complementary approximation spaces. Numerical experiments support the proposed framework and demonstrate the effectiveness of weak Gauss–Newton formulations and hybrid finite element–neural approximations.

Column subset selection aims to seek a small set of representative columns that accurately approximates a given matrix. We study this problem from the perspective of determinantal point processes (DPPs) and their relation to several existing algorithms. To support incremental column selection, we introduce conditional DPPs, which can be applied to any column selection strategy. In particular, we combine it with adaptive randomized pivoting(ARP), and develop a multi-stage ARP (MSARP) algorithm. We establish theoretical guarantees on the expected Frobenius-norm approximation errors for these methods. In addition, we propose a fixed-precision variant of ARP that adaptively determines the number of selected columns. Numerical experiments on column subset selection and Nyström approximation show that MSARP achieves accuracy comparable to standard ARP while enabling incremental selection.

We design and analyse a hybrid high-order (HHO) method for the linear elasticity problem with a linear slip interface, i.e.\ with interface conditions of spring type, in which the jump of the displacement across the interface is proportional to the traction through a compliancy tensor $\bK=\alpha\bI+(\beta-\alpha)\udl{n}\otimes\udl{n}$. The mesh is fitted to the interface, the discrete unknowns are polynomials of degree $k\ge 1$ on the mesh faces and of degree $k+1$ in the mesh cells, and general polytopal cells are allowed. The two novelties of the method lie in the local symmetric strain reconstruction, which incorporates the interface condition, and in the interface stabilisation, which is built from a regularised interface stiffness $\bS_h$ in the spirit of Hansbo and Hansbo \cite{HH04}. As a consequence, a single formulation covers the whole range of compliancies, from the perfectly bonded interface $\alpha=\beta=0$, where the interface stabilisation acts as a Nitsche-type penalty, to the traction-free interface obtained as $\alpha,\beta\to+\infty$, and no unknown is attached to the jump. We prove that the discrete bilinear form is coercive and that the errors converge as $h^{k+1}$ in the energy norm and as $h^{k+2}$ in the $L^2$ norm, with constants that are independent of the compliancy parameters and of the Lamé coefficient $\lambda$, so that the method is also locking free. Various numerical examples and comparisons are provided to confirm the theoretical results.