# On Incomplete Factorization Implicit Technique for 2D Elliptic FD Equations.

1 Introduction

The Incomplete Factorization iterative method (IF) for solving 2D and 3D elliptic finite-difference (FD) equations was suggested by Buleev [2]. In the method, schemes of factorization can be of Explicit (IFE) or Implicit (IFI) type, see [6]. The IFI scheme used here was proposed by Ginkin [5], and enhanced by Buleev [3].

Sabinin [8, 10] observed a similarity of iteration parameters for IFI and ones for Alternating Direction Implicit (ADI) method, and he proposed a cyclic set of parameters which permitted the effective solving boundary-value problems by IFI. He expanded this technique to 3D problems [9,11], and solved by it some mixed boundary-value problems of groundwater hydrodynamics with predominating of Neumann boundary conditions (e.g. [13]). Also, IFI was applied to a convective-dispersion equation [12], and to a problem of groundwater salt-heat transport [14]. Another modification of the iteration parameters was suggested in [7] and applied to some steady-state problems of underground water-oil flows.

Traditionally, IFE solvers are developed a lot (e.g. [1]) although IFI gives a better rate of convergence, as can be seen from the works [10], and [7].

In present paper, new variant of IFI technique is suggested which is different by applying the matrix algorithm of Thomas for enhancing stability. This is applied to as 5-point, as to high-order 9-point FD schemes of type "cross" for 2D elliptic problems. Also, it is suggested for IFI the cyclic set of iteration parameters, and the "hammer" sequence of parameters inside the cycle. Results of a numerical experiment are presented which demonstrate high effectiveness of IFI with such cyclic sets in boundary-value problems by Dirichlet, by Neumann, and in mixed ones for Poisson equation in rectangle, and in curvilinear (quasi-circular) areas. The results give a low difference between number of iterations in the Neumann problems, and in the Dirichlet problems, and a high effectiveness of IFI in solution of mixed boundary-values problems with predominating presence of Neumann boundary conditions.

The problem of obtaining the optimal set of iteration parameters for IFI is not yet analytically solved. Here, a numerical procedure for enhancing the ADI-type parameters is suggested.

The 9-point IFI solver is proposed to the high-order 9-point FD scheme of "cross" type. As this algorithm is complicated and bulky then it is suggested to use the 5-point IFI solver to solve the 9-point FD problems. The results of numerical experiment prove effectiveness of this idea.

2 Incomplete factorization implicit technique 2D

The 5-point FD scheme for 2D equations of type

div (K grad [phi]) + f = 0, (2.1)

where K = K(x, y, [phi]), and f = f(x, y, [phi]), gives for boundary-value problems a system of algebraic equations of the form A[phi] = f, where square matrix A is defined in grid nodes as

[mathematical expression not reproducible] (2.2)

where i = 0, ..., I; j = 0, ..., J; [mathematical expression not reproducible]; a, b, c, d [greater than or equal to] 0; [e.sub.ij] [greater than or equal to] [a.sub.ij] + [b.sub.ij] + [c.sub.ij] + [d.sub.ij].

Let the Dirichlet condition ([[phi].sub.ij] = const) be at least in one grid node with index i = [i.sub.0] (and j = [j.sub.0]), provided 0 [less than or equal to] [i.sub.0] [less than or equal to] I. The iterative process for solving (2.2) can be written as follows (s is the number of iteration):

[mathematical expression not reproducible] (2.3)

In general case 0 [less than or equal to] [i.sub.0] [less than or equal to] I, matrices L and U are defined as follows:

[mathematical expression not reproducible] (2.4)

where index i - 1 is for i < [i.sub.0], and index i + 1 is for i > [i.sub.0],

[Lv.sub.ij] [equivalent to] [v.sub.ij] for i = [i.sub.0]; (2.5)

[mathematical expression not reproducible] (2.6)

(index i + 1 is for i < [i.sub.0], and index i - 1 is for i > [i.sub.0]), and

[Uu.sub.ij] [equivalent to] ([A.sup.s] - [F.sup.s])[u.sub.ij] for i = [i.sub.0]. (2.7)

Equations for definition of unknown coefficients in matrices L and U are derived from the matrix equality LU = [A.sup.s] - [F.sup.s] + B, where matrix F is diagonal, with elements [F.sub.ij] = [([partial derivative]f/[partial derivative][phi]).sub.ij] [less than or equal to] 0, and matrix B is formed from two interpolation equations [mathematical expression not reproducible], where [omega] is the iteration parameter:

