Abstract This special issue is dedicated to the memory of Raytcho Lazarov, who passed away on November 30, 2024. The articles collected here reflect some of his principal research interests, particularly in numerical methods for partial differential equations and the development of algorithms for their efficient and reliable application. We briefly highlight his major scientific contributions and provide an overview of the results presented in this volume. The contributions are authored by his former students, colleagues, coauthors, and long-term collaborators – researchers whose work has been deeply influenced by Lazarov’s ideas and achievements. Over more than half a century of remarkably productive scholarship, he made fundamental contributions across a broad range of important areas in computational and applied mathematics.
In this paper, we study a reduced-order model based on proper orthogonal decomposition (POD) for a two-dimensional semilinear parabolic equation. The weak Galerkin finite element method (WG-FEM) is used for spatial discretization and the Crank-Nicolson (CN) scheme is employed for temporal discretization. We introduce the POD formulation and establish the optimal error estimates for both WG-FEM and WG-POD solutions, achieving second-order accuracy in time and O(h(k+1)) accuracy in space in the L-2 norm, and the numerical results are consistent with the theoretical conclusions. This indicates that the method can significantly reduce the degrees of freedom and save CPU time while maintaining high accuracy.
In this paper, we design two classes of numerical schemes for the time-fractional Allen-Cahn equation based on the discrete gradient method and the L1-type formula of the Riemann-Liouville fractional derivative. The suggested schemes are as accurate and efficient as a recent work [H. Liao, T. Tang and T. Zhou, An energy stable and maximum bound preserving scheme with variable time steps for time fractional Allen-Cahn equation, SIAM J. Sci. Comput. 43 2021, 5, A3503-A3526] for the time-fractional Allen-Cahn equation with a double-well potential function, but can address the time-fractional Allen-Cahn equation with a general potential function. Thanks to the orthogonal convolution technique on nonuniform time mesh: (i) the first type of schemes unconditionally preserves the discrete energy dissipation law; (ii) the second type of schemes is constructed by introducing the fixed-point iteration and stabilized factor method for the first scheme, and proven to be maximum bound principle preserving. Finally, extensive numerical comparisons are reported to verify the efficiency of the proposed schemes and the correctness of the theoretical results.
We analyze discrete approximations of the second fundamental form on graphs of functions that are piecewise affine on irregular meshes. Being related with the Morley finite element, the approximation in this precise form has first been suggested in [E. Grinspun, Y. Gingold, J. Reisman and D. Zorin, Computing discrete shape operators on general meshes, Comput. Graphics Forum 25 (2006), no. 3, 547-556]. We show how to use this framework to variationally approximate functionals of the form E-0 (M) :=. integral(M) f(x, n(M)(x), Dn(M)(x)) dH(2)(x), where n(M) denotes the normal of the surface M. Here the integrand f is not necessarily quadratic. This corresponds to nonlinear Euler-Lagrange equations. Our approximation is rigorously formulated in the framework of G-convergence: We combine an ansatz-free asymptotic lower bound for any uniform approximation and a recovery sequence consisting of any regular triangulation of the limit sequence and an almost optimal choice of edge director. We give numerical examples showing the efficiency and accuracy of the algorithm in nonlinear problems.
In this work, we present a multicontinuum homogenization approach for a two-phase model with linear transport, governed by coupled flow and transport equations with high-contrast coefficients. The method introduces multiple macroscopic variables to represent distinct physical characteristics of the solution, such as local averages over subregions. To capture fine-scale heterogeneities, multiscale basis functions are constructed by solving local constrained energy minimization problems. The resulting reduced-order model is defined on a coarse scale and involves the macroscopic variables and homogenized coefficients. It retains the essential features of the original fine-scale system while reducing computational costs. Numerical examples are provided to validate the proposed model. The results demonstrate that the solutions obtained from the proposed model accurately approximate those from the original model.
Roos and Stynes [Some open questions in the numerical analysis of singularly perturbed differential equations, Comput. Methods Appl. Math. 15 2015, 4, 531-550] posed several open problems, one of which is obtaining uniform convergence of higher-order finite element methods on Bakhvalov-type meshes. In this article, we propose a finite element analysis of any order on an eXp-Bakhvalov mesh to solve a two-dimensional singularly perturbed boundary value problem, whose solution exhibits exponential layers. A careful selection of the interpolation operator, considering the characteristics of the layers, allows the finite element method to achieve optimal-order convergence with respect to the singular perturbation parameter. Numerical results are presented to support the theoretical findings.
We develop multipoint stress mixed finite element methods for linear elasticity with weakly enforced stress symmetry on distorted quadrilateral grids, which can be reduced to positive definite cell-centered systems. The methods utilize the lowest-order Brezzi-Douglas-Marini finite element spaces for the stress and employ vertex quadrature rules to localize the interaction of degrees of freedom. This approach allows for local stress elimination around each vertex. We introduce two methods. The first method uses a piecewise constant rotation, resulting in a cell-centered system for the displacement and the rotation. The second method employs a continuous piecewise bilinear rotation, enabling further elimination of the rotation and resulting in a cell-centered system for the displacement only. The methods utilize a non-symmetric vertex quadrature rule for the stress bilinear form and both non-symmetric and symmetric vertex quadrature rules for the asymmetry bilinear forms. Stability and error analysis are performed for both methods. First-order convergence is established for all variables in the L 2 {L<^>{2}} -norm. Numerical results are presented that verify the theoretical results.
We present a method for the numerical approximation of an optimal control problem constrained by a full-space transmission problem. The problem can be controlled by the interior source term, the jump over the interface, the jump of the flux over the interface, or any combination thereof. We complete the first-order optimality condition by the Johnson-N & eacute;d & eacute;lec coupling, and a Lagrange multiplier which corresponds to the Bielak-MacCamy coupling. Our final formulation fulfills the LBB conditions on the continuous as well as discrete level, without restrictions. We develop a reliable a posteriori error estimator, propose an adaptive algorithm, and show its rate optimality. Numerical experiments confirm our theoretical findings.
This article presents a MATLAB software package for solving a three-dimensional (3D) second-order elliptic problem with mixed boundary conditions using the primal hybrid finite element method (FEM). First, we introduce a novel fast 3D uniform finite element mesh refinement technique implemented in MATLAB and establish the Nodes-to-Edge and Faces-to-Tetrahedron connectivity through an efficient and systematic approach. We then describe an efficient MATLAB assembly procedure for the 3D lowest-order primal hybrid finite element matrices. Furthermore, we develop a vectorized Schur complement solver, where the computational improvement over MATLAB's default direct solver (mldivide) arises from reducing the original block system to a substantially smaller Schur complement system. The run-time performance of the software is demonstrated through numerical experiments.
In this work, we present a semi-implicit structure-preserving nonstaggered central scheme (SP-NCS) that is both well-balanced and positivity-preserving, designed to solve the two-dimensional shallow water Exner equations with time-dependent bottom topography. The SP-NCS can be viewed as a two-dimensional extension of the numerical method presented in [D. Li and J. Dong, A robust hybrid unstaggered central and Godunov-type scheme for Saint-Venant-Exner equations with wet/dry fronts, Comput. & Fluids 235 2022, Article ID 105284], which is uncoupled for the shallow water Exner equations. The SP-NCS is also inspired by the nonstaggered central scheme presented in [J. Dong and X. Qian, Structure-preserving nonstaggered central difference schemes at wet-dry fronts for the shallow water equations, Commun. Appl. Math. Comput. 8 2026, 1, 366-410], but their method is significantly adapted to preserve the still-water steady state when the computational domain contains wetting and drying transitions. Retaining the stationary solution in the two-dimensional case, especially when the domain include wetting and drying transitions, presents a significant challenge. This is because the backward step and the discretization of the source term in the corrector step become more complex when dealing with such fronts. A key innovation is the introduction of a structure-preserving parameter that helps retain the stationary solution even in the presence of wet-dry fronts. This task is nontrivial due to the robust requirements involved in ensuring the stability and accuracy of the solution. In particular, the SP-NCS demonstrates robustness in handling the relatively strong interactions between the two coupled models even with a large time step. This is particularly relevant when the rate of change in the bed topography is much slower than the speed of the water surface waves. We rigorously prove the positivity-preserving and well-balanced properties of the SP-NCS. Finally, we conduct several numerical experiments to validate the theoretical results and demonstrate the effectiveness of the proposed method in preserving the stationary solution and in scenarios involving wet-dry fronts, especially for the relatively strong interactions between the two coupled models.
A linear BDF2 numerical scheme is proposed to solve the Boussinesq system. By using the exponential scalar auxiliary variable (E-SAV) approach, we explicitly deal with the nonlinear terms of the Boussinesq system, and decouple the velocity and temperature in the numerical simulation. These numerical scheme is unconditionally stable. We give rigorous error analysis for the velocity and temperature. Numerical experiment is performed to verify the proposed numerical scheme.
We introduce a framework for repurposing error estimators for source problems to compute an estimator for the gap between eigenspaces and their discretizations. Of interest are eigenspaces of finite clusters of eigenvalues of unbounded nonselfadjoint linear operators with compact resolvent. Eigenspaces and eigenvalues of rational functions of such operators are studied as a first step. Under an assumption of convergence of resolvent approximations in the operator norm and an assumption on global reliability of source problem error estimators, we show that the gap in eigenspace approximations can be bounded by a globally reliable and computable error estimator. Also included are applications of the theoretical framework to first-order system least squares (FOSLS) discretizations and discontinuous Petrov-Galerkin (DPG) discretizations, both yielding new estimators for the error gap. Numerical experiments with a selfadjoint model problem and with a leaky nonselfadjoint waveguide eigenproblem show that adaptive algorithms using the new estimators give refinement patterns that target the cluster as a whole instead of individual eigenfunctions.
This paper delves into a class of nonlinear nonlocal partial differential equations, characterized by a gradient flow structure, as previously outlined in the literature. Finite volume methods are employed to investigate the numerical solutions of this model in two-dimensional settings. It reveals that the semi-discrete numerical scheme satisfies the entropy dissipation and the fully discrete numerical scheme preserves the positivity through a meticulous combination of numerical scheme definitions, upwind numerical fluxes and sophisticated estimation techniques. Furthermore, a novel contribution is made by demonstrating that the semi-discrete numerical scheme preserves the positivity and the fully discrete numerical scheme satisfies the entropy dissipation. These fundamental properties, pivotal in guaranteeing the convergence of numerical solutions towards steady states, are rigorously proven utilizing the invariant region method, complemented by detailed norm estimations and rigorous analytical techniques. Two-dimensional numerical experiments validate the entropy dissipation and the long-time convergence of numerical solutions generated by the scheme. The numerical analysis methodology and computational findings in this work offer valuable insights for research on models with gradient flow structures.
In this paper we formulate and analyze adaptive (space-time) least-squares finite element methods for the solution of convection-diffusion equations. The convective derivative v . del u is considered as part of the total time derivative d/dt u = partial derivative(t)u + v . del u, and therefore we can use a rather standard stability and error analysis for related space-time finite element methods. For stationary problems we restrict the ansatz space H-0(1)(Omega) such that the convective derivative is considered as an element of the dual H-1(Omega) of the test space H-0(1)(Omega), which also allows unbounded velocities v. While the discrete finite element schemes are always unique solvable, the numerical solutions may suffer from a bad approximation property of the finite element space when considering convection dominated problems, i.e., small diffusion coefficients. Instead of adding suitable stabilization terms, we aim to resolve the solutions by using adaptive (space-time) finite element methods. For this we introduce a least-squares approach where the discrete adjoint defines local a posteriori error indicators to drive an adaptive scheme. Numerical examples illustrate the theoretical considerations.
Computational technologies for the approximate solution of multidimensional boundary value problems often rely on irregular computational meshes and finite-volume approximations. In this framework, the discrete problem represents the corresponding conservation law for control volumes associated with the nodes of the mesh. This approach is most naturally and consistently implemented using Delaunay triangulations together with Voronoi diagrams as control volumes. In this paper, we employ meshes with nodes located both at the vertices of Delaunay triangulations and at the generators of Voronoi partitions. The cells of the merged Voronoi-Delaunay mesh are orthodiagonal quadrilaterals. On such meshes, scalar and vector functions, as well as invariant gradient and divergence operators of vector calculus, can be conveniently approximated. We illustrate the capabilities of this approach by solving a steady-state diffusion-reaction problem in an anisotropic medium.
We present and discuss a generalization of the popular MINI mixed finite element for the 2D Stokes equation by means of conforming virtual elements on polygonal meshes. We prove optimal error estimates for both velocity and pressure. Theoretical results are confirmed by several numerical tests performed with different choices of polynomial accuracy and meshes.
The aim of this paper is to propose an efficient spectral-Galerkin method for the numerical approximation of the time-space fractional diffusion equation in an unbounded domain by using the fractional-order generalized Jacobi functions and the mapped Chebyshev functions, and theoretically prove the high convergence rate. The reason for using the fractional-order generalized Jacobi functions is to approximate the solution with singularity at the initial time. We prove that the proposed approximation scheme has a spectral convergence rate when the solution of a given problem satisfies a particular condition. We also establish the stability of the method.
We study the popular modularity matrix and respective functional used in connection with graph clustering and derive some properties useful when performing vertex aggregation of the associated graph. These properties are employed in the derivation of a multilevel parallel pairwise aggregation algorithm. Comparative performance results of the studied algorithm applied to graph clustering tested against the popular Louvain algorithm are presented. Some illustrative examples show that the resulting aggregates if used in an adaptive algebraic multigrid (AMG) are able to follow strong direction of anisotropy in finite element problems.
Large linear systems are ubiquitous in modern computational science and engineering. Their efficient solution is frequently challenging, especially in conjunction with problems described by parametric PDEs, where many such systems have to be solved. The primary approach for solving large linear systems is the use of Krylov subspace iterative methods with well-designed preconditioners. Recently, graph neural networks (GNNs) have been used in conjunction with solving parametric PDEs, and have been shown to be promising tool for designing preconditioners. Nevertheless, the GNN-based preconditioners reported in the literature have worse effect on system's spectrum than analogous preconditioners from classical linear algebra. Here we employ well-established preconditioners from linear algebra as starting point for training GNN to obtain preconditioners that yield a more substantial reduction in the condition number of the systems when compared to classical preconditioners. Numerical experiments demonstrate the efficiency of our approach in comparison to classical and known neural network-based methods for parametric PDEs. In addition, a heuristic justification for the loss function employed in this study is provided, and we demonstrate that preconditioners obtained by learning with this loss function reduce the condition number in a way that is more desirable for the conjugate gradient iterative method.
We consider quasi-polynomial spaces of differential forms defined as weighted (with a positive weight) spaces of differential forms with polynomial coefficients. We show that the unisolvent set of functionals for such spaces on a simplex in any spatial dimension is the same as the set of such functionals used for the polynomial spaces. The analysis in the quasi-polynomial spaces, however, is not standard and requires a novel approach. We are able to prove our results without the use of Stokes' Theorem, which is the standard tool in showing the unisolvence of functionals in polynomial spaces of differential forms. These new results provide tools for studying exponentially-fitted discretizations stable for general convection-diffusion problems in Hilbert differential complexes.
