Symmetric Positive Definite Systems Research Articles

Two linear energy-preserving compact finite difference schemes for coupled nonlinear wave equations

In this paper, we propose and analyze two highly efficient compact finite difference schemes for coupled nonlinear wave equations containing coupled sine-Gordon equations and coupled Klein-Gordon equations. To construct energy-preserving, high-order accurate and linear numerical methods, we first utilize the scalar auxiliary variable (SAV) approach and introduce three auxiliary functions to rewrite the original problem as a new equivalent system. Then we make use of the compact finite difference method and the Crank-Nicolson method to propose an efficient fully-discrete scheme (SAV-CFD-CN). The modified energy conservation and the convergence of the SAV-CFD-CN scheme are proved in detail, which has fourth-order convergence in space and second-order convergence in time. In order to preserve the discrete energy of original system, we further combine Lagrange multiplier approach, compact finite difference method and the Crank-Nicolson method to propose the second fully-discrete scheme (LM-CFD-CN). The proposed two schemes are high-order accurate, linear and highly efficient, only four symmetric positive definite systems with constant coefficients are required to be solved at each time level. Numerical experiments for the coupled nonlinear wave equations are given to confirm theoretical findings.

Avoiding breakdown in incomplete factorizations in low precision arithmetic

The emergence of low precision floating-point arithmetic in computer hardware has led to a resurgence of interest in the use of mixed precision numerical linear algebra. For linear systems of equations, there has been renewed enthusiasm for mixed precision variants of iterative refinement. We consider the iterative solution of large sparse systems using incomplete factorization preconditioners. The focus is on the robust computation of such preconditioners in half precision arithmetic and employing them to solve symmetric positive definite systems to higher precision accuracy; however, the proposed ideas can be applied more generally. Even for well-conditioned problems, incomplete factorizations can break down when small entries occur on the diagonal during the factorization. When using half precision arithmetic, overflows are an additional possible source of breakdown. We examine how breakdowns can be avoided and we implement our strategies within new half precision Fortran sparse incomplete Cholesky factorization software. Results are reported for a range of problems from practical applications. These demonstrate that, even for highly ill-conditioned problems, half precision preconditioners can potentially replace double precision preconditioners, although unsurprisingly this may be at the cost of additional iterations of a Krylov solver.

A scalable physically consistent particle method for high-viscous incompressible flows

AbstractA scalable matrix solver was developed for the moving particle hydrodynamics for incompressible flows (MPH-I) method. Since the MPH-I method can calculate both incompressible and highly viscous flows while ensuring stability through physical consistency, a wide range of industrial applications is expected. However, in its implicit calculation, both the pressure and velocity must be solved simultaneously via a linear equation with a nondefinite symmetric coefficient matrix. In this study, this nondefinite linear system was converted into a symmetric positive definite (SPD) system where only the velocity is unknown. This conversion enabled us to solve the system with well-known solvers such as the conjugated gradient (CG) and conjugated residual (CR) methods. For scalability, bucket-based multigrid preconditioned CG and CR solvers were developed for the SPD system. To handle multidimensionality during preconditioning, an extended Jacobi smoother that is even applicable in a nondiagonally dominant matrix system was proposed. The numerical efficiency was confirmed via a simple high-viscosity incompressible dam break calculation, and the scalability within the presented case was confirmed. In addition, the performance under shared memory parallel computations was studied.

Optimal polynomial smoothers for multigrid V‐cycles

AbstractThe idea of using polynomial methods to improve simple smoother iterations within a multigrid method for a symmetric positive definite system is revisited. A two‐level bound going back to Hackbusch is optimized by a very simple iteration, a close cousin of the Chebyshev semi‐iterative method, but based on the Chebyshev polynomials of the fourth instead of first kind. A full V‐cycle bound for general polynomial smoothers is derived using the V‐cycle theory of McCormick. The fourth‐kind Chebyshev iteration is quasi‐optimal for the V‐cycle bound. The optimal polynomials for the V‐cycle bound can be found numerically, achieving about an 18% lower error contraction factor bound than the fourth‐kind Chebyshev iteration, asymptotically as the number of smoothing steps . Implementation of the optimized iteration is discussed, and the performance of the polynomial smoothers is illustrated with numerical examples.

