CC BY
8
Comparative analysis of Jacobi and Gauss-Seidel
iterative methods
P. V. Khrapov, N. S. Volkov
Abstract — The paper presents a comparative analysis of iterative numerical methods of Jacobi and Gauss-Seidel for solving systems of linear algebraic equations (SLAEs) with complex and real matrices. The ranges of convergence for both methods for SLAEs in two and three unknowns, as well as the interrelationships of these ranges are obtained. An algorithm for determining the convergence of methods for SLAEs using the complex analog of the Hurwitz criterion is constructed, the realization of this algorithm in Python in the case of SLAEs in three unknowns is given. A statistical comparison of the convergence of both methods for SLAEs with a real matrices and the number of unknowns from two to five is carried out.
Keywords — iterative methods, system of linear algebraic equations, Jacobi method, Gauss-Seidel method, stable polynomials, Hurwitz criterion.
I. Introduction
In the modern world, a large number of both applied and theoretical problems in various fields of science and technologies are reduced to the problem of finding exact solutions of various SLAEs or solutions that maximally approximate the exact ones, numerical methods for solving which have been developing over the years due to the huge number of areas of their application [1], [2].
A special place in the theory of SLAEs' solutions is occupied by the simple iterative method, which is an alternative to direct methods of finding SLAEs' solutions. At the same time, based on the simple iterative method, new methods for solving SLAEs are developed, which are an improved version of the classical method [3], [4], [5].
Some of these, based on the simple iterative method, are the Jacobi and Gauss-Seidel iterative methods for solving SLAEs, the meaning of which is to allocate elements on, above and below the diagonal of the original SLAE's matrix as separate matrices and conduct the simple iterative method using them instead of the original, which often greatly simplifies the calculations [6], [7]. Iterative Jacobi and
Manuscript received July 19, 2023.
P. V. Khrapov - Bauman Moscow State Technical University (5/1 2-nd
Baumanskaya St., Moscow 105005, Russia),
ORCID: https://orcid.org/0000-0002-6269-0727 ,
e-mail: pvkhrapov@gmail.com , khrapov@bmstu.ru.
N. S. Volkov - Bauman Moscow State Technical University (5/1 2-nd
Baumanskaya St., Moscow 105005, Russia),
ORCID: https://orcid.org/0009-0008-7095-4092 ,
e-mail: nikita.volkov01@mail.ru, volkovns@student.bmstu.ru.
Gauss-Seidel methods, also being classical iterative methods of solving SLAEs, have recently undergone various improvements, some of which are described for example in [7], [10], [12], [13], [19], [20], [21], [22], [25]. Nevertheless, many modern alternatives to the classical Jacobi and Gauss-Seidel methods are based on sufficient condition of their convergence to an exact solution in the case of diagonal pre- dominance in the original SLAEs' matrices, without considering the cases without diagonal predominance when these methods can also converge to an exact solution, and are also described only for special types of matrices [23], [24].
The convergence of iterations to an exact solution is one of the main problems, since, as a consequence of the classical simple iterative method, the Jacobi and Gauss-Seidel methods not always converge to an exact solution, and have convergence criteria following from a similar criterion for the simple iterative method [6]. The search for convergence ranges and the theoretical comparison of the effectiveness of the methods based on it is the main task of this work.
The convergence criteria obtained in [6], according to which the eigenvalues of the matrices in the method should be less than one in absolute value, are reduced to the problem of finding the roots of the algebraic polynomials of degree n with complex coefficients inside the unit circle, various solutions of which are described for example in [8], [9], [11], [17], [27], [28], [29], [30], [31], and for polynomials of a special kind in [16], [18].
It can be solved by making a fractional linear transformation that translates the interior of the unit circle of the complex plane to the left half-plane and reduces it to the study of stability of the polynomial [11]. In [11] this problem is considered for polynomials with real coefficients of the second and third degree.
In this paper, a comparative analysis of two methods using the examples of SLAEs in two and three unknowns is carried out by considering the ranges of their convergence, which are obtained under the assumption that the boundary of each range is formed when at least one root of the corresponding equation has a unit absolute value, and the rest does not exceed one, and all points are contained inside the range, for which all roots have absolute values less than one.
There is described the general convergence criteria for each method in paragraph 2.
In paragraph 3, the convergence ranges of the methods for SLAEs with complex coefficients in two unknowns are obtained, and the conclusion of their comparison is given: for the Jacobi method the convergence range and its boundary are found when substituting roots with absolute values less
than one into the corresponding equation, and for the Gauss-Seidel method by directly solving the equation.
In paragraph 4, by a similar substitution of roots with absolute values not exceeding one, the boundary conditions of the methods for SLAEs with complex coefficients in three unknowns are obtained, and on their basis the convergence ranges in the real case are obtained, for which a comparative analysis is given.
There is described a general convergence check method for SLAEs with complex matrices based on [11] and [14], and a general comparison of both methods is made in paragraph 5.
In paragraph 6, a statistical comparison of convergence of both methods for SLAE with real matrix is carried out using mathematical modeling.
II. Convergence conditions of the Jacobi and Gauss-Seidel methods
When solving a system of linear algebraic equations
Ax = b (1)
in accordance with the Jacobi method, the matrix A of the original SLAE is represented as a sum: A = L + D + R det A * 0
where L, D, R are, respectively, the matrices with subdiagonal, diagonal, and overdiagonal elements of matrix A, and then there is a system obtained from the original SLAE (1):
x = - D-'(L + R) x + D-lb for which the simple iterative method converges if all roots of the equation
Xau
a21
a
a
12
a
1 n
Xa,,
a
a
Xa„
= 0
(2)
n1 n 2
have absolute values less than one [6]; aij - elements of the
original matrix A, a. e □ .
Similarly, the Gauss-Seidel method transforms the original SLAE (1) to a system:
X = -(L + D )-1 Rx + (L + D )-1 b for which the simple iterative method converges to an exact solution if all roots of equation
XaY Xan
a12 Xa.
22
a.
a,,
Xan1 Xan2 ■ " Xan
= 0
(3)
have absolute values less than one [6], a.. e □ .
i.
When the dimension of the original SLAE is small, we can find the convergence ranges of the methods by directly solving the equations (2) and (3). Let us show this for the
cases of SLAEs in two and three unknowns, which often arise in applied research.
III. SYSTEM OF LINEAR ALGEBRAIC EQUATIONS IN TWO UNKNOWNS
A. Jacobi method The equation (2) has the form:
Xau
a,
12
a,,
Xa..
= X a11a22 - a12a21 = 0
(4)
and for the convergence of the method, it is necessary that its roots lie inside the unit circle.
In the general case a., Xl2 e □ , and the system (2) is equivalent to equations
(X- rxeipp )(X - r2 eP) = 0
X2 +X( - re - r2 eP) + rre ei<P2 = 0 (5) where r1e"Pi, r2e1<P2 are the roots of equation (4), rx, r2 < 1.
Comparing (4) and (5), we obtain the system (ana22 * 0, since an, a22 are elements of the diagonal matrix D):
rè + r2el<P2 = 0
iœ iœ a12a21
r;r2e ^e =—
(6)
a11a22
from which follow:
r1 = r2
r1r2 = r1 =
a12 a21
a11a22
< 1
a12 a21 < a11a22
(7)
The condition (7) defines the convergence range for the Jacobi method in the general case: the absolute value of the product of the off-diagonal elements of the matrix A of the system (1) must be less than the absolute value of the product of its diagonal elements for the method to converge in the case of an SLAE in two unknowns.
B. Gauss-Seidel method The equation (3) has the form:
Xa,
a,
Xa21 Xa,
'22
and its roots
= Xana22 -Xa12a21 = 0 X = 0
n n
X =■
(8)
a11a22
must have an absolute value less than one ( a11a22 ^ 0, since a11, a22 are diagonal elements of triangular matrix L + D ).
Since one of them is zero, only the second root is checked for the convergence condition, for which, in order for its absolute value to be less than one, it is necessary to fulfill the condition
|а12a2l| < \aua22\ (9)
Thus, both the Jacobi method and the Gauss-Seidel method for SLAEs in two unknowns have the same range of convergence (9).
IV. System of linear algebraic equations in three
UNKNOWNS
A. Jacobi method The equation (2) has the form
Я aiia22a33 A(ai3a22a3i + a23aiia32 + ai2a33a21 ) +
+ (ai3a32a2i + ai2a23a3i) = 0 (i0)
and for convergence of the method, it is necessary that all its roots lie inside the unit circle.
Let's divide it by aiia22 a33 (there are no zero elements on the diagonal of matrix A, since matrix D must have an inverse):
Л3 +
^13^22^31 a23aiia32 a\2a33a2\
Л +
a\\a22a33 a\3a32a2\ + a\2 a23a3\
a\\a22a33
= 0
Denoting
P =
-a\3a22a3\ — a23a\\a32 — a\2a33a2\ a\\a22a33 a\3a32a2\ + a\2 a23a3\
q =
a\\a22a33
+r2 r3em em) - r2 r3em e
m —
= 0
r r < 1
'2' '3 — 1
Comparing (11) and (12), we obtain a system of equations for the first boundary of the convergence range:
+ r2em + r3eim = 0
г2е'и e'w + r3em e3 + r2 re e= p (i3) -r2 r3em em em = q from which follow:
|q| < \
arg (q) =n+m +m+m
(\4) (\5)
r3e1<p3 = -r2e1<p2 - e1(p1 (16)
The first equation of the system (13) has a geometric interpretation (fig. 1)
JmX
ReX
we obtain the canonical cubic equation:
A3 + pA + q = 0 (11)
In the general case its coefficients and roots are complex:
p, q, A1,2,3 e □ .
We find the convergence range of the method expressed in terms of p, q e □ by obtaining the equations of its boundaries and combining them.
We obtain the equations of the boundaries under the assumption that there is at least one root of the equation (11) on boundaries, the absolute value of which is equal to one, and the interior points of the range are those in which the absolute value of each root is less than one. The boundary is not included in the convergence range, since at least one of the roots has a unit absolute value on it, which contradicts the convergence condition [6]. Consider several cases.
Find the equation for the first boundary of the convergence range: let one of the roots of the equation (11) have a unit absolute value, and the other two roots have an absolute value not exceeding one:
(A- e9)(A- r2 e92)(A- r3 e9) = 0 (12) A3 + A2(-e* - re - r3e'9) + A(r2e9e9 + r3el<P1 e9 +
Fig. 1. Geometric interpretation of the first equation of the system (13).
Substituting (16) into the second and third equations of the system (13), we have expressions:
p = - e m e to - ^ e m e to - ^e to e to
q = r2 em e'9 em + r22em e'9 e'9
comparing which, we obtain the equation of the first boundary of the convergence range of the Jacobi method in the general case:
p = -qe-91 - e2m (17)
from which follow the relationships of absolute value and argument for p and q :
|p| = Vr2 + 1 + 2rq cos(9q - 39,) (18)
, , , Jq sin(9q -91) + sin(291)
arg(p) = arctan(—---) (19)
rq cos(9q -9) + cos(29) The relationships (18) and (19) show that the absolute value and argument forp depend on three parameters - the
absolute value rq, the argument 9, and the argument 9q, the last of which depends not only on 9 (15), so we will take the argument 9 as a parameter to visualize the absolute value and argument for p .
Let us take for example 9 = 0 and 9 =n (we use these parameters for further visualization of the special case of
SLAE with real matrix when at least one of the roots of the equation (11) is real): at p1 = 0, according to geometrical considerations (fig. 1) and (15)
(pq = œ
n
-n-2
n 2 - ^
p = -q -1
Similarly, when ( =n
p = p2 + p3 + n + n e
n n
(20) (21)
(22) (23)
p = q -1
The picture of the absolute value (18) in this case is as follows (fig. 2):
.ReX
Fig. 3. Location of roots of the equation (11) on the complex plane before rotation.
Fig. 4. Location of roots of the equation (11) on the complex plane after rotation.
(X- r1eipp )(X - el (p2+p))(X-<
,l'(-(2 +()\ —
) = 0 (24)
r < 1
Fig. 2. Dependence of the absolute value p| at the boundary (17) at p = 0 (Re(q) < 0) and p =n (Re(q) >
0).
To find the second boundary of the convergence range, consider the case when two roots of the equation (11) on the complex plane have axial symmetry with respect to the line passing through the vector of the third root (the case when two roots are complex-conjugate and the third is real is a special case of this case), and the roots located symmetrically have a unit absolute value, and the third root has an absolute value not exceeding one.
This case can be considered as a rotation of the system of vectors of roots of the equation on the complex plane from the zero angle by the angle p1, which is the argument of the first root: taking into account that before the rotation by the angle p one root was real, and the other two roots were
complex-conjugate with arguments p2 and -p2
respectively, after the rotation the picture on the complex plane will be as follows (fig. 3, fig. 4):
!m\
Opening the brackets and comparing (24) with equation (11), we obtain the system for the second boundary of the convergence range:
' ei p2 + P\) + ei (-P2 + P) + reP = 0
elip1 + r1ei (p+p2)ei<p1 + r/(p ep = p (25) y-r1eii(p1 = q from which follow:
\q\ ^ 1 (26)
•/p =-e (P2+P1) - e(-Pi+P1) (27)
Let's substitute (27) into the second and third equations of the system (25):
,2iP _ e2i(P P2) (28)
p = -e'2i((+(2) -e2i( -e2i((-(2)
q = e2iv (ei(Pi+P2) + ei-Pi))
qe~2iP = ei(P+P2) + ei(Pi (29)
Comparing the square of the expression (29) and the expression (28), we obtain the equation of the second boundary of the convergence range of the Jacobi method in general case:
p = —q 2e—ipp + e2ipp (30)
arg(q) = 3pi +n (3i)
from which we find the relationships of absolute value and argument for p and q subject to the condition (3 i):
|p| = I — r12 = i — |q|2 (32)
2 2 arg(p) = 2p = 3 arg( q) — 3 n (33)
When p1 G [—n;n],the picture for the absolute value (32) is as follows (fig. 5):
In the real case, when aij e A; p, q are real numbers,
Fig. 5. Dependence of the absolute value p| on q.
It follows from (31), (33) that when 9 = ±n, arguments arg(q) = arg( p) = 0 are in the same phase corresponding to the real case, and when 9 = 0
arg(q) = n arg(p) = 0
the arguments are in antiphase corresponding also to the real case (p is positive, q is negative).
In addition, for the first bound (17), it follows from (19)
that for 9q e {0, ±n} and 9 e {0, ±n} as in the case of
the second bound:
arg(p) e{0, ±n}
Hence, when at least one of the roots of the equation (11) is real, the cross section of the convergence range boundaries shown in the figures (fig. 2, fig. 5), at Im(q) = 0, can be combined into one general boundary of the convergence range of the Jacobi method in the real case (fig. 6) (by getting rid of absolute values, part of the boundaries, according to (21) and (23), moves to the area of negative values of p ):
(34)
and the roots of equation (11) are either all real or two of them are complex-conjugate, equations (17) and (30) form the following boundary of the convergence range (fig. 6):
p = -q -1 p = q -1
p = -q2 +1
-1 < q < 1
The boundary, due to the above assumptions that it contains at least one root of the equation (11) whose absolute value is equal to one, does not belong to the convergence range. The convergence range consists of the set of points bounded by the boundary (34), which does not belong to this range, so for the convergence range, given the conditions (14) and (26), it follows that:
-1 < q < 1
Let's show that the area in the figure (fig. 6) can be filled completely:
p =
a13a31 a11a33
a23a32 a22 a33
a12 a21 a11a22
a13a32 a21a22 a31
a11a22a33a21a32
a12 a23a31a11a32
a11a22 a33a12a31
a12a21
a11a22
Let's denote:
a,
x =
22
a
21
y = ■
a
a
12
t =
a
31
a
p =
a13a32a21 a11a22a33
xt -
32
a12 a23a31 y a11a22 a33 t
J_ xy
Assuming that:
Fig. 6. Boundary of the convergence range of the Jacobi method for the real case.
In particular, when 9 e {0, ±n}, the equation (30) takes the form of a parabola:
p = q2+1
Thus, the boundary of the convergence range of the Jacobi method in the general complex case is the union of the sets of points satisfying the equations (17) and (30) and the condition:
Iql < 1
y
xt = — = a t
t * 0
xy = a2
1
p = - aq —2
a
(35)
(36)
x, y, t are independent of each other, so it is always possible to choose the coefficients of the matrix A of the system (1) such that the condition (35) is satisfied and the line (36) is obtained. At the same time, p and q depend on three more parameters on which x, y, t do not depend, so it is possible to choose p and q such that they lie in the convergence range. Thus, we can construct any number of lines of the form (36), some of whose points lie inside the convergence range labeled in the figure (Fig. 6). The set of such lines completely intersects the convergence range. Accordingly, it is always possible to find a SLAE for which
p and q lie within the convergence range of the Jacobi method.
B. Gauss-Seidel method
The equation (3) has the form:
X a11a22a33 +X (a21a13a32 - a13a22a31 -
-a32ana23 -a^a33a12) + Xa^a23a31 = 0 (37) and the method converges if all its roots lie inside the unit circle.
One of the roots of the equation (37) is zero, and the other two roots are found from the quadratic equation:
X2 a + Xd + b = 0 (38)
where
d = a21a13a32 - a13a22a31 - a32a11a23 - a21a33a12 a = ana22 a33
b = ana23a31
For the Gauss-Seidel method, we find the convergence range, given that within it all roots of the equation (38) have an absolute value not exceeding one.
In the general case a, b, d, X1,2 e □ and the equation (38) is equivalent to equation:
aX2 - Xa(re + r2ep) + arlr2e,P1 e,tp2 = 0 0 < r1, r2 < 1 comparing it with (38), obtain the system:
-ar,em - are
i(2 —
2
( = b
= d
(39)
I ar1r2e " e
which defines the convergence range of the Gauss-Seidel method in the general case, and from which it follows that in
the convergence range ( a * 0, since aii are the diagonal elements of the triangular matrix L+D, i e {1,2,3}):
b
a
< 1
(40)
and the boundary of the convergence range (39), assuming that at least one of the roots of the equation (38) has a unit absolute value on it (let r1 = I ), is given by the conditions:
d = -aeip - be-*1
ar2eipp e( = b 0 < r2 < 1
(41)
Pb, =P +(2
|dJ = 41 + rh + 2rb CoS(( -(2)
arg(dj) = arctan(-
sin (j + rb sin (2
cos p1 + rb cos p2 In particular, when the matrix A of the system (i) contains real coefficients, a, b, d G □ , solving directly the quadratic equation (38) and applying the convergence criterion of the Gauss-Seidel method, taking into account the condition (40), we obtain the convergence range of the Gauss-Seidel method in the case of real roots of the equation (38):
Id < la + b|
< 1
(44)
and in the case of complex-conjugate roots of the equation (38):
0 < — < 1 (45)
a
Note that in the latter case, the condition
|d| < |a + b|
follows directly from the condition (45) and the negativity of the discriminant of the equation (38): d2 < 4ab. Therefore, the system (44) is a single range of convergence of the Gauss-Seidel method in the case of real matrix elements of the system (1).
The first condition of the system (44) is interpreted as a segment d on an infinite line.
Note that the conditions (44) are consistent with the boundary (41).
Unlike the Jacobi method, the convergence range of the Gauss-Seidel method in the case of real coefficients of the equation (38) is not constant, and the length of the above segment can vary depending on the parameters a and b.
Let's compare the convergence ranges of both methods in the case of real coefficients of the system (1). For this purpose, we construct the convergence range bounded by the boundary (34) and the range (44) on the same coordinate plane qOp. The parameters p and q for the Jacobi method and d, a, b for the Gauss-Seidel method are related by the relation:
d = (p + q)a - b substituting it into (44), we obtain:
\(p + q)a- b| < la + b|
, d b For dj = — and bj = — on the boundary (41), we can a a
also find the relationships between absolute values and arguments:
d =-e( - bxe~i(h d =- cos (1 - rh cos( -(1) - i(slnPl + rh sin((b1 -(1))
< 1
(46)
Expanding the absolute values in the first inequality of the system (46), we find that one of the boundaries of the convergence range of the Gauss-Seidel method is always the line
p = -q -1 (47)
which is also one of the boundaries (34) of the convergence range of the Jacobi method, and the second one is also a straight line, which has the following form:
p = - q+-
a + 2b
a
It also shows that the size of the convergence range of the Gauss-Seidel method depends on the parameters a, b. Moreover, at some values of these parameters the convergence range of the Gauss-Seidel method can partially pass through the convergence range of the Jacobi method, and at other values it can completely contain it.
From the second inequality of the system (46) follows:
, a + 2b „
-1 <-< 3
a
so, the convergence range of the Gauss-Seidel method on the plane qOp is a part of this plane, which is always bounded from below by the line (47), and, depending on the particular case, bounded from above by a line parallel to it, the uppermost of which is the line
p = -q+3
Thus, together the convergence ranges of each method on the same plane qOp are as follows (fig. 7):
Fig. 7. Convergence ranges of methods on the plane qOp. The band a is the maximum (with upper boundary
p = -q + 3 ) convergence range of the Gauss-Seidel method; the area ABCD is the convergence range of the Jacobi method.
According to fig. 7, the advantages of the Gauss-Seidel method over the Jacobi method when the system (1) has real matrix elements are obvious (in the case in fig. 7, the convergence range of the Jacobi method is entirely contained in the convergence range of the Gauss-Seidel method), especially when the parameters p and q have large absolute values - then the Jacobi method does not converge. Nevertheless, the upper bound of the range for the Gauss-Seidel method varies depending on the parameters a and b, so if the iterative process of the Jacobi method converges for the SLAE, it does not mean that the iterative process of the Gauss-Seidel method converges.
Let's give examples of constructing the convergence range of the Gauss-Seidel method in coordinates qOp to demonstrate how it varies depending on the parameters a and b, and in the same coordinates we construct the convergence range of the Jacobi method for clarity.
Example 1.
Let the parameters a = 2, b = 1, then the convergence range of the Gauss-Seidel method has the form: '| 2( p + q) -1 < 3
= 1 < 1
2
thus:
a + 2b
= 2
a
Then, by analogy with fig. 7, the convergence ranges for each method on the plane qOp look as follows (fig. 8):
-1.5 -1 -Ö.S
O.S 1 1.5
Fig. 8. Convergence ranges of Jacobi and Gauss-Seidel methods at parameters a = 2 and b = 1.
The figure 8 shows that the convergence range of the Jacobi method lies entirely within the convergence range of the Gauss-Seidel method, so in this particular case of parameters a, b for any SLAE for which the Jacobi method converges, the Gauss-Seidel method also converges, but the converse is not true.
Example 2.
Here is an example of a SLAE in three unknowns, for which the Jacobi method converges, but the Gauss-Seidel method does not converge:
C-8 6 -4 ^
A =
-9 8 4 -5
6 3
v ■ " •'J In this case the parameters are as follows: a = -192 b = 144
then the convergence range of the Gauss-Seidel method has the following form:
1-192(p + q) -144 < 48
< 1
< b 144
a 192
a + 2b
a
1 2
By analogy with fig. 7, we obtain the following picture of convergence ranges for both methods on the plane qOp
(fig 9):
«gtnct lanzr of
Fig. 9. Convergence ranges of Jacobi and Gauss-Seidel methods at parameters a = -192 and b = 144.
This range does not satisfy the parameters p, q, which in this particular example for matrix A are equal to:
50
P =
q _
192
36
192
From the figure 9 it is obvious that this point (p, q) does not belong to the convergence range of the Gauss-Seidel method on the plane qOp, but it belongs to the convergence range of the Jacobi method.
In addition, in this particular case we see that the Gauss-Seidel method does not converge in most of the convergence range of the Jacobi method, but it can converge at large values of p and q, while the Jacobi method does not converge at large values of p and q.
V. The general case of systems of linear algebraic
EQUATIONS WITH COMPLEX MATRICES
The convergence check of each method is an investigation to find all roots of a polynomial of degree n inside the unit circle, which can be transformed to a stability study problem [11].
In general, a polynomial of degree n with complex coefficients is obtained from the determinant equations (2) or (3):
f (A) = a0An + a1An-1 +... + an = 0, a0 * 0 (48)
For convergence of the method to which the given polynomial corresponds, it is necessary and sufficient that all its roots lie inside the unit circle, for which, in turn, it is necessary and sufficient that the polynomial:
f (z) = a0(z + 1)n + a,(z + 1)n-1(z-1) + a2(z + 1)"-2(z-1)2 +...
+an (z -1)n = 0
obtained from (48) be stable [11].
In general, to check its stability, we can use the complex analog of Hurwitz's stability criterion [14]: let there be an
arbitrary polynomial of degree n with complex coefficients, the stability of which should be investigated:
f ( z ) = d0zn + dizn-1 +... + dn it's equivalent, under the assumption that d0 ^ 0, to the polynomial whose first coefficient is equal to one:
f ( z ) = zn + d1 zn-1 +... + d- (49) d 0 d 0 Replacing in (49) the variable z by a purely imaginary number io, o e □ , we have the polynomial:
/(70) = (io)n + d(ia)n-1 +... + ^ d0 d0 which, by raising the multiplier io of each summand to the appropriate degree and separating the purely imaginary elements from the purely real ones, is represented as the sum of two polynomials with real coefficients
f (io) = g (o) + ih(o) for which, according to [14], if the degree of the polynomial (49) is n = 2m :
g = (-1) mg h = (-i)m-1 h
if n = 2m +1 :
g = (-1)mh h = (-1) mg
Let
B = b0 xn + bxn-1 +... + b
0 1 n
be an arbitrary polynomial of degree n with real coefficients with positive prime factor b0 , and let
C_c xn
-c1 x
■'n-1
be an arbitrary polynomial of degree at most n -1 with real coefficients. Definition. Square matrix
f b0 b1 b2 b3 • • 0 >
0 c0 c1 c2 • ■■ 0
0 b0 b1 b2 • ■■ 0
0 0 c0 c1 • 0
0 0
b0 b1 0
b
"0 "n-1 J
of order 2n is called the Hurwitz matrix of polynomials B and C, and its principal minors of even order are called the Hurwitz determinants of polynomials B and C.
The complex analog of the Hurwitz stability criterion: polynomial of degree n
f (z) = zn + zn-1 +... + ^ = 0 d0 d0 with complex coefficients and a unit (real and positive, but not necessarily unit) coefficient at the highest degree is stable
if and only if all Hurwitz determinants of polynomials g and h are positive.
In particular, when all coefficients of the resulting polynomial are real, we can use the classical Rouse-Hurwitz stability criterion for polynomials with real coefficients, or other similar [i5] criteria to check stability.
Thus, in general, to check the convergence of the Jacobi and Gauss-Seidel iterative methods, in order to avoid a direct search for the roots of a polynomial with complex or real coefficients, it is necessary to reduce it to a new one, which is checked for stability, which can be done using a computer by the above method. This method of checking convergence is especially relevant when the initial SLAEs have a large dimension, because of which we obtain equations of large powers, the solution of which is often very cumbersome.
Let's show that the range bounded by the boundary (34) is also obtained by applying the described method of checking convergence through the complex analog of the Hurwitz criterion for a polynomial with complex coefficients:
f (Я) = Я3 + pЯ + q k (z) = (z + i)3 + p( z + i)( z -1)2 + q( z -1)3 k(z) = (1 + p + q)z3 + (3 - p - 3q)z2 + +(3 - p + 3q) z + (1 + p - q) Assuming that 1 + p + q ^ 0, divide the last polynomial by this sum
k(z) = z3 +3 - p- 3q -2.3 - p+3q -.1+p- q
\ Im3 - p - 3q \ + p + q - Re3 - p + \ + p + q - Im \ + p - q \ + p + q 0 0
0 Re3 - p - 3q \ + p + q Im3 - p + 3q \ + p + q - Re \ + p -1 \ + p + q 0 0
0 \ Im3 - p - 3q \ + p + q - Re3 - p + 31 \ + p + q - Im \ + p - " \ + p + q 0
0 0 Re3 - p - 3q \ + p + q Im3 - p + 31 \ + p + q - Re \ + p - " \ + p + q 0
0 0 \ Im3 - p - 3q \ + p + q - Re3 - p + 31 \ + p + q - Im \ + p - " \ + p + q
0 0 0 Re3 - p - 3q \ + p + q Im3 - p + 31 \ + p + q - Re \ + p - " \ + p + q
\ + p + q
\ + p + q \ + p + q
When finding the convergence range in the real case, equating all imaginary elements in the obtained Hurwitz matrix to zero, we obtain the corresponding Hurwitz matrix, the principal minors of even order of which give the conditions we obtained above from the boundary (34).
In general, to check the convergence of the Jacobi method for a particular SLAE in three unknowns (in our case), we can program the described algorithm. For example, in the Python language:
import numpy as np
def complex_gurvic(p, q):
a = (3 - p - 3»q)/(l + p + q) b = (3 - p + 3*q)/(1 + p + q) c = (1 + p - q)/(l + p + q)
dl = np.array([[l, a.imag, -l#b.real, -lfc.imag, 0, 0], [0, a.real, b.imag, -l#c.real, 0, 0], [0, 1, a.imag, -l*b.real, c imag, 0], [0, 0, a.real, b.imag, -l*c.real, 0], [0, 0, 1, a.imag, -l*b.real, -l*c.imag], [0, 0, 0, a.real, b.imag, -l*c.real]]) deterl = np.linalg.det(dl[0:2, 0:2]) deter2 = np.linalg.det(dl[0:4, 0:4]) deter3 = np.linalg.det(dl) if deterl > 0 and deter2 > 0 and deter3 > 0:
print(1 The polynomial is stable. The method converges.1) else:
print (1 The polynomial is unstable. The method does not ■ . converge.')
~ . 3 3 - p - 3q 2 3 - p + 3q . 1 + p - q
k(iw) = -iw----- w + ---Liw + —-—-
1 + p + q 1 + p + q 1 + p + q
Let's separate the real and imaginary parts
7~/ - x / t^ 3 - p - 3q 2 t 3 - p + 3q 1 + p - q k(iw) = (-Re----w - Im----w + Re—-—-) +
\ + p + q
\ + p + <
\ + p + q
■t 3 t 3 - p -3q 2 n +i(-w -Im-—---w + Re
g (w) = - Re
i + p + q
3 - p - 3q 2 t ---—w - Im
i + p + q
3 - p - 3q, 2
3 - p + 3q i + p - q
w + Im-—-—-)
! + p + q 3 - p + 3q
h(w) = -w3 -1^^—— w2 + R^—r ' w + Im
i + p + q 3 - p + 3q
w + Re
\ + p + q \ + p - q
i + p + q
\+p - q
\+p+i
\+p+q
\+p+i
The degree of the polynomial k( z ) is odd, so
,, . 3 T 3 - p - 3q 2 „ 3 - p + 3q \ + p -t g = -h(w) = w +1^^-- w - R^^-- w - Im-—-—=
1 + p + q
1 + p + c
\ + p + i
~ „ 3 - p - 3q ^T3 - p + 3q 1 + p - q h = -g(w) = R^-—--- w +1^-—---w - Re-
\ + p + q
\ + p + q \ + p + q
The Hurwitz matrix for polynomials g, h has the form:
In general, the following conclusion can be made about the comparison of convergence of the two methods: in the equation (3) for the Gauss-Seidel method, it is always possible to take Я from the last line beyond the sign of the determinant, thus lowering the degree of the polynomial whose stability is to be investigated by one, which is not always possible for the Jacobi method according to the equation (2). Thus, for SLAEs in n > 2 unknowns with complex matrices, in general case, the polynomial, whose stability should be investigated, obtained for the Jacobi method, has degree by one more in contrast to the analogous polynomial for the Gauss-Seidel method.
VI. Statistical comparison of convergence of Jacobi
AND GAUSS-SEIDEL METHODS 100000 random matrices of SLAEs (1) with real matrix elements that are uniformly distributed random variables on the interval [-100; 100], with the number of unknowns from two to five, for each of them the well-known convergence criteria of each method were checked, then for each number of unknowns the number of cases in which both methods converge, only the Gauss-Seidel method converges, only the Jacobi method converges was determined. The obtained data are summarized in the table 1.
Table 1. Convergence results of Jacobi and Gauss-Seidel methods
Number of unknowns Both methods converge The Gauss-Seidel method converges, but the Jacobi method does not converge The Jacobi method converges, but the Gauss-Seidel method does not converge
2 49916 0 0
3 11818 7521 1095
4 1436 3411 528
5 111 726 76
The data obtained in the table for the number of unknowns n > 2 confirm the conclusions that, in general, the Gauss-Seidel method converges much more often than the Jacobi method, but the convergence of one of the methods cannot guarantee the convergence of the other. At the same time, we also see that as the number of unknowns in the SLAEs increases, both methods converge much less frequently, which is consistent with the above complex analog of the Hurwitz criterion.
Note also that in the case of SLAEs in two unknowns, the data from the table 1 confirm the conclusions that in this case both methods converge in the same way - if one converges, the other converges as well.
VII. Conclusion
The found boundary conditions in the complex case, as well as convergence ranges in the real case allowed us to see the picture of convergence conditions of Jacobi and Gauss-Seidel iterative methods and on this basis to make a comparative analysis of the effectiveness of each method: if in the case of square matrices of SLAEs in two unknowns both methods converge equally effectively, in the case of matrices of SLAEs in three and more unknowns methods have a noticeable difference in the convergence conditions -with increasing number of unknowns in SLAEs, the Gauss-Seidel method is noticeably more effective.
For example, in the case of an SLAEs' matrices in three unknowns, when the convergence ranges for both methods are plotted for the real case on the same coordinate plane, it can be seen that in the general case the Gauss-Seidel method has better convergence than the Jacobi method, since its convergence range is bounded by straight lines, but infinite in contrast to the convergence range of the Jacobi method, one of whose boundaries even enters the boundary of the convergence range of the Gauss-Seidel method. However, as it has been shown, the convergence range of the Gauss-Seidel method depends on the parameters that do not always give a full convergence range of the Jacobi method into the convergence range of the Gauss-Seidel method, because of which there may be situations when iterations converge to the exact solution by the Jacobi method but do not converge by the Gauss-Seidel method. Statistical comparison of convergence of both methods also confirms these conclusions.
When the number of unknowns over the field of complex numbers is large, the convergence of each method can be checked using the complex analog of the Hurwitz stability criterion, or using the classical Rouse-Hurwitz criterion in the real case.
References
[1] Bylina, J., & Bylina, B. (2008, October). Merging Jacobi and Gauss-Seidel methods for solving Markov chains on computer clusters. In 2008 International Multiconference on Computer Science and Information Technology (pp. 263-268). IEEE. DOI: 10.1109/IMCSIT.2008.4747250
[2] Nutzi, G., Schweizer, A., Moller, M., & Glocker, C. (2014, August). Projective jacobi and gauss-seidel on the gpu for non-smooth multi-body systems. In International Design Engineering Technical Conferences and Computers and Information in Engineering Conference (Vol. 46391, p. V006T10A013). American Society of Mechanical Engineers. DOI: 10.1115/DETC2014-34606
[3] Saad, Y., & Schultz, M. H. (1986). GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM Journal on scientific and statistical computing, 7(3), 856-869. DOI: 10.1137/0907058
[4] Tarigan, A. J. M., Mardiningsih, M., & Suwilo, S. (2022). The search for alternative algorithms of the iteration method on a system of linear equation. Sinkron: jurnal dan penelitian teknik informatika, 7(4), 2124-2424. DOI: 10.33395/sinkron.v7i4.11817
[5] Gunawardena, A. D., Jain, S. K., & Snyder, L. (1991). Modified iterative methods for consistent linear systems. Linear Algebra and Its Applications, 154, 123-143. DOI: 10.1016/0024-3795(91)90376-8
[6] Bagnara, R. (1995). A unified proof for the convergence of Jacobi and Gauss-Seidel methods. SIAM review, 37(1), 93-97. DOI: 10.1137/1037008
[7] Salkuyeh, D. K. (2007). Generalized Jacobi and Gauss-Seidel methods for solving linear system of equations. NUMERICAL MATHEMATICS-ENGLISH SERIES-, 16(2), 164.
[8] Chen, W. Y. (1995). On the polynomials with all their zeros on the unit circle. Journal of mathematical analysis and applications, 190(3), 714-724. DOI: 10.1006/jmaa.1995.1105
[9] Bharanedhar, S. V., Selvan, A. A., & Ghosh, R. (2023). Zeros of self-inversive polynomials with an application to sampling theory. Applied Mathematics and Computation, 439, 127547. DOI: 10.1016/j.amc.2022.127547
[10] Milaszewicz, J. P. (1987). Improving jacobi and gauss-seidel iterations. Linear Algebra and Its Applications, 93, 161-170. DOI: 10.1016/S0024-3795(87)90321-1
[11] Zadorozhniy, V. G. (2018). The conditions under which the roots of a polynomial lie inside the unit circle. Bulletin of VSU. Series: System Analysis and Information Technologies, 2, 22-25. https://www.elibrary.ru/item.asp?id=35449768
[12] Sun, L. Y. (2005). A comparison theorem for the SOR iterative method. Journal of computational and applied mathematics, 181(2), 336-341. DOI: 10.1016/j.cam.2004.12.007
[13] Ahmadi, A., Manganiello, F., Khademi, A., & Smith, M. C. (2021). A parallel Jacobiembedded Gauss-Seidel method. IEEE Transactions on Parallel and Distributed Systems, 32(6), 1452-1464. DOI: 10.1109/TPDS.2021.3052091
[14] Postnikov, M. M. (1981). Stable polynomials. Nauka", Moscow.
[15] Gantmakher, F. R. (2000). The theory of matrices (Vol. 131). American Mathematical Soc..
[16] Erd'elyi, T. (2001). On the zeros of polynomials with Littlewood-type coefficient constraints. Michigan Mathematical Journal, 49(1), 97-111. DOI: 10.1307/mmj/1008719037
[17] Konvalina, J., & Matache, V. (2004). Palindrome-polynomials with roots on the unit circle. Comptes Rendus Mathematiques, 26(2), 39.
[18] Mercer, I. D. (2006). Unimodular roots of special Littlewood polynomials. Canadian Mathematical Bulletin, 49(3), 438-447. DOI: 10.4153/CMB-2006-043-x
[19] Kohno, T., Kotakemori, H., Niki, H., & Usui, M. (1997). Improving the modified GaussSeidel method for Z-matrices. Linear Algebra and its Applications, 267, 113-123. DOI: 10.1016/S0024-3795(97)00063-3
[20] Li, W., & Sun, W. (2000). Modified Gauss-Seidel type methods and Jacobi type methods for Z-matrices. Linear Algebra and its Applications, 317(1-3), 227-240. DOI: 10.1016/S0024-3795(00)00140-3
[21] Shang, Y. (2009). A distributed memory parallel Gauss-Seidel algorithm for linear algebraic systems. Computers & Mathematics with Applications, 57(8), 1369-1376. DOI: 10.1016/j.camwa.2009.01.034
[22] Courtecuisse, H., & Allard, J. (2009, June). Parallel dense gauss-seidel algorithm on manycore processors. In 2009 11th IEEE International Conference on High Performance Computing and Communications (pp. 139-147). IEEE. DOI: 10.1109/HPCC.2009.51
[23] Koester, D. P., Ranka, S., & Fox, G. C. (1994, November). A parallel Gauss-Seidel algorithm for sparse power system matrices. In Supercomputing'94: Proceedings of the 1994 ACM/IEEE Conference
on Supercomputing (pp. 184-193). IEEE. DOI: 10.1145/602770.602806
[24] Amodio, P., & Mazzia, F. (1995). A parallel Gauss-Seidel method for block tridiagonal linear systems. SIAM Journal on Scientific Computing, 16(6), 1451-1461. DOI: 10.1137/0916084
[25] Tavakoli, R., & Davami, P. (2007). A new parallel Gauss-Seidel method based on alternating group explicit method and domain decomposition method. Applied mathematics and computation, 188(1), 713-719. DOI: 10.1016/j.amc.2006.10.023
[26] Karunanithi, S., Gajalakshmi, N., Malarvizhi, M., & Saileshwari, M. (2018). A Study on comparison of Jacobi, Gauss-Seidel and SOR methods for the solution in system of linear equations. Int. J. of Math. Trends and Technology, (I JMTT), 56(4). DOI: 10.14445/22315373/IJMTTV56P531
[27] Korsakov, G. F. (1973). The number of roots of a polynomial outside a circle. Mathematical notes of the Academy of Sciences of the USSR, 13, 3-8. DOI: 10.1007/BF01093620
[28] Biberdorf, E. A. D. (2000). A Criterion for the Dichotomy of Roots of a Polynomial on the Unit Circle. Sibirskii Zhurnal Industrial'noi Matematiki, 3(1), 16-32. https://www.elibrary.ru/item.asp?id=9484660
[29] Joyal, A., Labelle, G., & Rahman, Q. (1967). On the location of zeros of polynomials. Canadian mathematical bulletin, 10(1), 53-63. DOI: 10.4153/CMB-1967-006-3
[30] Dehmer, M. (2006). On the location of zeros of complex polynomials. Journal of Inequalities in Pure and Applied Mathematics, 7(1).
[31] Frank, E. (1946). On the zeros of polynomials with complex coefficients. DOI: 10.1090/S0002-9904-1946-08526-2
Сравнительный анализ итерационных методов Якоби и Гаусса-Зейделя
П. В. Храпов, Н. С. Волков
Аннотация — в работе дан сравнительный анализ итерационных численных методов Якоби и Гаусса-Зейделя решения систем линейных алгебраических уравнений (СЛАУ) с комплексными и действительными матрицами. Получены области сходимости для обоих методов для СЛАУ с двумя и тремя неизвестными, а также взаимосвязи данных областей. Построен алгоритм определения сходимости методов для СЛАУ с помощью комплексного аналога критерия Гурвица, дана реализация этого алгоритма на языке Python в случае СЛАУ с тремя неизвестными. Проведено статистическое сравнение сходимости обоих методов для СЛАУ с вещественной матрицей и количеством неизвестных от двух до пяти.
Ключевые слова — итерационные методы, система линейных алгебраических уравнений, метод Якоби, метод Гаусса-Зейделя, устойчивые многочлены, критерий Гурвица.
Об авторах:
Храпов Павел Васильевич,
доцент кафедры «Высшая математика», факультет «Фундаментальные науки», ФГБОУ ВО «Московский государственный технический университет имени Н.Э. Баумана (национальный исследовательский университет)» (105005, Российская Федерация, г. Москва, ул. 2-я Бауманская, д. 5, к. 1), кандидат физико-математических наук, ОЯСГО: https://orcid.org/0000-0002-6269-0727, khrapov@bmstu.ru , pvkhrapov@gmail.com .
Волков Никита Сергеевич,
студент кафедры «Высшая математика», факультет «Фундаментальные науки», ФГБОУ ВО «Московский государственный технический университет имени Н.Э. Баумана (национальный исследовательский университет)» (105005, Российская Федерация, г. Москва, ул. 2-я Бауманская, д. 5, к. 1),
ОЯСГО: https://orcid.org/0009-0008-7095-4092, nikita.volkov01@mail.ru , volkovns@student.bmstu.ru .