[mathematical expression not reproducible] (2.8)

where index i - 1 is for i < [i.sub.0], index i + 1 is for i > [i.sub.0].

As the result of these definitions, in the case of 0 < [i.sub.0] < I, one can obtain the recurrent formulas of opposite-currents algorithm for calculating the coefficients:

[mathematical expression not reproducible] ((2.9))

[mathematical expression not reproducible] ((2.10))

(index i - 1 is for i < [i.sub.0], and index i + 1 is for i > [i.sub.0]),

[mathematical expression not reproducible] (2.11)

In the special case [i.sub.0] = 0 (or [i.sub.0] = I), formulas for i < [i.sub.0] (or i > [i.sub.0]) are excluded from the algorithm.

The calculation of coefficients is made in following order. For all j = 0, ..., J, in the sequence of i = 0, 1, ..., [i.sub.0] - 1, I, I - 1, ..., [i.sub.0] + 1, the coefficients of L and U are calculated by (2.9)-(2.11), and the following equation is also solved (derived from (2.3)):

[Lv.sub.ij] = - [A.sup.s][[phi].sub.ij.sup.s] + [f.sub.ij.sup.s].

After that, one should solve equation [Uu.sub.ij] = [v.sub.ij], and calculate [[phi].sub.ij.sup.s+1] = [[phi].sub.ij.sup.s] + [u.sub.ij], as this follows.

In the case of 0 < [i.sub.0] < I, the system of three equations [Uu.sub.ij] = [v.sub.ij] for lines i = [i.sub.0] - 1, [i.sub.0], [i.sub.0] + 1 is solved by Thomas matrix algorithm, and after that, the remaining equations [Uu.sub.ij] = [v.sub.ij] are successively solved for values of i in the sequence of i = [i.sub.0] + 2, [i.sub.0] + 3, ..., I, [i.sub.0] - 2, [i.sub.0] - 3, ..., 0 by Thomas scalar algorithm.

For the case of [i.sub.0] = 0, the system of two equations [Uu.sub.ij] = [v.sub.ij] for i = [i.sub.0], [i.sub.0] + 1 is solved by the matrix algorithm of Thomas, and then the equations [Uu.sub.ij] = [v.sub.ij] are solved by the scalar Thomas algorithm for the rest values of i in the sequence of i = [i.sub.0] + 2, [i.sub.0] + 3, ..., I.

For the case of [i.sub.0] = I, the system of two equations [Uu.sub.ij] = [v.sub.ij] for i = [i.sub.0] - 1, [i.sub.0] is solved by the matrix algorithm of Thomas, and then the equations [Uu.sub.ij] = [v.sub.ij] are solved by the scalar algorithm of Thomas for the rest values of i in the sequence of i = [i.sub.0] - 2, [i.sub.0] - 3, ..., 0.

In the previous works (e.g. [7, 10]), I used another form of equation (2.7) which did not require the use of the matrix algorithm of Thomas, but was less effective.

It is obvious that one can apply the matrix algorithm of Thomas not only to three lines, as above, but to more lines, for example to five lines (i = [i.sub.0] - 2, [i.sub.0] - 1, [i.sub.0], [i.sub.0] + 1, [i.sub.0] + 2). Naturally, in this case, definitions (2.5) and (2.7) should be valid for lines i = [i.sub.0] - 1, [i.sub.0], [i.sub.0] + 1, and definitions (2.4) and (2.6) should be valid for lines i < [i.sub.0] - 1, and i > [i.sub.0] + 1. This 5-lines modification will be used here for solving the finite-difference problem in quasi-circle.

3 On the iteration parameter [omega]

The transition matrix T for transition to the next iteration ([[phi].sup.s+1] = [T.sup.s] [[phi].sup.s]) can be expressed as

T = [(A - F + B).sup.-1] B. (3.1)

Let us consider for simplicity a case of Poisson equation in a square-mesh FD grid, without dependence of function f on [phi], and provided [i.sub.0] = I for certainty. Then, one can define eigenvalues of T by substitution of exp (i[pi](m/I i + n/J j)) into (2.2) and (2.8), where i = [square root of -1], and estimate the absolute value [theta] of them. The result can be presented as

[theta] [less than or equal to] | [[lambda.sub.n]] - [OMEGA]/[[lambda.sub.n]] - [OMEGA] + R([mu] + [[lambda.sub.n]/[mu])]|.

