Preconditioned iterative method for nonsymmetric saddle point linear systems

Abstract

In this paper, a new preconditioned iterative method is presented to solve a class of nonsymmetric nonsingular or singular saddle point problems. The implementation of the proposed preconditioned Krylov subspace method avoids solving inverse of Schur complement and only needs to solve one linear sub-system at each step, which implies that it may save considerable costs. Theoretical convergence analysis, including the bounds of eigenvalues and eigenvectors, the degree of the minimal polynomial of the preconditioned matrix, are discussed in details. Moreover, a novel algebraic estimation technique for finding a practical iteration parameter is presented, which is very effective and practical even for large scale problems. At last, some numerical examples are carried, showing that the theoretical results are valid and convincing.

Introduction

In this paper, we focus ourselves on the solution of the following saddle point linear equations

$$\mathcal{A}x=\left(\begin{array}{cc}\hfill A\hfill & \hfill B\hfill \\ \hfill -B^T\hfill & \hfill 0\hfill \end{array}\right)\left(\begin{array}{c}\hfill u\hfill \\ \hfill p\hfill \end{array}\right)=\left(\begin{array}{c}\hfill f\hfill \\ \hfill -g\hfill \end{array}\right)=b,$$

where $A\in {\mathbb{R}}^{n\times n}$ is a positive real matrix, i.e., its symmetric part is positive definite, $B\in {\mathbb{R}}^{n\times m}$ is a rectangular matrix with $n\gg m$, $u,f\in {\mathbb{R}}^n$ and $p,g\in {\mathbb{R}}^m$. Frequently, B is full rank, which means $\mathcal{A}$ is nonsingular. Whereas B may be rank deficient, in this case the coefficient matrix $\mathcal{A}$ is a singular matrix. In recent years, linear system (1.1) has got great attention owing to its widely applications in science and engineering problems [1], [2], [3], [4], [5]. To be more specific, we consider the following incompressible Navier-Stokes equations