Pixel-based boundary element method for computing effective thermal conductivity of heterogeneous materials

In the present article, we propose an efficient 2D Pixel-based Boundary Element formulation to compute the effective thermal conductivity of heterogeneous materials. In the proposed formulation, each pixel of the digital image is represented by a subdomain with four boundary elements. We develop the formulation within the context of the Symmetric Galerkin Boundary Element Method and, apart from the Boundary Element Method itself, the other main ingredients are the Domain Decomposition technique, to treat each pixel as a subdomain, and the Element-by-Element technique used in conjunction with the Preconditioned Conjugate Gradient method to solve the resulting symmetric positive definite large linear system. In addition, Periodic Boundary Conditions are employed. We validate our formulation and implementation with the aid of analytical results for the case of a plate with periodic circular inclusions. Further, we illustrate the application of the proposed formulation by computing the effective thermal conductivity of a plate with periodic square inclusions, a woven-like material, and a cast iron sample, and compare the obtained results with those obtained via finite elements. Our results indicate that the Pixel-based Boundary Element formulation performs, at least, as good as the Pixel-based Finite Element formulation.

Analysis of (shifted) piecewise quadratic polynomial collocation for nonlocal diffusion model

The piecewise quadratic polynomial collocation is used to approximate the nonlocal model, which generally leads to a nonsymmetric indefinite system (Chen et al. (2021) [5]). In this case, the discrete maximum principle is not satisfied, which might be trickier for the stability analysis of the high-order numerical schemes (D'Elia et al. (2020) [10]; Leng et al. (2021) [26]). Here, we present a modified (shifted-symmetric) piecewise quadratic polynomial collocation for solving the linear nonlocal diffusion model, which leads to a symmetric positive definite system and satisfies the discrete maximum principle. Using Faulhaber's formula and Riemann zeta function, the perturbation error for symmetric positive definite system and nonsymmetric indefinite system are given. Then rigorous convergence analysis for the nonlocal models are provided under the general horizon parameter δ=O(hβ), with β≥0. More concretely, the global error is O(hmin⁡{2,1+β}) if δ is not set as a grid point, while it recovers O(hmax⁡{2,4−2β}) when δ is set as a grid point. We also prove that the shifted-symmetric scheme is asymptotically compatible, which has the global error O(hmin⁡{2,2β}) as δ,h→0. The numerical experiments (including two-dimensional case) are performed to verify the convergence.

A formal construction of a divergence-free basis in the nonconforming virtual element method for the Stokes problem

We develop a formal construction of a pointwise divergence-free basis in the nonconforming virtual element method of arbitrary order for the Stokes problem introduced in Zhao et al. (SIAM J. Numer. Anal. 57(6):2730–2759, 2019). The proposed construction can be seen as a generalization of the divergence-free basis in Crouzeix-Raviart finite element space (Brenner, Math. Comp. 55(192):411–437, 1990; Thomasset, 1981) to the virtual element space. Using the divergence-free basis obtained from our construction, we can eliminate the pressure variable from the mixed system and obtain a symmetric positive definite system. Several numerical tests are presented to confirm the efficiency and the accuracy of our construction.

Efficient implementation of mixed finite element methods for parabolic problems

It is well-known that conventional mixed finite element method (FEM) results in a large saddle point system or non-symmetric system. A splitting mixed FEM is introduced by Yang (Numer. Methods Partial Differential Equations 17: 229–249, 2001), which computes the flux variable by solving a symmetric positive definite system. Then, the original unknown is solved based on the computed flux. In this paper, we prove that the splitting mixed FEM proposed by Yang is equivalent to the conventional mixed FEM provided that the initial data is carefully taken. Moreover, we show that the equivalence also holds for more complicated models, such as splitting mixed FEM for the compressible flow in porous media. Therefore, all available theoretical results of conventional mixed FEMs hold for the splitting ones. Numerical comparisons are provided for the conventional mixed FEM and the splitting mixed FEM. Efficiency of the splitting mixed method is clearly demonstrated. Corresponding results can be extended to many other mixed FEMs, such as Brezzi–Douglas–Marini mixed FEM.

Continuation Model Predictive Control Based on a Regularized and Smoothed Fischer-Burmeister Method

In this paper, the regularized and smoothed Fischer-Burmeister (FBRS) method is introduced in model predictive control based on the continuation method. By introducing the FBRS method, an optimal control problem formed as a convex problem was regularized and smoothed to satisfy regularity of an equation derived by using a continuation method. In addition, the equation was reduced to a general linear system and a symmetric positive definite system to minimize the computational burden. Furthermore, we discussed the effect of the regularizing and the smoothing parameters on condition number related to numerical accuracy in computation. A guiding principle for the parameters was obtained via the discussion. Finally, an example for a constrained multi-input multi-output system, control design of an engine airpath system was demonstrated to confirm the reasonability of the method.

Least‐squares formulation to solve non‐Newtonian fluid flow and application of data assimilation in 2D

AbstractIn this contribution, we present an approach to model steady flow of incompressible non‐Newtonian fluids with data assimilation into the numerical solution based on the least‐squares finite element method (LSFEM). The assimilation of data (e.g. experimental or analytical) into numerical simulations offers promising possibilities when examining complex problems. Potential applications include the enhancement of numerical models using measured data or the completion of experimental data using numerical methods, e.g. the determination of non‐measured quantities such as pressure. In particular for the field of fluid mechanics, the implemented LSFEM has some theoretical advantages compared to the well‐known (mixed) Galerkin finite element method, since it is not restricted to the LBB‐condition. Additionally, it results in a minimization problem with symmetric positive (semi‐)definite equation systems also for differential equations with non‐self‐adjoint operators. A further advantage in this context is that the assimilation of data can easily be performed by adding a term to the least‐squares (LS) functional such that it does not significantly increase the computational cost. The approach of data assimilation is shown by solving steady flow of a non‐Newtonian fluid through a channel with a smooth contraction.

A least-squares finite element method for steady flows across an unconfined square cylinder placed symmetrically in a plane channel

This paper develops a least-squares finite element method (LSFEM) for non-Newtonian flows of power-law fluid across an unconfined square cylinder at two inclination angles, placed symmetrically in a two-dimensional channel. We used Newton's method to linearize the equation of motion and employed the least-squares approach with systematic mesh refinements to flow past the square cylinder. The LSFEM offers a direct approximation of the extra stress tensor components, a symmetric positive definite system, and the openness of choosing finite element spaces. We proved that the least-squares approximation converges to linearized solutions of non-Newtonian fluid flows. We demonstrated that the numerical results using low-order basis functions agree with the theoretical estimates and presented the effects of different physical parameters on the physical attributes of the power-law model at a low Reynolds number, such as drag coefficients, viscosity, and velocity. These results agree with others published in the literature. Hence, the LSFEM in non-Newtonian fluid simulations can accomplish the identification of shear-thinning and shear-thickening characteristics.

A numerical scheme with divergence free H-div triangular finite element for the Stokes equations

In this short paper, we introduce an efficient finite element method for solving the Stokes equations based on a divergence free H(div) element. The strong curl of C0-Pk+1 basis functions forms a basis of the divergence-free subspace of the Pk-H(div) (i.e. BDMk) finite element space. The [Pk+1]2×2 weak gradients are employed on these non-H1 functions to obtain a conforming weak formulation. A strength of this approach is that the saddle point problem is replaced by a symmetric positive definite system with much less unknowns, which can be solved efficiently. Numerical examples are investigated to test accuracy and robustness of the method.

An efficient multigrid solver for two-dimensional spatial fractional diffusion equations with variable coefficients

Extrapolation cascadic multigrid (EXCMG) method with the conjugate gradient smoother is shown to be an efficient solver for large sparse symmetric positive definite systems resulting from linear finite element discretization of second-order elliptic boundary value problems [Pan et al. J. Comput. Phys. 344 (2017) 499–515]. In this paper, we generalize the EXCMG method to solve a class of spatial fractional diffusion equations (SFDEs) with variable coefficients. Both steady-state and time-dependent problems are considered. First of all, space-fractional derivatives defined in Riemann–Liouville sense are discretized by using the weighted average of shifted Grünwald formula, which results in a fully dense and nonsymmetric linear system for the steady-state problem, or a semi-discretized ODE system for the time-dependent problem. Then, to solve the former problem, we propose the EXCMG method with the biconjugate gradient stabilized smoother to deal with the dense and nonsymmetric linear system. Next, such technique is extended to solve the latter problem since it becomes fully discrete when the Crank-Nicolson scheme is introduced to handle the temporal derivative. Finally, several numerical examples are reported to show that the EXCMG method is an efficient solver for both steady-state and time-dependent SFDEs, and performs much better than the V-cycle multigrid method with banded-splitting smoother for time-dependent SFDEs [Lin et al. J. Comput. Phys. 336 (2017) 69–86].

Preconditioning for symmetric positive definite systems in balanced fractional diffusion equations

In this paper, we study the finite volume discretization method for balanced fractional diffusion equations where the fractional differential operators are comprised of both Riemann-Liouville and Caputo fractional derivatives. The main advantage of this approach is that new symmetric positive definite Toeplitz-like linear systems can be constructed for solving balanced fractional diffusion equations when diffusion functions are non-constant. It is different from non-symmetric Toeplitz-like linear systems usually obtained by existing numerical methods for fractional diffusion equations. The preconditioned conjugate gradient method with circulant and banded preconditioners can be applied to solve the proposed symmetric positive definite Toeplitz-like linear systems. Numerical examples, for both of one- and two- dimensional cases, are given to demonstrate the good accuracy of the finite volume discretization method and the fast convergence of the preconditioned conjugate gradient method.

A least‐squares finite element formulation to solve incompressible non‐Newtonian fluid flow

AbstractThe Modeling of non‐Newtonian fluids is an essential issue for many industrial applications, for example in the field of chemistry, bioengineering and medical science. In this contribution we present a least‐squares (LS) finite element approach to model steady flow of incompressible non‐Newtonian fluids based on the Navier‐Stokes equations. The application of the least‐squares finite element method (LSFEM) especially in the case of fluid mechanics is motivated by some theoretical advantages compared to the well‐known (mixed) Galerkin method. The LSFEM is not restricted to the LBB‐condition and results in a minimization problem with symmetric positive definite equation systems also for differential equations with non self‐adjoint operators. Furthermore, the LSFEM provides an inherent a posteriori error estimator without additional costs, which enables an efficiency enhancement through adaptive mesh refinements.In this contribution the first‐order LS formulation in terms of stresses and velocities, as introduced in [4], is extended to consider the nonlinear dependence of the viscosity on the shear rate of the fluid. Therefore, the well‐known Carreau model for generalized Newtonian fluids, see for instance [1], is investigated. For the approximation of the independent variables, velocities and stresses, we apply conforming discretizations in H1 and H(div) using standard Lagrange interpolation polynomials and vector‐valued Raviart‐Thomas interpolations functions. A flow through a square domain with exact solution is considered to validate the proposed schemes with respect to accuracy and efficiency.

Parallel gaussian elimination of symmetric positive definite band matrices for shared-memory multicore architectures

This study presents a new parallel Gaussian elimination approach for symmetric positive definite band systems. For each task, the appropriate start time and adequate processor are determined. Unnecessary dependencies between tasks are eliminated. Simultaneously, all processors perform their associated tasks with precedence constraints under consideration. Our main goal is to obtain a high degree of parallelism by balancing the load of processors and reducing the total idle and parallel execution times. The theoretical lower bounds for parallel execution time and number of processors required to execute the precedence graph at an optimal time are also computed. The validity of our investigation is confirmed by carrying out several experiments on a shared-memory multicore architecture using OpenMP. Practical results prove the efficiency of the proposed method.

Two-Level Nyström--Schur Preconditioner for Sparse Symmetric Positive Definite Matrices

Randomized methods are becoming increasingly popular in numerical linear algebra. However, few attempts have been made to use them in developing preconditioners. Our interest lies in solving large-scale sparse symmetric positive definite linear systems of equations where the system matrix is preordered to doubly bordered block diagonal form (for example, using a nested dissection ordering). We investigate the use of randomized methods to construct high quality preconditioners. In particular, we propose a new and efficient approach that employs Nystr\"om's method for computing low rank approximations to develop robust algebraic two-level preconditioners. Construction of the new preconditioners involves iteratively solving a smaller but denser symmetric positive definite Schur complement system with multiple right-hand sides. Numerical experiments on problems coming from a range of application areas demonstrate that this inner system can be solved cheaply using block conjugate gradients and that using a large convergence tolerance to limit the cost does not adversely affect the quality of the resulting Nystr\"om--Schur two-level preconditioner.

A least-squares virtual element method for second-order elliptic problems

In this paper, a least-squares virtual element method is presented for approximating the vector and scalar variables of second-order elliptic problems. The H(div)-conforming and scalar-conforming virtual elements are used to approximate the vector and scalar variables, respectively. The method allows the use of very general polygonal meshes and leads to a symmetric positive definite system. The optimal a priori error estimates are established for the vector variable in H(div) norm and the scalar variable in H1 norm. A simple a posteriori error estimator is also presented and proved to be reliable and efficient. The virtual element method handles the hanging nodes naturally, thus the local mesh post-processing to remove hanging nodes is not required. Numerical experiments are conducted to verify the accuracy of the method, and show the effectiveness and flexibility of the adaptive strategy driven by the proposed estimator and suitable mesh refinement strategy.

Spectral estimates for saddle point matrices arising in weak constraint four‐dimensional variational data assimilation

SummaryWe consider the large sparse symmetric linear systems of equations that arise in the solution of weak constraint four‐dimensional variational data assimilation, a method of high interest for numerical weather prediction. These systems can be written as saddle point systems with a 3 × 3 block structure but block eliminations can be performed to reduce them to saddle point systems with a 2 × 2 block structure, or further to symmetric positive definite systems. In this article, we analyse how sensitive the spectra of these matrices are to the number of observations of the underlying dynamical system. We also obtain bounds on the eigenvalues of the matrices. Numerical experiments are used to confirm the theoretical analysis and bounds.

Switching preconditioners using a hybrid approach for linear systems arising from interior point methods for linear programming

In general, interior point methods are successful in solving large-scale linear programming problems. Their effectiveness is determined by how fast they calculate each solution of a linear system. When solving large-scale linear systems, iterative methods are useful options since they require slightly more computational memory in each iteration and preserve the sparsity pattern of the matrix system. In this context, a common choice is the conjugate gradient method using a preconditioning strategy. Choosing a single preconditioner that fits well during all iterations of the optimization method is not an easy task, because the distribution of the eigenvalues of the system matrix may vary significantly from the first to the last iterations. In order to simplify this task, one can use different preconditioners in different iterations, in a hybrid preconditioner strategy. In the case of symmetric positive definite systems, the Controlled Cholesky Factorization achieves excellent performance for the first interior point iterations, whereas the Splitting preconditioner is very useful in the last iterations. However, since we apply a hybrid approach combining both preconditioners, it is necessary to decide when using each preconditioner. This paper addresses the critical issue of choosing between the two preconditioners of the hybrid strategy by using the condition number of the matrix system. The Ritz values obtained from the conjugate gradient method provide an approximation to the eigenvalues, which offers an estimate of the condition number. The main contribution of this research is a new heuristic to switching the preconditioners based on the estimated condition number. Numerical results for large-scale problems show that our choice to change the preconditioners adds both speed and robustness to a hybrid approach that combines the Controlled Cholesky Factorization and the Splitting preconditioner.