Here, [mathematical expression not reproducible] is the absolute value of [mathematical expression not reproducible]. Using definitions (2.9) and (2.10), one can derive the continued fraction 2R + 1 = [rho] - 1/[rho]-..., where [rho] = 4[OMEGA] + 2, from which one can obtain the asymptotic value of R = [OMEGA] + [square root of [OMEGA] + [[OMEGA].sub.2]].

The eigenvalues [[lambda].sub.m] of second FD derivative respect to i are independent on the eigenvalues [[lambda].sub.n], therefore, let us evaluate them empirically by expression [mu] = (1 + [square root of 2])[OMEGA]/R, which gives 0 [less than or equal to] [mu] [less than or equal to] 1 if 0 [less than or equal to] [OMEGA] [less than or equal to] 1.

Then, as the result, one can obtain

([mathematical expression not reproducible] (3.2))

The right-hand side of (3.2) is the same which is used in the ADI method for obtaining optimal set of iteration parameters in the case of commuting operators [15]. Consequently, one can suppose that optimal parameters for IFI are linked with those for ADI, and use ADI-parameters for definition of [OMEGA] in IFI. As advantage, here in IFI is no requirement of commutation of operators, as it is in ADI. This lets to IFI be more effective in mixed boundary-value problems.

Following [15], the set of ADI parameters by Jordan [16] gives for IFI the cyclical set of increasing values [omega] (S is the period of cycle):

[[omega].sub.s] = 1 - 2[[OMEGA].sub.s], s = 0, ..., S - 1; [[omega].sub.s] = [[omega].sub.s]-S, s [greater than or equal to] S, (3.3)

where [mathematical expression not reproducible], where q = [[eta].sup.2] (1 + [[eta].sup.2]/2)/16, [sigma] = (2s + 1)/(2S), index c is the number of cycle.

The period of cycle can be estimated by empirical formula: S = [2ln(J)].

A solution of problem changes its spectrum after each applying the cycle. Here, this is taken into consideration by empirical parameter [b.sub.c]. For instance, in the experiment below, the set of values [b.sub.c]={1, 1/2, 2, 1/4, 4, 1/8, 8, ...} was used.

As [[alpha].sub.i-1] and [theta] are less when [omega] is less, then the sequence of values [[omega].sub.s] inside the cycle should be from minimum to maximum, in general. However, empirically, the best sequence proved be following the "hammer" principle: [[omega].sub.0], [[omega].sub.S-1], [[omega].sub.k], [[omega].sub.k+1], [[omega].sub.1], [[omega].sub.S-2], [[omega].sub.k-1], [[omega].sub.k+2], [[omega].sub.2], ..., where k is the integer part of S/2.

4 Iterative procedure for estimating optimal iteration parameters

One can estimate the optimal set of iteration parameters in following iterative procedure by using the least-squares method. Choosing a representative set of eigenvalues [[lambda].sub.n], and [[lambda].sub.m], one can write a functional of error from (3.1) as:

[mathematical expression not reproducible] (4.1)

For obtaining minimum of the functional, one should set the first derivatives from F respect to [[omega].sub.s] equal to zero. Resulting non-linear system of equations is transformed to the linear system for using in the iterative process (k = 0, 1, ... is the index of iteration):

([mathematical expression not reproducible] (4.2))

The iterative process (4.2) begins from initial values [[omega].sub.s.sup.0] = [[omega].sub.s] calculated by equations (3.3).

For example, values [omega] calculated by formulas (3.3), and by method (4.1)-(4.2) for S =10 are compared in Table 1. These values do not differ significantly.

This iterative procedure spends much time for obtaining the values of parameters and therefore is not recommended for practical using.

5 IFI technique for 9-point finite-difference scheme

The FD scheme (2.2) is constructed by applying a 3-point stencil in 1D to FD approximation of equation (2.1) respect to each coordinate. If applying the high-order 5-point stencil in 1D (see e.g. [4]) to this, one obtains the 9-point FD scheme "cross" which gives a system of algebraic equations of the form A[phi] = f for 2D elliptic boundary-value problems, where the square matrix A is defined in the grid nodes as

[mathematical expression not reproducible] (5.1)

where i = 0, ..., I; j = 0, ..., J; [mathematical expression not reproducible].

Let the Dirichlet condition ([[phi].sub.ij] = const) be at least in one grid node with index i = [i.sub.0], and 0 [less than or equal to] [i.sub.0] < I. The iterative process for solving (5.1) can be written as follows (s is the number of iteration):

