CC BY
9
DOI: 10.14529/cmsel90201
THE USE OF THE LINE-BY-LINE RECURRENT METHOD FOR SOLVING SYSTEMS OF DIFFERENCE ELLIPTIC EQUATIONS WITH NINE-DIAGONAL MATRICES
© 2019 A.A. Fomin1, L.N. Fomina2
1T.F. Gorbachev Kuzbass State Technical University (Vesennaya 28, Kemerovo, 650000 Russia),
2Kemerovo State University (Krasnaya 6, Kemerovo, 650043 Russia)
E-mail: fomin_aa@mail.ru, lubafomina@mail.ru Received: 21.11.2018
The applying of the line-by-line recurrent method for solving systems of difference elliptic equations with nine-diagonal matrices is the subject of the article. Such matrices take place in the case of difference approximation of 2D differential problems of a higher order of accuracy on a regular grid covering the area under consideration. The technology of the so-called compensatory transform which allows replacing the initial nine-diagonal matrix of the system with the five-diagonal one is offered in the article, due to the fact that originally the line-by-line recurrent method was designed for solving systems of difference equations with a five-diagonal matrix. The efficiency of this technology is analyzed by comparing the solutions of the test boundary value problem in a unit square. The solutions are found both with the help of different implementations of the compensatory transform technology and by other modern highly efficient iterative methods for solving the systems of difference equations. The problem is solved on the sequence of grids from coarse (501 x 501) to fine (4001 x 4001) nodes. The accuracy of the solution convergence is determined by the relative norm of the residual, which is equal to 10“12 in the present work. It is shown that the line-by-line recurrent method retains its high efficiency over the entire range of the grids under consideration despite the use of the intermediate technology of the compensatory transform.
Keywords: grid method, system of difference elliptic equations, iterative method, convergence of solution.
FOR CITATION
Fomin A.A., Fomina L.N. The Use of Line-by-Line Recurrent Method for Solving Systems of Difference Elliptic Equations with Nine-Diagonal Matrices. Bulletin of the South Ural State University. Series: Computational Mathematics and Software Engineering. 2019. vol. 8, no. 2. pp. 5-21. (in Russian) DOI: 10.14529/cmsel90201.
Introduction
As is well known, many “standard” finite-difference approximation technologies for two-dimensional differential equations of problems of fluid dynamics and heat transfer on regular grids are based on five-point stencil. The Patankar scheme [1] can be mentioned here as an appropriate example of this kind of technologies. This technology yields an algorithm which strictly provides monotonicity and central-point dominance of the finite-difference five-point scheme. And as a result, a system of linear algebraic equations (SLAE) with a five-diagonal matrix of positive type arises [2]. The Patankar scheme has shown itself well in numerous studies. However, the need for higher-order accuracy schemes arises in a number of cases. For example, the five-point central-difference scheme provides a smaller error than the Patankar scheme, despite the fact that these schemes are formally of the same order of approximation. Unfortunately, the SLAE matrix obtained on the basis of the central-difference scheme can lose its positive type in certain conditions.
One can use the nine-point stencil (five points along each direction) to avoid this difficulty and increase the approximation order at the same time. In this connection a variety of higher-order
numerically stable schemes was developed using the nine-point stencil [3-8]. As a rule, the convective flux discretization is the focus of a higher-order difference approximation, as it was done in the previously cited works. While, the diffusion part of the transport equations is usually approximated by the standard five-point scheme of the second order accuracy. The works that use a higher-order difference approximation of diffusion flux are found much less often [9]. Difference approximation of differential equations with derivatives of higher order (up to the fourth one) is another origin of difference linear equations on five points along each independent coordinate [10]. In this case, it is evident that for 2D problems the difference stencil will consist of nine nodes.
It is remarkable that the absence or cursory mention of an iterative method for solving arising nine-diagonal SLAE is a common feature of the above works. Probably the researchers believe that the expansion of existing methods for solving systems with five-diagonal matrices to methods for solving systems with nine-diagonal ones is not an insuperable barrier. Indeed, the transition of such well-known methods as block successive over-relaxation (BSOR) or line-by-line method [1] from five-point to nine-point versions do not cause much difficulty. But on the other side, there are effective computational technologies for the SLAE solution, which don’t have an easy expansion up to nine-point stencils (see, for example, [11, 12]).
A deferred correction procedure was developed to overcome this problem. Its idea is that a nine-diagonal matrix of the system of equations is replaced by a five-diagonal one which is obtained with the use of a difference approximation of lower order in a special way. After that one can use any available method for solution of SLAE with a five-diagonal matrix. Interestingly, the procedure occurs in two variants in the literature. The first variant is used when there is an independent difference approximation of the convective and diffusive terms of the transport equation [6]. So, it operates with convective terms of the equation on the stage of their difference approximation. The second variant of the procedure is applied to an already formed algebraic difference equation, regardless of how it was obtained and, therefore, it is more universal [13]. The article will consider the second version of the deferred correction procedure. It should be noted that the procedure has a weak point, namely: if there is no difference approximation of some terms of the equation of lower order on a reduced number of nodes, then it is not applicable. For example, such situation takes place when there are fourth-order derivatives in a differential equation that cannot be approximated using less than five nodes. Therefore, in this case one have to design other procedure to replace a nine-diagonal matrix of SLAE with the five-diagonal one, which does not have this shortcoming.
Recently, a highly efficient line-by-line recurrent method was developed and successfully used for solving systems of difference elliptic equations in some problems of computational fluid dynamics and heat transfer [14]. This method is applicable to systems of equations with five-diagonal matrices due to its design features in the case of two-dimensional problems [15]. However, attempts to modify this method for the case of systems of equations with nine-diagonal matrices faced great difficulties. Therefore, the way out is to use the line-by-line recurrent method in cooperation with an intermediate technology of matrices transformation like the deferred correction procedure, for example. In the light of the foregoing, the objective of the work is to develop a universal technology of matrix transformation and to investigate the effectiveness of the line-by-line recurrent method in solving systems of difference elliptic equations with nine-diagonal matrices.
The paper is organized as follows. The mathematical statement of a problem, namely: definition area, differential transport equation, boundary conditions, closing formulas — are described in Section 1. In the following Section 2, the details of the numerical technique including high order numerical discretization, deferred correction procedure, compensatory transform technology are given. The comparisons of solutions of the test problem obtained by different methods are presented in Section 3. The conclusions are drawn in the final Section.
1. Statement of problem
Let О = {(ж, у) : 0 < ж < 1, 0 < у < 1} be a unit square in Cartesian coordinates (Fig. 1) as the definition area of unknown Ф(ж, у) which is governed by a differential transport equation.
Fig. 1. Scheme of the problem area
In this case the formulation of a test 2D boundary value problem in О can be used as an origin for obtaining a system of difference elliptic equations with the sparse matrix with the help of difference approximation of the initial differential problem. The generalized steady-state convection-diffusion transport equation written for Ф(x,y) can be stated as [1]
дФ дФ д / <9Ф\ д / <9Ф
дх ду дх у дх) ду \ ду
(1)
where U(x,y), V(x,y) — flow velocity components, Г(ж,у) — transfer coefficient, S(x,y) — source. Dirichlet conditions take place on the area boundaries. Let velocity components and transfer coefficient be as follows:
у3
U{x,y) = -3y2arctanx, V{x,y) = ——5-; Г(ж, y) = exp(-l2), l2 = x2 + y2.
1 + жz
It should be noted that the velocity field is solenoidal one. Lastly, let the solution of the test problem be the function u(x,y) = exp(—Wl2) cos(87tI2), then the substitution of u, U, V, and Г in equation (1) makes it possible to define the expression for the source S(x, y). It is not difficult to see that Dirichlet conditions at the area boundaries in the case of the u(x, y) are written as:
0 < у < 1 : Ф(0,у) = exp(-10y2) cos(87ry2), Ф(1,у) = exp(-10(l + y2)) cos(87t(1 + y2));
0 < x < 1 : Ф(х, 0) = exp(—10ж2) cos(87ra:2), Ф(ж, 1) = exp(—10(1 + ж2)) cos(87t(1 + ж2)).
So, the test 2D boundary value problem is defined and one can begin to solve it numerically.
2. Numerical technique
2.1. High-resolution numerical discretization
The problem area is covered by a uniform orthogonal mesh which nodes one can separate into three groups: internal, near-boundary, and boundary nodes. The nine-point stencil with an internal central node is presented in Fig. 2a. The so-called SMART scheme [7] is used for higher-order approximation of (1) in internal mesh nodes. Briefly, the technology of the scheme
obtaining is as follows [16]. Integrating the equation (1) over the control volume (Fig. 2a) and using the divergence theorem for a Cartesian coordinate system allows getting the following discrete equation:
Je Jyo + Jn Js — Qi (2)
where Je, Jw, Jn, Js represent the total fluxes of unknown Ф across faces e, w, n, s of the control volume, and Q is the volume integral of the source term S. Each of the surface fluxes J contains convective and diffusive contributions. It is expressed, for example, for the face e, as follows:
where
= ил, - Г,
Ф p — Фр
6x,
А У,
(3)
% + (Фе-%) /(Фр), Ue> О,
Фее + (Фр — Фее) /(Фе), Ue < 0;
0 < Ф < 1/6,
1/6 < Ф < 5/6,
5/6 < Ф < 1,
Ф elsewhere.
The tilde above unknown Ф denotes a so-called normalized variable which is defined as Ф = (Ф — Фр)/(Фр — Фр). Here, subscripts U, D mean upflow and downflow nodes relative the central point, correspondingly. For example, for the face e the central node will be point P in the case Ue > 0 while conversely - point E if Ue < 0. So, Ф = Фр = (Фр — Фгу)/(Фр — Фw) for Ue > 0, and Ф = Фр = (Фр — Фрр)/(Фгу — Фее) for Ue < 0.
The fluxes through the w, n and s faces can be found in a similar manner.
Ue±Up r 2JjTp_ ф
2 ’ Гр + Гр’
/(Ф) =
ЗФ,
3/8 +3/4 Ф,
1,
Ф,
^5xWv Me Ik
nn\ 4
J ^ 7 1 у Nf J
J ^ 7 5 уп j ~7 5 ys 7 + 1 > V X t ww гЦ J n
У- m i % m e< У
s S i
7-UO A J
Ax
п- -3 n .s+Y -2 n—l 7
л:
b) near-boundary mesh node
Fig. 2. Finite-difference stencil for higher-order discrete approximation
The eight-point stencil with a near-boundary central node is presented in Fig. 2b. For example, the right boundary of the problem area H is chosen. For this case the general approximation scheme described above can be applied if Ue > 0. In contrary, if Ue < 0 (as in the figure) it is necessary to follow special practices. In this context, it should be kept in mind
that a lower order approximation takes place as a result of the “windward rule” in any case. In literature the first order upwind scheme for near-boundary node is used as a rule [6]. But in the present work the more complex technology of Patankar scheme is adapted with a view to obtain a second order upwind difference scheme. For this case it is not difficult to get the approximation formula for Фе based on this methodology applying the profile of the fifth degree for unknown Ф, namely:
Фе = ФP + (Фе - Фр) v(Pe),
where, in the context of the grid uniformity
<P(Pe) =
Ф(Ре
2Ф(Ре/2)
Pe =
(Гр + TE)(Up + UE)
4ГРГр
Sxe]
ЩРе) =
Pet Pe < -10,
(l + 0,lPe)5-Pe, -10<Pe<0,
(1 - 0, LPe)5, 0 < Pe < 10,
0, Pe > 10.
Indeed, as is well known, the solution of the problem “convection and diffusion”
(4)
(5)
(6)
d_
dx
(рИФ)
±(т^.
dx V dx
for a domain 0 < x < L with boundary conditions: Ф = Фд at x = 0, and Ф
Фр at x = L is
Ф = Ф0 + (Фр
Фо)
exp(Px/L) — 1 exp(P) — 1
(7)
on the assumption with Г and pU are constants [1]. Here P is a Peclet number defined by P = pUL/T.
The value of Фе in (4) is calculated according to the solution profile (7), i. e. it is assumed that Фо = Фр, Фр = Фр, L = dxe, x = 8xe/2, P = Pe, and Ф = Фе in the formula (7). So,
Фе
Фр + (Фр — Фр)
ехр(.Ре/2) - 1
exp(Ре) ~ 1
Фр + (Фр — Фр)
1 Ре ехр(Ре/2) - 1
2 ехр(Р>е) — 1 Ре/2
or
Фе
Фр + (Фр — Фр)
Ф(Ре)
2Ф(Ре/2)’
(8)
where Ф(г) = z/(exp(z) — 1). Because an exponential function is very expensive to compute, Ф(г) is approximated by Patankar’s power-law scheme (see formulas (5.27) and (5.33) in [1]) which is represented by the complex formula (6) in the present work. In other words, Ф ~ Ф. As a result, it is easy to see that in this case the formulas (8) and (4) are almost identical taking into account the formula (5). What was required to show.
Finally, the trivial “approximation” takes place for the third group of the mesh (i. e. for the boundary lines of Q) because of Dirichlet boundary conditions in the problem.
2.2. Deferred correction procedure
The deferred correction (DC) method is a simple and proven procedure that enables the use of high order approximation schemes in codes initially written for low order schemes. Let, in
general case, there be a difference scheme as a result of approximation of the original differential equation (1) on the nine-point stencil (Fig. 2a) of the following kind
аРФр = чеФе + aw&w + aN^N + as^s + чееФее + aww^ww + aNN^NN + assess + b. (9)
In turn, let the difference scheme of lower order approximation for the same equation (1) be as follows
арФр = ОеФе + о,цгФ\у + адтФдг + Q's'hs + b. (19)
It is easy to see, adding to both sides of equation (9) the combination of арФр — ^)а^ьФпь,
nb
composed of the terms of equation (10), the DC procedure results in a five-point equation
арФр+1 = ^ + ’S^ianb ~ апъ)ФпЬ + ^ аппъФппЪ + (aP ~ «р)Фр + b, (11)
nb nb nnb
where к is number of iteration, nb = {E, W, N, S}, nnb = {EE, WW, NN, SS}. It is clear that the solution of equation (11) tends to the solution of equation (9) with the convergence of the
iterative process (i. e., with Ф^+1 -у Ф^). At the same time, one can use any previously
k^oo
created methods for solving SLAE with five-diagonal matrices to solve modified system on the base of equation (11).
Further in the article the DC procedure will be denoted as DSPt5, since the Patankar scheme with the profile of unknown Ф of the fifth degree is used to the lower order approximation.
2.3. Compensatory transform technology
As was mentioned above, the DC method is usable when a lower order approximation on a truncated five-point stencil takes place. Otherwise, one must apply other more general technology to transform a nine-diagonal matrix of SLAE to a five-diagonal one. Precisely that kind of a procedure, the so-called compensatory transform technology, is offered in the work. The major idea of the compensatory transform technology is to express the iterative increment of the sought-for solution in the “extreme” nodes of the stencil (in the Fig. 3 they are marked in cyan) through the increment in the “internal” nodes (white and black nodes of the stencil). So, the “extreme” nodes of the stencil are excluded and the matrix of the system of equations is transformed from nine-diagonal to five-diagonal one.
j+2
7+1
7-1 7-2
i—2 i— 1 i г+l i+2
/
V, к /: J к
V. W С к
ч J S
Fig. 3. The scheme of the compensatory transform of the nine-point stencil into five-point one
The transformation formula has the first or second order of accuracy depending on the number of the “internal” nodes of the difference stencil used in the expression. For example, in the case of uniform grid the formula of the first order accuracy for node EE is as follows
Дф|У = в (2ДФ^+1 - Дф*+1) ,
and the formula of the second order of accuracy for the same node is as follows
3 (ДФ|+1 - ДФр+1) - ДФ^1 .
Here ДФ^+1 = ф/г+1 — Ф^ - is increment of the sought-for solution, 9 is a parameter of compensation, which should be in the range 0 < 9 < 1 [2]. It is easy to verify that the application of the above formulas will lead out to the following expressions for the transformed coefficient in the nearby point E
aE = dE + 2 9ciee,
dE = cie + 0 (3йЕЕ + dww)
for the first and second order of accuracy respectively. Transformed coefficients for other nearby points W,N,S are written in a like manner. As a result, the transformed five-point difference equation is arrived as follows
Дф|^ = в
аРФкр+1 = Y^dnb<S>k+1+Ъ, (12)
nb
where for the first order of accuracy
dp = ap + 9 (dEE + dwW + dNN + ass) 5
b = b + aEE [ФкЕЕ ~ 9 (2Ф1 - фр)] + o,ww [ф^т " 9 (2фт “ фр)] +
+ aNN [Ф*^ - 9 (2Ф% - Фкр)] + aSS [Ф|<? " 9 (2Ф! " фр)] ,
and for the second order of accuracy, respectively
ap = ap + 39 (apE + dww + dNN + ass),
b = b + aEE {ФкЕЕ - 9 [3 (ф| - фр) + фт] } + aww {$kww ~ 9 [3 (фт - фр) + ф1] } +
+ aNN {Ф^я " 9 [3 (фя - фр) + ф|] } + ass {Ф|я " 9 [3 (ф! - фр) + фя] } •
In the further, the compensatory transform technique of the first order of accuracy will be denoted as Cl and of the second order - as C2.
3. Computed results and discussion
3.1. Nomenclature of methods and the research strategy
In general, eight different methods are used to solve the problem formulated in the first section of the article. Nomenclature of methods (abbreviations and their expansions) is presented in Tab. 1. The solution of the problem is calculated with five uniform grids of different resolution: 501 x 501, 1001 x 1001, 2001 x 2001, 3001 x 3001, 4001 x 4001. Thus, the number of unknowns in generating SLAE varies from about 25 x 104 (coarse mesh) to 16 x 106 (fine mesh).
The research strategy is a comparative analysis of the characteristics of the convergence of the methods for solving SLAE (see Tab. 1) for each of the mesh partitions of the problem domain H.
Table 1 Iterative methods for SLAE solutions construct
No. Type of transform Method Abbreviation expansion
1 DCPt5 LR2sK Deferred Correction Procedure with profile of the fifth degree polynomial + Line-by-Line Recurrent Method of the second order, accelerated in Krylov subspaces [19]
2 C2 LR2sK Compensatory Transform Technology of the second order accuracy + Line-by-Line Recurrent Method of the second order, accelerated in Krylov subspaces
3 Cl LRlsK Compensatory Transform Technology of the first order accuracy + Line-by-Line Recurrent Method of the first order, accelerated in Krylov subspaces [19]
4 Cl LR1 Compensatory Transform Technology of the first order accuracy + Line-by-Line Recurrent Method of the first order
5 C2 LR2 Compensatory Transform Technology of the second order accuracy + Line-by-Line Recurrent Method of the second order
6 BCGSt9 В Bi-Conjugate Gradient Stabilized Method [17] for nine-diagonal matrix of SLAE with preconditioner on the base of explicit Buleev method [2, 18]
7 DCPt5 BCGSt В Deferred Correction Procedure with profile of the fifth degree polynomial + Bi-Conjugate Gradient Stabilized Method for five-diagonal matrix of SLAE on the base of explicit Buleev method
8 — BSOR9 Block Successive Over Relaxation Method [2] for nine-diagonal matrix of SLAE
The maximum effective value of the iteration parameter was selected for each method in each calculation, since all methods use the iteration parameters. In other words, an upper estimate of the effectiveness was made for each method. This approach made it possible to correctly identify the advantages of one methods in relation to others because all methods were placed in the same conditions.
3.2. Results: coarse and fine meshes
The most interesting for the analysis is the behaviour of the convergence curves which are the dependencies of the \\Rk || / ||i?° || value on the iteration number or the CPU time of the problem. Here ||-Rfc||2 is Euclidean norm of the residual error at the kth iteration. Such convergence curves as functions of the iteration number are plotted in Fig. 4 for the coarse and fine meshes. It is not difficult to see that accuracy of the solution convergence is 10-12. The same value of accuracy is applied in all other results of the work.
Analysis of the curves in Fig. 4a allows the following conclusions. First, the classical block SOR method (curve 8) is not usable due to a huge number of iteration to converge the method -more than one thousand. Naturally, there is not enough place for such curve on the graph. Second, the combination DCPt5 + LR2sK (curve 1) is most effective both in reducing the initial residual
8
J_____I__I___I_1—1—1______________I________L
■ ■ ■ I 100*
10
a) mesh 501x501
10
100
к
b) mesh 4001x4001
Fig. 4. Behavior of the convergence curves depending on the iteration number. Methods: 1 - DCPt5 + LR2sK, 2 - C2 + LR2sK, 3 - Cl + LRlsK, 4 - Cl + LR1, 5 - C2 + LR2, 6 - BCGSt9 B, 7 - DCPt5 + BCGSt B, 8 - BSOR9
error on the first iteration and the total number of iterations for the method convergence. And finally, third, in whole, the versions with line-by-line algorithm are more powerful with respect to the variants of the bi-conjugate gradient method. As to calculations with the fine mesh (see Fig. 4b), here the results coincide qualitatively with the ones on the coarse grid, but, as a rule, the number of iterations for solution convergence is several times greater. The absence of the convergence curve of the BSOR9 method is explained by the lack of convergence of the method - the relative residual error was more than 5 x 10-8 after 3000 iterations.
It is obvious that different methods require different amounts of mathematical operations and, accordingly, different amounts of a CPU time to perform calculations of one iteration. Again, a researcher is ultimately interested in the time spent by a computer for working out a solution. Therefore, a comparison of computation times is also of research interest. Yet it is clear, that only the relative CPU times have a sense here. In other words, only the times of calculations performed on the same computer can be compared. Precisely such results in the form of convergence curves are shown in Fig. 5 for coarse and fine meshes. From now on, CPU time is presented in seconds.
It is easy to see that the combination of DCPt5 + LR2sK methods (curve 1) has not even got into the “top three winners”. The reason is quite clear: recalculation of the right part of the system of linear equations by DC procedure at each iteration is a time-consuming activity. Owing to similar arguments the CPU time of the DCPt5 + LR2sK combination is almost equal to CPU time of the C2 + LR2sK one for the 4001 x 4001 mesh (see Fig. 5b). In all other respects, the relative behaviours of the curves in the graphs of Fig. 4 and Fig. 5 coincide qualitatively.
As is known, the average rate of convergence is one of the most indicative characteristics of the efficiency of an iterative method which is formulated as follows
(13)
Fig. 5. Behavior of the convergence curves depending on CPU time. Methods: 1 - DCPt5 + LR2sK, 2 - C2 + LR2sK, 3 - Cl + LRlsK, 4 - Cl + LR1, 5 - C2 + LR2, 6 - BCGSt9 B, 7 - DCPt5 + BCGSt B, 8 - BSOR9
where \\Zk
is Euclidean norm of the solution error, и — is analytical
solution of the problem. The curves of the average convergence rate for the performed calculations with coarse and fine meshes are shown in Fig. 6. It goes without saying that the higher the curve the more efficient the method is. The low-lying curve 8 once again confirms the relative inefficiency of the classical method BSOR9. The productivities of the other methods are comparable with each other. And as expected from the previous graphs, the highest curve 1 corresponds to the most effective method - the combination DCPt5 + LR2sK. The almost direct behavior of the curves in the logarithmic system of coordinates indicates a power dependence of the average rate of convergence Qk on the iteration number k.
Fig. 6. Behavior of the average convergence rate curves depending on iteration number. Methods: 1 - DCPt5 + LR2sK, 2 - C2 + LR2sK, 3 - Cl + LRlsK, 4 - Cl + LR1, 5 - C2 + LR2, 6 - BCGSt9 B, 7 - DCPt5 + BCGSt B, 8 - BSOR9
Some quantitative results of solving the problem on the coarse and fine grids with the use of the methods under consideration are presented in Tab. 2. The data in brackets for BSOR9 method for grid 4001 x 4001 emphasize the lack of the solution convergence with the required accuracy. Here ||Zfc|| value is an infinity norm of the error at the moment of solution convergence. One
Table 2
The results of solving the problem by various methods with various meshes
Mesh Type of transform Method 9 Number of iterations CPU time, s \\zk\\ II lloo
- BSOR9 1.9870 1 331 73.3 2.85E-05
- BCGSt9 В 0.99980 63 6.2 2.85E-05
DCPt5 BCGSt В 0.999922 92 32.0 2.88E-05
501 x 501 DCPt5 LR2sK 0.99999945 8 3.8 2.88E-05
Cl LR1 0.999720 36 3.1 2.85E-05
C2 LR2 0.99999350 92 7.4 2.85E-05
Cl LRlsK 0.99930 13 2.1 2.85E-05
C2 LR2sK 0.9999950 13 2.1 2.85E-05
- BSOR9 1.9980 (3 000) (29 833.2) (4.13E-03)
- BCGSt9 В 0.9999979 182 1 655.2 4.46E-07
DCPt5 BCGSt В 0.99999942 343 8 425.9 4.82E-06
4001 x 4001 DCPt5 LR2sK 0.9999999995 7 288.0 4.81E-06
Cl LR1 0.999972 204 1 448.7 4.42E-07
C2 LR2 0.999999953 292 2 120.5 4.42E-07
Cl LRlsK 0.999932 38 463.8 4.46E-07
C2 LR2sK 0.999999920 22 279.2 4.45E-07
can see the norm is reduced by two orders of magnitude with a decrease in the value of the grid step by only an order of magnitude. It should also be noted that the number of iterations in this case increases by less than an order of magnitude, while the CPU time is increased by as much as two orders of magnitude on average. Special attention should be paid to the fact that 8 iterations were required for the method convergence on the 501 x 501 grid, and only 7 iterations — on the 4001 x 4001 grid, while using the combination DCPt5 + LR2sK. The explanation for this fact will be presented a little later.
3.3. Influence of the mesh resolution
The influence of the grid resolution on the number of iterations and the CPU time of the solution convergence is presented in Fig. 7. It is known that an increase in the dimension of SLAE (a decrease in the magnitude of the grid step), other things being equal, leads to an increase in the number of iterations [20]. Exactly such regularity takes place in Fig. 7a, except for curve 1 (the combination DCPt5 + LR2sK), which demonstrates the decrease in the number of iterations with increasing the grid dimensionality. The reason is that the original line-by-line recurrent method has a fundamental individuality: as the grid step decreases,
the number of iterations required for convergence decreases [15]. This feature of the method has been manifested only in the combination DCPt5 + LR2sK. Yet other combinations with line-by-line recurrent method (curves 2-5) do not demonstrate the features because the presence
of the approximate compensatory transform technology in the combinations suppresses this fundamental individuality.
Fig. 7. The number of iterations (Kj) and CPU time (T/) required for the convergence of the method, depending on the grid resolution. Methods: 1 — DCPt5 + LR2sK, 2 — C2 + LR2sK, 3 - Cl + LRlsK, 4 - Cl + LR1, 5 - C2 + LR2, 6 - BCGSt9 B, 7 - DCPt5 + BCGSt В
In the general case, line-by-line recurrent method of the second order LR2 is more efficient than the one of the first order LR1, regardless of whether this method was accelerated in Krylov subspaces or not, and the relationship of the curves 2 and 3 confirms this thesis. However, the locations of the curves 4 and 5 demonstrate the opposite. The reason is that additional approximate compensatory transformation technology decreases the method stability. It is necessary to lower the value of compensation parameter в in relation to its optimum value to maintain stability. As an effect, the more parameter в differs from its optimum, the slower the method carries out. And one has to make в lower for LR2 than for LR1 due to lower stability of LR2, which in turn leads to a greater slowdown LR2 in relation to LR1.
And finally, it is not difficult to see that power dependences of the number of iterations Ki and CPU time Tj against a mesh resolution take place because graphs are almost direct in the logarithmic coordinates.
Conclusions
The technology of expanding the use of the line-by-line recurrent method on the case of SLAE with nine-diagonal matrices arising from the difference approximation of 2D boundary value problems of higher order was considered in the article. Approximation of the government differential equation have been carried out using the SMART scheme of the third order of accuracy. Also, approximation of the second order of accuracy in the near-boundary nodes using the Patankar scheme instead of classical upwind one of the first order of accuracy was proposed in the paper. Both the known deferred correction method and the original compensatory transform technology were used in the work to replace an initial nine-diagonal matrix of SLAE with a five-diagonal one. The comparative analysis of several modern methods and their combinations with algorithms for replacement of nine-diagonal matrices with five-diagonal ones to solve SLAE has been performed to reveal their efficiency in relation to each other.
Based on the conducted study, one can draw the following conclusions:
1. The compensatory transform technology does not require recalculation of the right part of the system of equations and, at least, is not inferior in efficiency to the deferred correction method.
2. The high efficiency of the line-by-line recurrent method is also conserved in the solving of systems of linear equations with the nine-diagonal matrix when considering two-dimensional boundary value problems.
3. The number of iterations and, accordingly, the CPU time required for the solution convergence has a power dependence against the grid resolution for all previously explored methods.
This paper is distributed under the terms of the Creative Commons Attribution-Non Commercial 3.0 License which permits non-commercial use, reproduction and distribution of the work without further permission provided the original work is properly cited.
References
1. Patankar S.V. Numerical Heat Transfer and Fluid Flow. Hemisphere Publishing Corporation. New York, 1980. 197 p.
2. Il’in V.P. Iterative Incomplete Factorization Methods. Singapore: World Scientific Publishing Co., 1992. 300 p.
3. Leonard B.P. A Stable and Accurate Convective Modelling Procedure Based on Quadratic Upstream Interpolation. Computer Methods in Applied Mechanics and Engineering. 1979. vol. 19, no. 1. pp. 59-98. DOI: 10.1016/0045-7825(79)90034-3.
4. Gaskell P.H., Lau A.K.C. Curvature-Compensated Convective Transport: SMART, A New Boundedness - Preserving Transport Algorithm. International Journal for Numerical Methods in Fluids. 1988. vol. 8. pp. 617-641. DOI: 10.1002/fld.l650080602.
5. Leonard B.P. The ULTIMATE Conservative Difference Scheme Applied to Unsteady One-Dimensional Advection. Computer Methods in Applied Mechanics and Engineering. 1991. vol. 88, no. 1. pp. 17-74. DOI: 10.1016/0045-7825(91)90232-U.
6. Darwish M.S. A New High-Resolution Scheme Based on the Normalized Variable Formulation.
Numerical Heat Transfer, Part B: Fundamentals. 1993. vol. 24, iss. 3. pp. 353-371.
DOI: 10.1080/10407799308955898.
7. Darwish M.S., Moukalled F. The Normalized Weighting Factor Method: a Novel
Technique for Accelerating the Convergence of High-Resolution Convective Schemes.
Numerical Heat Transfer, Part B: Fundamentals. 1996. vol. 30, iss. 2. pp. 217-237.
DOI: 10.1080/10407799608915080.
8. Chirkov D.V., Chernyi S.G. Comparison of Accuracy and Convergence of Some TVD-Schemes. Vychislitelnye tekhnologii [Computational Technologies]. 2000. vol. 5, no. 5. pp. 86-107. (in Russian).
9. Sukhinov A.I., Chistakov A.E., Yakobovskii M.V. Accuracy of the Numerical Solution of the Equations of Diffusion-Convection Using the Difference Schemes of Second and Fourth Order Approximation Error. Vestnik Yuzhno-Uralskogo gosudarstvennogo universiteta. Seriya: “Vychislitelnaya matematika i informatika” [Bulletin of the South Ural State University. Series:
Computational Mathematics and Software Engineering], 2016. vol. 5, no. 1. pp. 47-62. (in Russian) DOI: 10.14529/cmsel60105.
10. Prokudina L.A., Yaparova N.M., Vikhirev M.P Numerical Simulation of the Oscillations of
the Elements of the Pipe with the Flow of an Incompressible Fluid. Vestnik Yuzhno-JJralskogo gosudarstvennogo universiteta. Seriya: “Vychislitelnaya matematika i informatika” [Bulletin of the South Ural State University. Series: Computational Mathematics and Software
Engineering], 2018. vol. 7, no. 3. pp. 55-64. (in Russian) DOI: 10.14529/cmsel80304.
11. Schneider G.E., Zedan M. A Modified Strongly Implicit Procedure for the Numerical Solution of Field Problems. Numerical Heat Transfer. 1981. vol. 4, no. 1. pp. 1-19. DOI: 10.1080/01495728108961775.
12. Zverev V.G. Modified Line-by-Line Method for Difference Elliptic Equations. Zhurnal vychislitelnoi matematiki i matematicheskoi fiziki [Computational Mathematics and Mathematical Physics]. 1998. vol. 38, no. 9. pp. 1490-1498. (in Russian).
13. Sikovskii D.F. Metody vychislitelnoi teplofiziki [Computational Thermal Physics Methods]. Novosibirsk: NSU, 2013. 98 p. (in Russian).
14. Fomin A.A., Fomina L.N. On the Solution of Fluid Flow and Heat Transfer Problem in a 2D Channel with Backward-Facing Step. Vestnik Samarskogo gosudarstvennogo tekhnicheskogo universiteta. Seriya: “Fiziko-matematicheskie nauki” [Journal of Samara State Technical University. Series: Physical and Mathematical Sciences]. 2017. vol. 21, no. 2. pp. 362-375. DOI: 10.14498/vsgtul545.
15. Fomin A.A., Fomina L.N. On the Convergence of the Implicit Iterative Line-by-Line Recurrence Method for Solving Difference Elliptical Equations. Kompyuternye issledovaniya i modelirovanie [Computer Research and Modeling], 2017. vol. 9, no. 6. pp. 857-880. (in Russian) DOI: 10.20537/2076-7633-2017-9-6-857-880.
16. Darwish M.S., Moukalled F.H. Normalized Variable and Space Formulation Methodology for High-Resolution Schemes. Numerical Heat Transfer, Part B: Fundamentals. 1994. vol. 26, iss. 1. pp. 79-96. DOI: 10.1080/10407799408914918.
17. Van der Vorst H.A. BI-CGSTAB: a Fast and Smoothly Converging Variant of BI-CG for the Solution of Nonsymmetric Linear Systems. SIAM Journal on Scientific and Statistical Computing. 1992. vol. 13, iss. 2. pp. 631-644. DOI: 10.1137/0913035.
18. Starchenko A.V. Comparative Analysis of Some Iterative Methods for the Numerical Solution of a Spatial Boundary Value Problem for Elliptic Equations. Vestnik Tomskogo gosudarstvennogo universiteta. Byulleten operativnoi nauchnoi informatsii [Bulletin of the Tomsk State University. The Bulletin of Operational Scientific Information], Tomsk: TSU, 2003. no. 10. pp. 70-80. (in Russian).
19. Fomin A.A., Fomina L.N. Acceleration of the Line-by-Line Recurrent Method in Krylov Subspaces. Vestnik Tomskogo gosudarstvennogo universiteta. Matematika i mekhanika [Tomsk State University Journal of Mathematics and Mechanics]. 2011. no. 2. pp. 45-54. (in Russian).
20. Faddeeva V.N. Computational Methods of Linear Algebra. N.Y.: Dover Publications, 1959. 252 p.
УДК 519.63 DOI: 10.14529/cmsel90201
ПРИМЕНЕНИЕ НЕЯВНОГО ИТЕРАЦИОННОГО ПОЛИНЕЙНОГО РЕКУРРЕНТНОГО МЕТОДА ПРИ РЕШЕНИИ СИСТЕМ РАЗНОСТНЫХ ЭЛЛИПТИЧЕСКИХ УРАВНЕНИЙ С ДЕВЯТИДИАГОНАЛЬНЫМИ МАТРИЦАМИ
© 2019 А.А. Фомин1, Л.Н. Фомина2
1 Кузбасский государственный технический университет имени Т.Ф. Горбачева (650000 Кемерово, ул. Весенняя, д. 28),
2 Кемеровский государственный университет (650043 Кемерово, ул. Красная, д. 6)
E-mail: fomin_aa@mail.ru, lubafomina@mail.ru Поступила в редакцию: 21.11.2018
В статье исследуется применение неявного итерационного полинейного рекуррентного метода для решения систем линейных разностных уравнений с девятидиагональными матрицами, которые возникают при разностной аппроксимации двумерных задач повышенного порядка точности на регулярном сеточном покрытии области решения. Поскольку изначально неявный итерационный полинейный рекуррентный метод разработан для решения систем уравнений с пятидиагональной матрицей, в работе предлагается технология так называемой компенсационной трансформации, позволяющая заменить исходную девятидиагональную матрицу системы уравнений на пятидиагональную. Эффективность подобного подхода анализируется путем сравнения параметров сходимости решения модельной краевой задачи в единичном квадрате как различными вариантами предложенного метода, так и другими современными высокоэффективными итерационными методами решения систем разностных уравнений. Задача решается на последовательности сеток от грубой в 501 х 501 узлов до подробной в 4001 х 4001 узлов. Точность сходимости решения определяется по относительной норме невязки, которая в настоящей работе равняется 10-12. Показано, что несмотря на использование промежуточной технологии компенсационной трансформации, неявный итерационный полинейный рекуррентный метод сохраняет свои высокие скоростные и разрешающие способности во всем диапазоне сеточного разбиения области решения задачи.
Ключевые слова: метод сеток, система разностных эллиптических уравнений, итерационный метод, сходимость решения.
ОБРАЗЕЦ ЦИТИРОВАНИЯ
Fomin A. A., Fomina L. N. The Use of Line-by-Line Recurrent Method for Solving Systems of Difference Elliptic Equations with Nine-Diagonal Matrices / / Вестник ЮУрГУ. Серия: Вычислительная математика и информатика. 2019. Т. 8, № 2. С. 5-21.
DOI: 10.14529/cmsel90201.
Литература
1. Патанкар С. Численные методы решения задач теплообмена и динамики жидкости. М.: Энергоатомиздат, 1984. 152 с.
2. Ильин В.П. Методы неполной факторизации для решения алгебраических систем. М.: Физматлит, 1995. 288 с.
3. Leonard В.Р. A Stable and Accurate Convective Modelling Procedure Based on Quadratic Upstream Interpolation // Computer Methods in Applied Mechanics and Engineering. 1979. Vol. 19, No. 1. P. 59-98. DOI: 10.1016/0045-7825(79)90034-3.
4. Gaskell P.H., Lau A.K.C. Curvature-Compensated Convective Transport: SMART, A New Boundedness - Preserving Transport Algorithm // International Journal for Numerical Methods in Fluids. 1988. Vol. 8. P. 617-641. DOI: 10.1002/fld.l650080602.
5. Leonard B.P. The ULTIMATE Conservative Difference Scheme Applied to Unsteady One-Dimensional Advection // Computer Methods in Applied Mechanics and Engineering. 1991. Vol. 88, No. 1. P. 17-74. DOI: 10.1016/0045-7825(91)90232-U.
6. Darwish M.S. A New High-Resolution Scheme Based on the Normalized Variable Formulation // Numerical Heat Transfer, Part B: Fundamentals. 1993. Vol. 24, Iss. 3. P. 353-371. DOI: 10.1080/10407799308955898.
7. Darwish M.S., Moukalled F. The Normalized Weighting Factor Method: a Novel
Technique for Accelerating the Convergence of High-Resolution Convective Schemes // Numerical Heat Transfer, Part B: Fundamentals. 1996. Vol. 30, Iss. 2. P. 217-237. DOI: 10.1080/10407799608915080.
8. Чирков Д.В., Черный С.Г. Сравнение точности и сходимости некоторых ТVD-схем // Вычислительные технологии. 2000. Т. 5, № 5. С. 86-107.
9. Сухинов А.И., Чистяков А.Е., Якобовский М.В. Точность численного решения уравнения диффузии-конвекции на основе разностных схем второго и четвертого порядков погрешности аппроксимации // Вестник ЮУрГУ. Серия: Вычислительная математика и информатика. 2016. Т. 5, № 1. С. 47-62. DOI: 10.14529/cmsel60105.
10. Прокудина Л.А., Япарова Н.М., Вихирев М.П. Численное моделирование колебаний элементов трубы с потоком несжимаемой жидкости / / Вестник ЮУрГУ. Серия: Вычислительная математика и информатика. 2018. Т. 7, № 3. С. 55-64. DOI: 10.14529/cmsel80304.
11. Schneider G.E., Zedan М. A Modified Strongly Implicit Procedure for the Numerical Solution of Field Problems // Numerical Heat Transfer. 1981. Vol. 4, No. 1. P. 1-19. DOI: 10.1080/01495728108961775.
12. Зверев В.Г. Модифицированный полинейный метод решения разностных эллиптических уравнений // Журнал вычислительной математики и математической физики. 1998. Т. 38, № 9. С. 1553-1562.
13. Сиковский Д.Ф. Методы вычислительной теплофизики: Учеб, пособие. Новосибирск: Новосибирский государственный университет, 2013. 98 с.
14. Fomin A.A., Fomina L.N. On the Solution of Fluid Flow and Heat Transfer Problem in a 2D Channel with Backward-Facing Step // Вестник Самарского государственного технического университета. Серия: Физико-математические науки. 2017. Т. 21, № 2. С. 362-375. DOI: 10.14498/vsgtul545.
15. Фомин А.А., Фомина Л.Н. О сходимости неявного итерационного полинейного рекуррентного метода решения систем разностных эллиптических уравнений / / Компьютерные исследования и моделирование. 2017. Т. 9, № 6. С. 857-880. DOI: 10.20537/2076-7633-2017-9-6-857-880.
16. Darwish M.S., Moukalled F.H. Normalized Variable and Space Formulation Methodology for High-Resolution Schemes // Numerical Heat Transfer, Part B: Fundamentals. 1994. Vol. 26, Iss. 1. P. 79-96. DOI: 10.1080/10407799408914918.
17. Van der Vorst H.A. BI-CGSTAB: a Fast and Smoothly Converging Variant of BI-CG for the Solution of Nonsymmetric Linear Systems // SIAM Journal on Scientific and Statistical Computing. 1992. Vol. 13, Iss. 2. P. 631-644. DOI: 10.1137/0913035.
18. Старченко А.В. Сравнительный анализ некоторых итерационных методов для численного решения пространственной краевой задачи для уравнений эллиптического типа // Вестник ТГУ. Бюллетень оперативной научной информации. Томск: ТГУ, 2003. № 10. С. 70-80.
19. Фомин А.А., Фомина Л.Н. Ускорение полинейного рекуррентного метода в подпространствах Крылова / / Вестник Томского государственного университета. Математика и механика. 2011. № 2. С. 45-54.
20. Фаддеев Д.К., Фаддеева В.Н. Вычислительные методы линейной алгебры. М.: Физматгиз, 1963. 656 с.
Фомин Александр Аркадьевич, к.ф.-м.н., ст.н.с., отдел развития и международного сотрудничества, Кузбасский государственный технический университет имени Т.Ф. Горбачева (Кемерово, Российская Федерация)
Фомина Любовь Николаевна, к.ф.-м.н., доцент, кафедра ЮНЕСКО по информационным вычислительным технологиям, институт фундаментальных наук, Кемеровский государственный университет (Кемерово, Российская Федерация)
