Abstract This paper proposes a universal algorithm for convex minimization problems of the composite form $g_{0}(x)+h(g_{1}(x),\dots , g_{m}(x)) + u(x)$. We allow each $g_{j}$ to independently range from being nonsmooth Lipschitz to smooth, and from convex to strongly convex, as described by notions of Hölder continuous gradients and uniform convexity. Note that, although the objective is built from a heterogeneous combination of such structured components, it does not necessarily possess smoothness, Lipschitzness, or any favorable structure overall other than convexity. Regardless, we provide a universal optimal method in terms of oracle access to the (sub)gradients of each $g_{j}$. The key insight enabling our optimal universal analysis and a core technical contribution is the construction of two new constants, the approximate dualized aggregate smoothness and strong convexity, which combine the benefits of each heterogeneous structure into single quantities amenable to analysis. As a key application, by fixing $h$ as the nonpositive indicator function, this model readily captures functionally constrained minimization of $g_{0}(x)+u(x)$ subject to $g_{j}(x)\leq 0$. In particular, our algorithm and analysis are directly inspired by the smooth constrained minimization method of Zhang and Lan, and consequently recover and generalize their accelerated guarantees.
Abstract We develop a Galerkin discretization of the electric field integral equation based on a lowest-order virtual element approximation for surface meshes composed of polygons. This new boundary element discretization relies on the divergence-conforming virtual element counterparts of the classical Raviart–Thomas (RT) finite elements on simplices and is particularly suited to handle hanging nodes. We prove the well-posedness of the resulting numerical scheme by establishing the stable-uniform character of the discrete inf-sup condition in the natural norm for polyhedral surfaces. Moreover, we demonstrate through an a priori error analysis the quasi-optimal convergence of the scheme, leading to the same convergence rate as that of the classical RT boundary element scheme. Finally, numerical experiments involving scattering problems are presented in order to give more insight into the behavior of the virtual boundary element scheme in terms of $h$-convergence and accuracy as a function of the regularity of solutions and meshes.
The Epstein zeta function generalizes the classical Riemann zeta function to oscillatory lattice sums in higher dimensions and has recently emerged as a key tool in the simulation of long-range interacting classical and quantum many-body systems. Its computation and analytic properties are therefore of significant interest, yet a rigorous and comprehensive treatment has been lacking. We address this gap by introducing a superexponentially convergent algorithm, complete with error bounds, for computing the Epstein zeta function in any dimension with arbitrary real parameters. Our approach is accompanied by a detailed analysis of the analytic properties of the Epstein zeta function. We first present a concise reformulation of its meromorphic continuation, functional equation and symmetries. We then establish, for the first time, its joint holomorphic continuation in all parameters and offer a complete characterization of the resulting complex singularity structure, which governs convergence rates in numerical algorithms based on the function. Recognizing that the function can be decomposed into power-law singularities and a regularized analytic part we provide an algorithm for removing singularities without cancellation error. This facilitates the evaluation of integrals involving the Epstein zeta function and enables fast precomputations through interpolation methods. It also enables the robust treatment of general Wood-type anomalies in wave scattering problems while avoiding catastrophic cancellation. Finally, it serves as the foundation for a new algorithm for efficiently computing magnetic interactions between arrays of solid bodies. We present the first high-performance implementation for arbitrary real arguments in EpsteinLib, a C library with Python and Julia bindings, and rigorously benchmark its performance and accuracy, achieving full-precision evaluation against known analytic results in dimensions 1, 2, 3, 4, 6 and 8 and against an arbitrary precision implementation across the entire parameter range. Finally, we apply our methods to the computation of quantum dispersion relations in three-dimensional spin systems and Casimir energies in three-dimensional geometries.
Many of the standard matrix factorizations have analogs that are structure-preserving. We present symplectic versions of the LDU, QR, singular value decomposition, and polar decompositions, as well as a symplectic (and constructive) version of the Cartan-Dieudonne theorem. These symplectic factorizations have the additional property that their symplectic factors all transparently have a determinant of +1. Thus, they are 'determinant-revealing', in the sense that they demonstrate the nontrivial fact that the determinant of every symplectic matrix is +1.
Abstract A $C^{0}$-continuous nonconforming virtual element method (VEM) is developed for a boundary value problem arising from strain gradient elasticity (SGE) in two dimensions. The SGE is a fourth-order singular perturbation problem with a small size parameter. Our proposed method is robust with respect to both the size parameter and the Lamé constant. The stability condition of the VEMs is derived by establishing Korn-type inequalities and inverse inequalities. Some crucial commutative relations for locking-free analysis as in elastic problems are derived. The sharp and uniform error estimates with respect to both the microscopic parameter and the Lamé coefficient are achieved in the lowest-order case, which is also verified by numerical results.
Abstract For iteratively computing the smallest eigenpair of a huge-scale symmetric matrix, we construct a randomized admissible block coordinate descent (BCD) method by first partitioning the matrix into a number of blocks with respect to its columns, then computing its next iterate through updating the current iterate along with a randomly selected block sub-vector of the affine coordinate direction, and finally obtaining the step-length through minimizing the Rayleigh quotient of the next iterate. This iteration method is indeed a blockwise variant of the admissibly randomized coordinate descent (CD) method proposed and analyzed recently by Bai & Chen (2025, Admissibly randomized coordinate descent methods for computing extreme eigenpairs of symmetric matrices. Numer. Linear Algebra Appl., 32, e70016:1–15), and it can also be considered as a randomized variant of the block CD method. For this class of iteration methods, we rigorously analyze its local and semilocal convergence properties, and solidly demonstrate its computational advantages over the admissibly randomized CD method, as well as the BCD method by numerical experiments.
We investigate numerical methods for the long-time dynamics of the weakly nonlinear Klein-Gordon equation (NKGE) with a quadratic nonlinearity ($\varepsilon <^>{p}u<^>{p+1}$, $p=1$), where the nonlinear strength is characterized by a dimensionless parameter $\varepsilon \in (0, 1]$. Different from previous studies, we consider a Gautschi-type exponential wave integrator Fourier pseudo-spectral method for the NKGE, which is time symmetric and energy-preserving. For the quadratic power nonlinearity ($p=1$), we establish improved error bounds for the proposed numerical scheme at $O(h<^>{m} + \varepsilon au <^>{2})$ up to the long time $O(1/\varepsilon )$, where $m$ depends on the regularity of the exact solution. Here and below, $h$ is the spatial mesh size and $ au$ is the time step size. This is in contrast to the classic cubic nonlinearity case $p=2$, where only uniform error bounds $O(h<^>{m} + au <^>{2})$ hold up to the long-time $O(1/\varepsilon <^>{p})$. The regularity compensation oscillation technique is highly involved in the error analysis, which has been developed recently to analyze the accumulation of errors carefully. For the first time, we report that the improved error bounds hold for Gautschi-type methods, while only time splitting methods and Deuflhard-type (or Lawson-type) exponential wave integrators are known to enjoy improved error bounds (w.r.t. $\varepsilon$) in literature. Our results indicate that the improved error bounds are not only related to numerical discretizations, but also related to the nonlinear structure. Numerical experiments are presented to verify theoretical findings.
This work concerns the numerical approximations of the boundary conditions for linear hyperbolic systems. Discrete boundary procedure computations and an artificial viscosity are proposed to corrected any three-point finite volume schemes. The fully discrete stability of the resulting numerical schemes is established. Numerical test cases are performed to illustrate the relevancy of the proposed procedure.
In this paper we consider the conforming finite element (FE) approximation of Maxwell's problem and analyse the prescription of essential boundary conditions in a weak sense using Nitsche's method. To avoid indefiniteness of the problem, the original equations are augmented with the gradient of a scalar field that allows one to impose the zero divergence of the magnetic induction, even if the exact solution for this scalar field is zero. Two FE approximations are considered, namely, one in which the approximation spaces are assumed to satisfy the appropriate inf-sup condition that render the standard Galerkin method stable, and another augmented and stabilised one that permits the use of FE interpolations of arbitrary order. Stability and convergence results are provided for the two FE formulations considered.
In this article, we analyze a semi-implicit noniterative decoupling method for the unsteady Navier-Stokes-Darcy system, based on multi-step backward differentiation schemes for the time discretization and finite elements for the spatial discretization. The results in the previous time steps are utilized to directly predict the interface information for decomposing the Navier-Stokes and Darcy sub-domains without any iteration at each time step. A semi-implicit scheme is proposed to linearize the nonlinear convection. In order to analyze the convergence of the finite-element solution of the proposed method, we rigorously prove the error estimate in the $L<^>{2}$ norm for the joint Stokes-Darcy Ritz-projection, notably without relying on the $H<^>{2}$ regularity assumption of the elliptic problem. The general $k$-step backward differentiation formulae ($1\le k\le 5$) are analyzed within a general framework utilizing multiplier techniques, and mathematical induction is employed to address the terms arising from nonlinear advection in the error estimate. Numerical examples are provided to verify the theoretical conclusions and illustrate the proposed method.
We present a finite element method based on a predictor-corrector scheme for the numerical approximation of the Ericksen-Leslie equations, a model for nematic liquid crystal flow including a nonconvex unit-sphere constraint. As a predictor step, we propose a linear, semi-implicit finite element discretization that naturally induces nodal orthogonality between the approximate director field and its time derivative. Afterwards, an explicit projection onto the unit-sphere is applied at every node without increasing the energy. Discrete well-posedness results and energy laws are established. Conditional convergence of the approximate solutions to energy-variational solutions to the Ericksen-Leslie equations is shown for a time-step restriction. Computational studies indicate the efficiency of the proposed linearization and the improved accuracy resulting from the inclusion of a projection step in the algorithm.
This paper develops a family of any order finite element solvers for both linear and nonlinear poroelasticity problems. The backward differentiation formulas (BDFs) are used for temporal discretizations, whereas the weak Galerkin (WG) finite elements on quadrilaterals are utilized for spatial discretizations. For the latter, the discrete weak gradients of shape functions are established in the vector- or matrix- Arbogast-Correa spaces. Such combinations of BDF and WG discretizations produce numerical solutions that have well-balanced spatial and temporal errors in all six quantities, namely, displacement, dilation (divergence of displacement), stress, pressure, velocity and normal flux. For nonlinear poroelasticity problems in which permeability depends on solid dilation or the mean stress, Picard iterations are employed to solve the resulting nonlinear algebraic systems. Rigorous analysis together with numerical experiments demonstrate that these new solvers have optimal-order convergence and are free of locking.
The generalized minimal residual methods (GMRES) for the solution of general square linear systems is a class of Krylov-based iterative solvers for which there exist backward error analyses that guarantee the computed solution in inexact arithmetic to reach certain attainable accuracies. Unfortunately, these existing backward error analyses cover a relatively small subset of the possible GMRES variants and cannot be used straightforwardly in general to derive new backward error analyses for variants that do not yet have one. We propose a backward error analysis framework for GMRES that simplifies the process of determining error bounds of many existing and future variants of GMRES. This framework describes modular bounds for the attainable normwise backward and forward errors of the computed solution that can be specialized for a given GMRES variant under minimal assumptions. To assess the relevance of our framework we first show that it is compatible with the previous rounding error analyses of GMRES in the sense that it delivers (almost) the same error bounds under (almost) the same conditions. Second, we explain how to use this framework to determine new error bounds for GMRES algorithms that do not have yet or have an incomplete backward error analysis, such as simpler GMRES, CGS2-GMRES, and mixed precision GMRES.
We build a finite volume scheme for the scalar conservation law $\partial _{t} u + \partial _{x} (H(x, u)) = 0$ with initial condition $u_{o} \in \mathbf{L}<^>{\infty }(\mathbb{R}, \mathbb{R})$ for a wide class of flux function $H$, convex with respect to the second variable. The main idea for the construction of the scheme is to use the theory of discontinuous flux. We prove that the resulting approximating sequence converges a.e. boundedly in $\mathopen ] 0, +\infty \mathclose [ imes \mathbb{R}$, to the entropy solution.
We consider the standard finite element method combined with the implicit Euler scheme to approximate second-order semilinear parabolic stochastic partial differential equations (SPDEs) with additive noise. The nonlinearity is of Nemytskii type, satisfies a one-sided Lipschitz condition, exhibits polynomial growth and includes irregular components. Such SPDEs serve as suitable models for various phenomena, including advection-reaction-diffusion processes in a convex polyhedral domain $\varLambda \subset \mathbb{R}<^>{d}$ ($d\in \{1,2, 3\})$. We prove strong convergence of the fully discrete scheme towards the mild solution, achieving an almost first-order temporal rate. We obtain a spatial convergence rate that, in the case $d=3$, depends on the spatial dimension and the polynomial growth order of the nonlinearity. The analysis is challenging because of the irregularities of the nonlinear drift function and the absence of a global Lipschitz condition. Numerical experiments are provided to illustrate the theoretical results.
We conduct a comprehensive global-in-time energy stability of temporally first- to third-order accurate exponential-time-differencing Runge-Kutta schemes for the Cahn-Hilliard-Ohta-Kawasaki equation modeling the microphase separation of diblock copolymer melts. The theoretical challenge lies in the energy stability analysis without assuming global Lipschitz nonlinearity or $an\ L<^>\infty $ bound on solutions. To address this, we propose a general framework for estimating a uniform bound on all stage solutions, which in turn determines the stabilization parameter required to ensure energy stability. This approach eliminates the need for traditional assumptions, as well as any restriction on the time step related to the spatial discretization or the convergence constant. Because of the absence of a high-order Sobolev norm in the energy functional, preliminary $H<^>{1}$ and $H<^>{2}$ estimates are established for each stage solution, with rough $L<^>\infty $ bounds obtained via the log-interpolation-embedding inequality in two dimensions. Consequently, by constructing an inverse operator, we establish a uniform-in-time $H<^>{2}$ estimate for the final stage solution, subject to an $\mathscr{O}(\varepsilon <^>{6} |\ln \varepsilon |<^>{-4})$ time-step constraint. This leads to a uniform-in-time $L<^>\infty $ bound for all stage solutions through the log-interpolation-embedding inequality. Consequently, we achieve an $\mathscr{O}(\varepsilon <^>{-2} |\ln \varepsilon |<^>{2})$ stabilization parameter, which guarantees global-in-time energy stability. The proposed framework is quite general and can be extended to other single-step schemes, such as implicit-explicit Runge-Kutta and exponential Runge-Kutta schemes, as well as to phase field models, including the classic Allen-Cahn and Cahn-Hilliard-type equations.
The Landau-Lifshitz-Bloch equation (LLBE) describes the evolution of the magnetic spin field in ferromagnets at high temperatures. In this paper we study the numerical approximation of the LLBE in the regime above the Curie temperature on bounded polytopal domains in $\mathbb{R}<^>{d}$, with $d\in \{1,2,3\}$, allowing for nonconvex domains when $d=2$. We first establish the existence and uniqueness of strong solutions to the LLBE and propose a linear, fully discrete, conforming finite element scheme for its approximation. While this scheme is shown to converge, the obtained rate is suboptimal. To address this shortcoming we introduce a viscous (pseudo-parabolic) regularization of the LLBE, which we call the $\epsilon $-LLBE. For this regularized problem we prove the unique existence of strong solutions and establish a rate of convergence of the solution $\boldsymbol{u}<^>\epsilon $ of the $\epsilon $-LLBE to the solution $\boldsymbol{u}$ of the LLBE as $\epsilon o 0<^>{+}$. Furthermore, we propose a linear, fully discrete, conforming finite element scheme to approximate the solution of the $\epsilon $-LLBE. Given sufficiently smooth initial data error analysis is performed to show stability and uniform-in-time convergence of the scheme. Finally, several numerical simulations are presented to corroborate our theoretical results.
Numerical approximations to the invariant measures of a class of stochastic differential equations (SDEs) with periodic coefficients are studied in this work. Compared with those existing fruitful results on the numerical studies on the invariant measures of autonomous SDEs, to our best knowledge, this paper is the first one to devote to the case of nonautonomous SDEs that are extensively recognized to characterize various real-world problems in finance and biology. Technical challenges including time-inhomogeneity and periodicity make this paper a challenging and nontrivial work. The existence and uniqueness of the invariant measure of the numerical solution and its convergence to the underlying one in the Wasserstein distance are proved in this paper. Moreover, we obtain the exponential decay rate of the numerical solution to its invariant measure and the polynomial convergence rate of the numerical invariant measure to the underlying one. Numerical simulations are provided to illustrate those theoretical results.
This paper presents a theory of nonconforming finite element exterior calculus based on a unified family of nonconforming finite element spaces for $H\varLambda <^>{k}$ in ${\mathbb{R}<^>{\boldsymbol{{n}}}}$ ($0\leqslant k\leqslant n$, $n\geqslant 1$), which are constructed in this paper by a novel approach that seeks to mimic the dual connections between adjoint operators. The family each employs piecewise Whitney forms as shape functions, including the lowest-degree Crouzeix-Raviart element space for $H\varLambda <^>{0}$, and optimal approximations and uniform discrete Poincar & eacute; inequalities are presented. Further, with these newly constructed finite element spaces, discrete de Rham complexes with commutative diagrams, and the discrete Helmholtz decomposition and Hodge decomposition for piecewise constant spaces are established, based on which the Poincar & eacute;-Leftschetz duality can be reconstructed discretely as an equality. The consequent framework of nonconforming finite element exterior calculus is naturally connected to the classical conforming one, but significantly different. Notably, all discrete operators involved are local, namely acting cell by cell separately. The newly constructed finite element spaces do not fit Ciarlet's finite element definition, though they admit locally supported basis functions, each spanning at most two adjacent cells, which makes the computation of the local stiffness matrices and the assembling of the global stiffness matrices implementable by following the standard procedure. Some numerical experiments are given to show the implementability and the performance of the new kind of spaces. The cooperation of conforming and nonconforming finite element spaces leads to new discretization schemes of the Hodge-Laplace problem.
This paper develops and analyses the mixed finite element methods (MFEMs) for It & ocirc;-type stochastic partial differential equations (SPDEs) driven by gradient-dependent multiplicative noise. Solutions to such SPDEs may lose regularity rapidly or potentially blow up in finite time due to gradient-dependent stochastic effects. By treating the gradient of the solution as an independent variable in MFEM, we explicitly track the influence of the gradient $\nabla u$ on the proposed schemes by setting the specified coefficient parameters relationship in the proof. We rigorously prove that the proposed semidiscrete and fully discrete schemes can theoretically achieve optimal strong convergence rates in space and time. Numerical tests are also presented to validate the theoretical results and demonstrate the effectiveness of our proposed method for SPDEs with gradient-dependent stochasticity.