LU([[phi].sub.ij.sup.s+1] - [[phi].sub.ij.sup.s]) = - [A.sup.s][[phi].sub.ij.sup.s] + [f.sub.ij.sup.s], s = 0, 1, ....

In the general case 0 [less than or equal to] [i.sub.0] [less than or equal to] I, matrices L and U can be defined by some several ways. Here, we present the most simple and symmetric variant:

[mathematical expression not reproducible]

where indexes i - 1 and i - 2 are for i < [i.sub.0], and indexes i + 1 and i + 2 are for

[mathematical expression not reproducible]

indexes i+1 and i+2 are for i < [i.sub.0], and indexes i - 1 and i - 2 are for i > [i.sub.0]; [Uu.sub.ij] [equivalent to] ([A.sup.s] - [F.sup.s])[u.sub.ij] for i = [i.sub.0].

Equations for definition of the unknown coefficients in L and U are derived, as above, from the matrix equality LU = [A.sup.s] - [F.sup.s] + B, where matrix F is diagonal, with elements [F.sub.ij] = [([partial derivative]f/[partial derivative][phi].sub.ij] [less than or equal to] 0, and matrix B is formed from eight interpolation equations [[phi].sub.i[+ or -]k,j[+ or -];l] + [[omega][phi].sub.ij] =[[phi].sub.i,j[+ or -]l] + [[omega][phi].sub.i[+ or -]k,j], where [omega] is the iteration parameter, and k, l = 1, 2.

As the result, one can obtain the recurrent formulas:

[mathematical expression not reproducible] (5.2)

indexes i - 1 and i - 2 are for i < [i.sub.0], and indexes i+1 and i + 2 are for i > [i.sub.0];

[mathematical expression not reproducible] (5.3)

The calculations by algorithm are made in the following order. For all j = 0, ..., J, in sequence of i = I, I - 1, ..., [i.sub.0] + 2, [i.sub.0] + 1, 0, 1, 2, ..., [i.sub.0] - 2, [i.sub.0] - 1, the coefficients of L and U are calculated by (5.2)-(5.3), and the following equation is solved, too:

[Lv.sub.ij] = -[A.sup.s][[phi].sub.ij.sup.s] + [f.sub.ij.sup,s]. (5.4)

Then, the system of five equations [Uu.sub.ij] = [v.sub.ij] for i = [i.sub.0] - 2, [i.sub.0] - 1, [i.sub.0], [i.sub.0] + 1, [i.sub.0] + 2 with index j = 0, ..., J is solved by the matrix five-diagonal Thomas algorithm along j to obtain [mathematical expression not reproducible], and [mathematical expression not reproducible].

Naturally, the system consists of five equations only if 1 < [i.sub.0] < I - 1. For [i.sub.0] = 1 or [i.sub.0] = I - 1, it becomes the system of four equations, and for [i.sub.0] = 0 or [i.sub.0] = I, it is of three equations.

Then, for all remaining lines i in the sequence of i = [i.sub.0] - 3, [i.sub.0] - 4, ..., 0, [i.sub.0] + 3, [i.sub.0] + 4, ..., I, the equation [Uu.sub.ij] = [v.sub.ij] is solved by the scalar five diagonal Thomas algorithm along j = 0, ..., J. Finally, [[phi].sub.ij.sup.s+1] = [[phi].sub.ij.sup.s] + [u.sub.ij] is calculated.

Iteration parameters should be corrected for the new scheme. This can be made by using in (3.3) the other definition [eta] = [sin.sup.2] ([pi]/[2b.sub.c]J)[(1 - 1/4 cos([pi]/2J))].sup.2].

6 Using 5-point IFI solver for 9-point FD scheme

Equation (5.4), and coupled to it [Uu.sub.ij] = [v.sub.ij] can be also solved with definitions (2.4)-(2.11) for L and U which use the FD scheme (2.2) instead of (5.1). It is possible if eigenvalues of operators B, and A which defined by (2.8), and (2.2) are close to those defined by (5.1). In such case, the transition operator (3.1) is distorted insignificant, and the using (3.3) for iteration parameters gives a convergence of iterative process. As shown in the experiment (see Section 9), this is so for Poisson equation.

7 Poisson equation in a rectangle area

In the numerical experiment below, it is made a one test of IFI techniques for the solution of Neumann-type boundary value problem for Poisson equation in a square.

