Programs for Programmers

Proof that Simple Colsum modified ILU Factorizations of Matrices arising in Fluid Flow Problems are not Singular

by John Appleyard B.A., Ph D, Polyhedron Software Ltd, May 2004

 

 

Incomplete LU factorizations are often used to precondition the conjugate gradient, GMRES or similar methods for iterative solution of the linear equations A.x=b. If A is divided into diagonal, lower and upper parts

 

A = D + L + U............................................................................. (1)

 

then, in the simplest case, the preconditioning matrix, B, is defined by

 

B = (γ +L) γ^-1(γ+U)...................................................................... (2)

B = γ +L + U + L γ^-1U.................................................................. (3)

 

where γ is a diagonal matrix. The rate of convergence of to a solution depends critically on how well B approximates A, and that, in turn, depends to an appreciable degree on how γ is defined. Common choices for γ are

 

γ = D – diag(L γ^-1U)..................................................................... (4)

γ = D – rowsum(L γ^-1U)................................................................ (5)

and γ = D – colsum(L γ^-1U)................................................................. (6)

 

where diag(X) is the diagonal matrix formed by ignoring all off-diagonal coefficients of X, rowsum(X) is the diagonal matrix formed by summing the coefficients in each row of X, and colsum(X) is the diagonal matrix formed by summing the coefficients in each column of X.

 

Approximation (4) is known to be stable, but relatively slow. The “modified ILU” method, represented by approximation (5) is usually faster, but sometimes prone to numerical problems. The third choice, (6) was proposed by Appleyard, Cheshire and Pollard, who noted that this choice results in the elimination of material balance errors in many fluid flow problems. In this note, I will show that, in those problems, preconditioners formed using approximation (5) may be singular, whereas those formed using (6) are not.

 

In fluid flow problems, the coefficient matrix consists of a diagonal matrix with positive coefficients, typically representing the fluid compressibility within a calculation element, and a series of “dipoles” representing the transmissibility for fluid flow between elements. Each dipole comprises a positive coefficient added to the diagonal of A, and, to conserve material, an equal negative coefficient at an off-diagonal position in the same column of A. Usually, there are two dipoles for each pair of calculation elements, making a symmetric pattern of 4 added coefficients. However, the coefficients are not, in general, numerically symmetric, and it is this asymmetry that makes approximations (5) and (6) quite different in practice.

 

It is easy to see that A is diagonally dominant, and therefore non-singular, because the dipoles sum to zero in each column, leaving only the positive diagonal terms.

 

If the coefficients of D, L and U are d_ii, l_ij and u_ji for and , the colsum constraint is applied by calculating the modified diagonal coefficients, γ_i , using

........................................................................ (7)

where ............................................................................. (8)

 

For i=1, it is clear that γ₁ = d₁₁. For i=2, we need to evaluate η₁ - the sum of the coefficients below the diagonal in the first column divided by the first diagonal coefficient. Because of the way the matrix is constructed, we know that

.................................................................... (9)

 

so η₁ must take a value between 0 and -1. Then combining (7) and (9), we obtain

......................................................... (10)

 

so that γ₂ must satisfy

................................................................................. (11)

 

Next we see from the definition of η_k (8) that η₂ must also be between 0 and -1. The proof proceeds in this fashion from one column to the next, so that we obtain the general result

................................................................................. (12)

 

Thus γ is always positive, and always greater than the negative sum of the coefficients in the column below it. It follows that B is non-singular.

 

For the rowsum constraint, we use

......................................................................... (13)

where ............................................................................. (14)

 

and the arguments described above cannot by used, as there is no equivalent to inequality (9) for rows of coefficients.

A concrete example illustrates the point: given the matrix the colsum constraint produces whereas rowsum produces, and B is singular.