$$\{\begin{array}{c}-\nu \mathrm{\Delta }\mathbf{\text{u}}+\mathbf{\text{u}}\cdot \mathrm{\nabla }\mathbf{\text{u}}+\mathrm{\nabla }p=\mathbf{\text{f}},\text{in}\mathrm{\Omega },\hfill \\ \mathrm{\nabla }\cdot \mathbf{\text{u}}=0,\text{in}\mathrm{\Omega },\hfill \\ \mathbf{\text{u}}=\mathbf{\text{g}},\text{on}\partial \mathrm{\Omega }\hfill \end{array}$$

in an open bounded domain $\mathrm{\Omega }\subset {\mathbb{R}}^2$ (or ${\mathbb{R}}^3$) with boundary ∂Ω, $\nu >0$ is the viscosity parameter. Linearization of the Navier-Stokes system (1.2) by Picard fixed-point iteration methods will involve in solving a sequence of Oseen problems of the form

$$\{\begin{array}{c}-\nu \mathrm{\Delta }\mathbf{\text{u}}+\mathbf{\text{w}}\cdot \mathrm{\nabla }\mathbf{\text{u}}+\mathrm{\nabla }p=\mathbf{\text{f}},\text{in}\mathrm{\Omega },\hfill \\ \mathrm{\nabla }\cdot \mathbf{\text{u}}=0,\text{in}\mathrm{\Omega },\hfill \\ \mathbf{\text{u}}=\mathbf{\text{g}},\text{on}\partial \mathrm{\Omega },\hfill \end{array}$$

where the divergence free field w is the velocity field obtained from the previous Picard iteration step. We use finite element methods to discretize the Oseen equations (1.3). In order to guarantee a unique solution, we let the finite element discretization satisfy the Ladyzhenskaya-Babuša-Brezzi (LBB) condition. If we use the $Q2-Q1$ element or $Q2-P1$ element [5], then we will obtain the following large sparse systems of the form

$${\mathcal{A}}_{+}x=\left(\begin{array}{cc}\hfill A\hfill & \hfill B\hfill \\ \hfill B^T\hfill & \hfill 0\hfill \end{array}\right)\left(\begin{array}{c}\hfill u\hfill \\ \hfill p\hfill \end{array}\right)=\left(\begin{array}{c}\hfill f\hfill \\ \hfill g\hfill \end{array}\right)\equiv b_{+},$$

where A is a convection-diffusion-type matrix, B and $B^T$ are the discrete divergence and gradient matrices. u and p represent discrete velocity and pressure, respectively, and f and g contain forcing and boundary terms. By rearranging (1.4), it is easy to get (1.1). As shown in [5], if a finite element does not satisfy the LBB condition, we need to apply stabilization techniques to obtain a stable discretization scheme, which results in solving a generalized saddle point linear system of the following form

$${\mathcal{A}}_{+}x=\left(\begin{array}{cc}\hfill A\hfill & \hfill B\hfill \\ \hfill B^T\hfill & \hfill -C\hfill \end{array}\right)\left(\begin{array}{c}\hfill u\hfill \\ \hfill p\hfill \end{array}\right)=\left(\begin{array}{c}\hfill f\hfill \\ \hfill g\hfill \end{array}\right)\equiv {\tilde{b}}_{+}.$$

However, the solution of linear system (1.5) is not within our consideration. Thus, in our experiments we have limited ourselves to finite elements that satisfy the BB condition.

Efficient solution of system (1.1) necessitates rapidly convergent iterative methods. Krylov subspace methods with a good preconditioner are considered for such large linear systems since they are cheap to implement and are able to exploit the sparsity of the coefficient matrix. It is well known that spectral distribution of the preconditioned matrix with a clustering of most of the eigenvalues relates closely to favorable convergence rate of Krylov subspace iterative methods, and the application of the preconditioner within a Krylov subspace method involves repeated solution of the linear system with the preconditioner as the coefficient matrix. Therefore, good preconditioner should be easily to implement [4].

In recent years, many effective preconditioners for Krylov subspace methods have been proposed. For examples, block preconditioners [2], [6], [7], constrained block preconditioners [8], [9], HSS preconditioners and its variants [10], [11], [12], [13], [14], [15], [16], [17], [18], DS(dimensional split)-like preconditioners [19], [20], other splitting preconditioners [21], [22], [23], [24], [25], [26], [27], [28], [29], [30], [31] and so on. For a survey of preconditioners for the saddle point problems, see [4] and the references therein.

Because of the advantage of the relaxation technique used in matrix splitting preconditioners, many existing preconditioners have been improved [14], [32], [33], [34], [35], [36], [37]. In fact, the relaxation technique was originally studied in [32] for the dimensional split preconditioner. Then this idea was extended to the HSS preconditioner [34], [35], the DPSS preconditioner [14], [33] and so on. For example, by taking use of the relaxation technique, Cao et al. proposed a simplified HSS (SHSS) preconditioner and proved that the SHSS preconditioner ${\mathcal{P}}_{\text{SHSS}}$ is a better approximation to the matrix $\mathcal{A}$ and is easier to implement than the HSS preconditioner ${\mathcal{P}}_{\text{HSS}}$ [34].

In this paper, motivated by relaxation technique and the preconditioner introduced in [27], we propose a new preconditioner for solving nonsingular or singular nonsymmetric saddle point problems. Convergence and semi-convergence for the corresponding iteration method are studied. The implementation of the preconditioned method is given, which shows that the proposed preconditioned GMRES method may save considerable costs. The eigen-properties of the preconditioned matrix are discussed and the dimension of the Krylov subspace methods for the preconditioned iteration method is obtained. Furthermore, a novel algebraic estimation technique for finding a practical iteration parameter is presented. Numerical experiments from Navier-Stokes equation show that the choice strategy of the practical iteration parameter is feasible, leading to the preconditioner' performance achieving faster convergence in terms of CPU time and iteration steps and fairly insensitive to the iteration parameter.

The remainder of the paper is organized as follows. In Section 2, new preconditioner and its implementation are given. In Section 3, for nonsingular saddle point problems, the convergence properties of the corresponding iterative method are studied, spectral properties of the preconditioned matrix are also discussed. Besides, selection strategy for the iteration parameter is studied. In Section 4, the new preconditioner is extended to solve singular linear system (1.1), semi-convergence of the corresponding iteration method is discussed. Numerical results are given in Section 5 to show the effectiveness of the new preconditioner. Finally, in Section 6, we draw some conclusions to end the paper.

Throughout this paper, $\sigma (\cdot ),\rho (\cdot )$, null(⋅), rank(⋅) and index(⋅) denote spectral set, spectral radius, null space, rank and index of a matrix, respectively.