Let [phi] = [x.sup.2][z.sup.2], then as the result, normal derivatives at the boundaries are [partial derivative][phi]/[partial derivative]/x = 2x[z.sup.2] and = 2[x.sup.2]z, and the Poisson equation is

[[partial derivative].sup.2][phi]/[partial derivative][x.sup.2] + [[partial derivative].sup.2][phi]/[partial derivative][z.sup.2]=2([x.sup.2] + [z.sup.2]). (7.1)

The 5-point FD scheme (2.2) and the 9-point FD scheme (5.1) are precise schemes for (7.1), that is [[phi].sub.ij] [equivalent to] [phi], or the FD solution is equal to the preset solution in all points of FD area including the boundary.

Let the grid size I = J. At one boundary node, the FD value is fixed: [[phi]i.sub.0],[j.sub.0] = [phi]. As an initial distribution (s = 0) in all other nodes, a step-wise function is used as the mostly inconvenient one: [[phi].sub.ij.sup.0] = [phi] + 1 for i + j < I, and [[phi].sub.ij.sup.0] = [phi] - 1 for i + j [greater than or equal to] I. The fastest convergence is observed when [i.sub.0] = [I/2].

Another test is made for the solution of Dirichlet boundary value problem for Poisson equation in the square. Here, the FD values are fixed at all boundary nodes. In this case, the 5-point and 9-point FD schemes are precise also for a cubic solution: [phi] = [x.sup.3][z.sup.3], provided the Poisson equation is [[partial derivative].sup.2][phi]/[partial derivative][x.sup.2] + [[partial derivative].sup.2][phi]/[partial derivative][z.sup.2] = 6xz([x.sup.2] + [z.sup.2]).

8 Poisson equation in a quasi-circular area

As an example of application of IFI iterative technique to non-rectangular areas, the Neumann-type boundary value problem for Poisson equation (7.1) is considered for an area which is a step-wise approximation to a circle, with [partial derivative][phi]/[partial derivative]x = 2x[z.sup.2] at every vertical segment of the boundary, and [partial derivative][phi]/[partial derivative]z = 2[x.sup.2]z at every horizontal segment of the boundary. The segment of step-wise boundary is the grid piece which is nearest to the circle in both i and j directions. At the one boundary node shared with the circumscribed square, the boundary value [[phi].sub.ij] is fixed.

9 Experiment. Comparison of solvers

The quality of iteration process is observed here on two values: the discrepancy of function [mathematical expression not reproducible], and the specific residual of equation [mathematical expression not reproducible].

In Figures 1- 3, curves of quality, that are graphics of the residual function - ln(r) in dependence on s/ln(J), are presented for the Dirichlet boundary-value problem for Poisson equation in square. An ideal graphic seems to be a straight line. Figure 1 corresponds to the variant of the 5-point solver and the 5-point FD scheme (2.2). Figure 2 corresponds to the variant of the 9-point solver and the 9-point FD scheme (5.1).

In Table 2, the number of IFI iterations is presented for variant of Figure 1, and different error levels [epsilon] (i.e. until it becomes r [less than or equal to] [epsilon]), in comparison with the theoretical number of iterations taken from [15] for ADI and A-T (Alternate-Triangular) methods. In Table 3, the number of IFI iterations is shown for two error levels, and the 9-point solver for Dirichlet problem, variant of Figure 2.

Figure 3 corresponds to the variant of the 5-point solver and the 9-point FD scheme. In Figure 4, graphics of the discrepancy function - ln(d) in dependence on s/ln(J) are shown for the case of Figure 3.

The quality of solvers in Neumann boundary-value problems for Poisson equation in square generally depends on position of [i.sub.0]. Results for the best value [i.sub.0] = I/2 are presented in Figures 5-8. Dependence on [i.sub.0] is considered in Section 11. Figures 5-7 show graphics of function - ln(r) in dependence on s/ln(J), and Figure 8 shows graphics of function - ln(d) in dependence on s/ln(J).

Figure 5 corresponds to the variant of the 5-point solver and the 5-point FD scheme (2.2). Figure 6 corresponds to the variant of the 9-point solver and the 9-point FD scheme (5.1).

In Table 4, the number of iterations is shown for different error levels, variant of Figure 5. In Table 5, the number of iterations is shown for different error levels, variant of Figure 6.

Figure 7 corresponds to the variant of the 5-point solver and the 9-point FD scheme. In Figure 8, graphics of function - ln(d) in dependence on s/ln(J) are shown for the case of Figure 7.

In Table 6, the number of iterations is shown for different error levels, variant of Figure 7. The comparison of Tables 5 and 6 shows that the 5-point solver behaves better than the 9-point one. The comparison of Tables 2 and 4 shows that the 5-point solver practically doesn't depend on the type of boundary-value problem.

From comparison of Figures 1-3, 5-7 and Tables 2-6, one can conclude that the quality of iteration processes, and number of iterations do not significantly differ on Dirichlet and Neumann boundary-value problems, and on the 5-point and 9-point IFI solvers in the case of rectangle area. It should be also noted that convergence by the discrepancy function is slower than by the residual one, especially for large J.

10 Experiment. Quasi-circle area

In Figure 9, the quasi-circle area for J = 500, and equidistant isolines of calculated solution [[phi].sub.ij] are shown.

In Figures 10-11, graphics of quality functions - ln(r), and - ln(d) on s/ln(J) are shown for the 5-point solver and the 5-point FD scheme in the Neumann boundary value problem for Poisson equation in quasi-circle, [i.sub.0] = I/2, [j.sub.0] = 0. It was used the 5-lines modification of 5-point solver because the more simple 3-lines algorithm was considerably slower. The parameter [b.sub.c] was set in [pi]/2 times more than [b.sub.c] for the square.

In Table 7, the number of iterations is shown for different error levels, variant of Figure 10. Data of Table 7 and Figure 10 show small slowing the iterative process in comparing with Table 4 and Figure 5 for the case of rectangle area.

11 Experiment. Mixed boundary-value problem

Close results of the Dirichlet and Neumann problems let us suppose that solution of the mixed boundary-value problems (especially which are close to the Dirichlet or Neumann problems) by IFI must be also economical. However, the high efficiency is obtained only at high predominance of Dirichlet or Neumann areas in the all boundary. This can be explained by that the spectra of FD mixed problems have not the homogeneity of the pure Dirichlet or Neumann problems.

The similar spectral situation is seemed to be in Neumann problems when the parameter [i.sub.0] takes a non-symmetrical position: [i.sub.0] [not equal to] I/2. It is really so, the number of iterations depends on [i.sub.0]. I tried to decrease the dependence by a choice of sequence for [b.sub.c] in formula (3.3). In Figure 12, the number of iterations is presented for reaching r < [10.sup.-10], and r < [10.sup.-6] as the function of position ([i.sub.0]/I) in percents. It was set [b.sub.c] = {4, 8, 1/4, 1/2, 1, 1, ...}. The solution was obtained by the 5-point solver at the 5-point FD scheme.

One can notice from comparison with Section 9 that achieving a relative constancy of results on [i.sub.0] entails increasing number of iterations at small values of J.

Let us consider now the mixed boundary-value problem for Poisson equation in the square, as above, at three sides of which have been set Neumann conditions, and at the fourth one (j = J) has been set the Dirichlet condition at [i.sub.0] < i < [i.sub.00], and the Neumann condition at the rest. Let's [i.sub.0] = I/2, and the parameter [i.sub.00] changes its value from [i.sub.0] till I. We will solve the 5-point FD scheme by the 5-point solver. It is set [b.sub.c] = {1, 3/4, 1, 0.3, 1/2, 1/4, 1/4, ...}.

In Figure 13, one can see the number of iterations as the function of the Dirichlet boundary ([i.sub.00]/I) in percents. The curves are far from constants, especially for large values of J. However, the increasing number of iterations in area of the nearest to [i.sub.0] values of [i.sub.00] (up to 10%) is not large. The left side of graphics corresponds to the case of Neumann boundary-value problem ([i.sub.00] = [i.sub.0]).

12 Discussion and conclusions

The suggested IFI technique is powerful for reducing computational expense due to small number of iterations, and weak dependence on type of boundary conditions. It is clearly seen from Figures 1-8, and Tables 2-3.

The 5-point IFI solver effectively solves also 9-points FD equations of type (5.1).

The IFI technique is suitable also for non-rectangle areas. However, it can be slower than in the rectangular areas, what can require increasing number of lines used in the matrix algorithm of Thomas.

IFI techniques are at least as fast as the best-known iterative methods (ADI, A-T), and can be useful for solution of different mixed boundary-value problems. The mixed problems with the predominant presence of Neumann boundary conditions are characterized by slow convergence of many iterative methods of solution. IFI slows down its convergence non-significantly.

Unfortunately, the proposed method for calculating the iteration parameter [omega] is not universal, and requires a presence of empirical parameters [b.sub.c]. This limits the effectiveness of applying IFI to complex problems, and requires further work to improve the IFI technique.

The interesting problem is how many lines for matrix Thomas algorithm is optimal: 3, 5, 7 or more. The increase in number of lines improves the convergence but increases computational time. This could be a topic for further IFI research.

The formulas of IFI algorithm in Section 2, and 5 are written for quasi-linear FD equations. Linearization of matrix A is done by Picard's method, and of the right part f - by Newton's method. It is the easiest, stable, and most natural way. It was successfully applied to the groundwater problems [8, 9, 11]. Application of Newton's method to matrix A is unstable.

As in quasi-linear equations, IFI is successfully applied to problems with significantly varying coefficients, and to complicated ones. One can see examples in [7, 8, 9, 11, 12, 13, 14]. In such problems, behavior of IFI depends on how lucky setting empirical parameters [b.sub.c].

The suggested IFI technique can be expanded to 3D problems by a way of works [9] and [11]. As the factorization here is made by three coordinates then there are two ways to introduce the lines for matrix Thomas algorithm: 1) in the third stage of factorization, or 2) in the second and third stages. In the first case, an area of the lines for matrix Thomas algorithm is a plane; in the second case, it is a volume. The 3D algorithms could be a topic for further IFI research, too.

Comment to references

Articles [2, 3, 5, 8, 9, 10, 11, 12, 14] are published in Russian. Scarce articles of the author [8, 9, 10, 11, 12, 13] can be found at the web-address https://yadi.sk/d/plXWGEUGJRIkWQ.

References

[1] H. Anzt, E. Chow, J. Saak and J. Dongarra. Updating incomplete factorization preconditioners for model order reduction. Numerical Algorithms, 73, 02 2016. https://doi.org/10.1007/s11075-016-0110-2.

[2] N.I. Buleev. The incomplete factorization method for solution of 2D and 3D finite-difference equations of diffusion type. JVM and MF, 10(4):1042-1044, 1970 (in Russian). https://doi.org/10.1016/0041-5553(70)90027-3.

[3] N.I. Buleev. The new variant of incomplete factorization method for solution of 2D finite-difference equations of diffusion. Novosibirsk, Chislennye metody mechaniki sploshnoi sredy, 9(1):5-19, 1978 (in Russian).

[4] F. Gibou and R. Fedkiw. A fourth order accurate discretization for the Laplace and heat equations on arbitrary domains, with applications to the Stefan problem. Journal of Computational Physics, 202.2:577-601, 2005. https://doi.org/10.1016/j.jcp.2004.07.018.

[5] V.P. Ginkin. The h-factorization method for solution of 2D finite-difference equations of diffusion. Novosibirsk, Vychislitelnye metody lineynoi algebry, pp. 123-132, 1977 (in Russian).

[6] V.P. Il'in. Iterative Incomplete Factorization Methods. World Scientific, 1992. https://doi.org/10.1142/1677.

[7] V. Sabinin. An incomplete factorization technique for fast numerical solution of steady-state ground-water flow problems. Geofisica Internacional, 44(3):275-282, 2005.

[8] V.I. Sabinin. The numerical solution of the horizontal systematic drainage problem with zone of incomplete saturation. Novosibirsk, Dinamika sploshnoi sredy, 46:122-136, 1980 (in Russian).

[9] V.I. Sabinin. The numerical solution of the 3D filtration problem with incomplete saturation. Novosibirsk, Dinamika sploshnoi sredy, 51:129-141, 1981 (in Russian).

[10] V.I. Sabinin. On one algorithm of the method of incomplete factorization. Novosibirsk, Chislennye metody mechaniki sploshnoi sredy, 16(2):103-117, 1985 (in Russian).

[11] V.I. Sabinin. The problem of computation of drainage of rice fields. Novosibirsk, Dinamika sploshnoi sredy, 81:103-116, 1987 (in Russian).

[12] V.I. Sabinin. Computer prognosis of ground-water contaminant transport. Novosibirsk, Dinamika sploshnoi sredy, 108:51-62, 1994 (in Russian).

[13] V.I. Sabinin. Ground water and slope surface flows numerical modelling. Modern Approaches to Flows in Porous Media, Int. Conf. dedicated to P. Ya. Polubarinova-Kochina (1899-1999). Abstracts. Moscow, pp. II-36-II-38, 1999.

[14] V.I. Sabinin. The solution to a problem of ground-water salt-heat transport by an incomplete factorization technique. Siberian J. of Numer. Mathem. / Sib. Branch of Rus. Acad. of Sci., Novosibirsk, Russia, 2(1):69-80, 1999 (in Russian).

[15] A.A. Samarsky and E.S. Nikolaev. Numerical Methods for Grid Equations. Volume II. Iterative Methods. Birkhauser Verlag, 1989. https://doi.org/10.1007/978-3-0348-9142-4.

[16] E.L. Wachspress. Optimum alternating-direction-implicit iteration parameters for a model problem. J. Soc. Indust. Appl. Math., 10:339-350, 1962. https://doi.org/10.1137/0110025.

Instituto Mexicano del Petroleo Eje Central Lazaro Cardenas 152, Col. San Bartolo Atepehuacan,

Gustavo A.Madero. C.P. 07730, Mexico D.F., Mexico.

E-mail: vsabinin@yahoo.com

Received February 18, 2019; revised October 20, 2019; accepted November 1, 2019

https://doi.org/10.3846/mma.2020.8485

Caption: Figure 1. 5-solver, 5-scheme, Dirichlet.

Caption: Figure 2. 9-solver, 9-scheme, Dirichlet.

Caption: Figure 3. 5-solver, 9-scheme, Dirichlet.

Caption: Figure 4. 5-solver, 9-scheme, Dirichlet, d-function.

Caption: Figure 5. 5-solver, 5-scheme, Neumann.

Caption: Figure 6. 9-solver, 9-scheme, Neumann.

Caption: Figure 7. 5-solver, 9-scheme, Neumann.

Caption: Figure 8. 5-solver, 9-scheme, Neumann, d-function.

Caption: Figure 9. Equidistant lines of equal head for the quasi-circular area.

Caption: Figure 10. 5-solver, 5-scheme, Neumann for circle.

Caption: Figure 11. 5-solver, 5-scheme, Neumann, circle, d-function.

Caption: Figure 12. Number of iterations on position of [i.sub.0] in the Neumann boundary-value problem.

Caption: Figure 13. Number of iterations on the size of Dirichlet area in the mixed boundary-value problem.
```Table 1. Optimal sets of iteration parameters for J=200, and S=10.

Eq.(3.3)  -0.7280  0.99097  0.2673  0.9970  0.7503  0.9990  0.9173
Eq.(4.2)  -0.7278  0.99098  0.2679  0.9970  0.7445  0.9990  0.9173

Eq.(3.3)  0.9997  0.9727  0.99986
Eq.(4.2)  0.9997  0.9726  0.99986

Table 2. Number of iterations on J. Dirichlet problem. In IFI, the
5-point solver, 5-point scheme is applied.

J                     50  200  500  2000

IFI, [epsilon] = [10.sup.-10]  25   37   44    54
IFI, [epsilon] = [10.sup.-6]   15   21   23    29
ADI, [epsilon] = [10.sup.-6]   13   17   20    24
A - T,[epsilon] = [10.sup.-6]  29   58   91   182

Table 3. Number of iterations on J. Dirichlet problem, the 9-point
solver, 9-point scheme.

J  50  200  500  2000

[epsilon] = [10.sup.-10]  27   39   44    56
[epsilon] = [10.sup.-6]   15   21   24    31

Table 4. Number of iterations on J. Neumann problem, the 5-point
solver, 5-point scheme.

J  50  200  500  2000

[epsilon] = [10.sup.-10]  27   39   44    55
[epsilon] = [10.sup.-6]   15   21   23    29

Table 5. Number of iterations on J. Neumann problem, the 9-point
solver, 9-point scheme.

J  50  200  500  2000

[epsilon] = [10.sup.-10]  28   40   46    57
[epsilon] = [10.sup.-6]   16   21   25    31

Table 6. Number of iterations on J. Neumann problem, the 5-point
solver, 9-point scheme.

J  50  200  500  2000

[epsilon] = [10.sup.-10]  27   39   44    55
[epsilon] = [10.sup.-6]   15   21   24    29

Table 7. Number of iterations for Neumann problem in quasi-circle.

J  50  200  500  2000

[epsilon] = [10.sup.-10]  27   40   47    57
[epsilon] = [10.sup.-6]  16   22   25    31
```