# Comparative analysis of Jacobi and Gauss-Seidel iterative methods

Pavel Khrapov  Pavel Khrapov
Department of Mathematics
Bauman Moscow State Technical University
ul. Baumanskaya 2-ya, 5/1, Moscow
105005, Moscow, Russian Federation<EMAIL_ADDRESS><EMAIL_ADDRESS>and  Nikita Volkov  Nikita Volkov
Department of Mathematics
Bauman Moscow State Technical University
ul. Baumanskaya 2-ya, 5/1, Moscow
105005, Moscow, Russian Federation<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### 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.

###### Key words and phrases:

iterative methods, system of linear algebraic equations, Jacobi method, Gauss-
Seidel method, stable polynomials, Hurwitz criterion.

###### 2010 Mathematics Subject Classification:

65F10

##  1\. 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 being 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 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 predominance 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 ranges of convergence 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
matrix 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.

##  2\. Convergence conditions of the Jacobi and Gauss-Seidel methods

When solving a system of linear algebraic equations

$A\vec{x}=\vec{b}$ (2.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\neq 0$

where $L,D,R,$ are, respectively, the matrices with subdiagonal, diagonal, and
overdiagonal elements of matrix $A$, and then there is a system is obtained
from the original SLAE (2.1):

$\vec{x}=-D^{-1}(L+R)\vec{x}+D^{-1}\vec{b}$

for which the simple iterative method converges if all roots of the equation

$\begin{vmatrix}\lambda a_{11}&a_{12}&...&a_{1n}\\\ a_{21}&\lambda
a_{22}&...&a_{2n}\\\ ...&...&...&...\\\ a_{n1}&a_{n2}&...&\lambda
a_{nn}\end{vmatrix}=0$ (2.2)

have absolute values less than one [6]; $a_{ij}$ \- elements of the original
matrix $A$, $a_{ij}\in\mathbb{C}$.
Similarly, the Gauss-Seidel method transforms the original SLAE (2.1) to a
system:

$\vec{x}=-(L+D)^{-1}R\vec{x}+(L+D)^{-1}\vec{b}$

for which the simple iterative method converges to an exact solution if all
roots of equation

$\begin{vmatrix}\lambda a_{11}&a_{12}&...&a_{1n}\\\ \lambda a_{21}&\lambda
a_{22}&...&a_{2n}\\\ ...&...&...&...\\\ \lambda a_{n1}&\lambda
a_{n2}&...&\lambda a_{nn}\end{vmatrix}=0$ (2.3)

have absolute values less than one [6], $a_{ij}\in\mathbb{C}$.
When the dimension of the original SLAE is small, we can find the convergence
ranges of the methods by directly solving the equations (2.2) and (2.3). Let
us show this for the cases of SLAEs in two and three unknowns, which often
arise in applied research.

##  3\. System of linear algebraic equations in two unknowns

###  3.1. Jacobi method

The equation (2.2) has the form:

$\begin{vmatrix}\lambda a_{11}&a_{12}\\\ a_{21}&\lambda
a_{22}\end{vmatrix}=\lambda^{2}a_{11}a_{22}-a_{12}a_{21}=0,$ (3.1)

and for the convergence of the method it is necessary that its roots lie
inside the unit circle.
In the general case $a_{ij},\lambda_{1,2}\in\mathbb{C}$, and the system (2.2)
is equivalent to equations

$(\lambda-r_{1}e^{i\varphi_{1}})(\lambda-r_{2}e^{i\varphi_{2}})=0$
$\lambda^{2}+\lambda(-r_{1}e^{i\varphi_{1}}-r_{2}e^{i\varphi_{2}})+r_{1}r_{2}e^{i\varphi_{1}}e^{i\varphi_{2}}=0,$
(3.2)

where $r_{1}e^{i\varphi_{1}},r_{2}e^{i\varphi_{2}}$ are the roots of equation
(3.1), $r_{1},r_{2}<1$.
Comparing (3.1) and (3.2), we obtain the system ($a_{11}a_{22}\neq 0$, since
$a_{11},a_{22}$ are elements of the diagonal matrix $D$):

$\begin{cases}r_{1}e^{i\varphi_{1}}+r_{2}e^{i\varphi_{2}}=0\\\
r_{1}r_{2}e^{i\varphi_{1}}e^{i\varphi_{2}}=-\frac{a_{12}a_{21}}{a_{11}a_{22}}\end{cases}$
(3.3)

from which follow:

$r_{1}=r_{2}$ $r_{1}r_{2}=r_{1}^{2}=|\frac{a_{12}a_{21}}{a_{11}a_{22}}|<1$
$|a_{12}a_{21}|<|a_{11}a_{22}|$ (3.4)

The condition (3.4) 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 (2.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.

###  3.2. Gauss-Seidel method

The equation (2.3) has the form:

$\begin{vmatrix}\lambda a_{11}&a_{12}\\\ \lambda a_{21}&\lambda
a_{22}\end{vmatrix}=\lambda^{2}a_{11}a_{22}-\lambda a_{12}a_{21}=0$ (3.5)

and its roots

$\lambda_{1}=0$ $\lambda_{2}=\frac{a_{12}a_{21}}{a_{11}a_{22}}.$

must have an absolute value less than one ($a_{11}a_{22}\neq 0$, since
$a_{11},a_{22}$ are diagonal elements of the 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

$|a_{12}a_{21}|<|a_{11}a_{22}|$ (3.6)

Thus, both the Jacobi method and the Gauss-Seidel method for SLAEs in two
unknowns have the same range of convergence (3.6).

##  4\. System of linear algebraic equations with three unknowns

###  4.1. Jacobi method

The equation (2.2) has the form:

$\lambda^{3}a_{11}a_{22}a_{33}-\lambda(a_{13}a_{22}a_{31}+a_{23}a_{11}a_{32}+a_{12}a_{33}a_{21})+(a_{13}a_{32}a_{21}+a_{12}a_{23}a_{31})=0$
(4.1)

and for convergence of the method it is necessary that all its roots lie
inside the unit circle.
Let’s divide it by $a_{11}a_{22}a_{33}$ (there are no zero elements on the
diagonal of matrix $A$, since matrix $D$ must have an inverse):

$\lambda^{3}+\frac{-a_{13}a_{22}a_{31}-a_{23}a_{11}a_{32}-a_{12}a_{33}a_{21}}{a_{11}a_{22}a_{33}}\lambda+\frac{a_{13}a_{32}a_{21}+a_{12}a_{23}a_{31}}{a_{11}a_{22}a_{33}}=0.$

Denoting

$p=\frac{-a_{13}a_{22}a_{31}-a_{23}a_{11}a_{32}-a_{12}a_{33}a_{21}}{a_{11}a_{22}a_{33}}$
$q=\frac{a_{13}a_{32}a_{21}+a_{12}a_{23}a_{31}}{a_{11}a_{22}a_{33}},$

we obtain the canonical cubic equation:

$\lambda^{3}+p\lambda+q=0,$ (4.2)

In the general case its coefficients and roots are complex:
$p,q,\lambda_{1,2,3}\in\mathbb{C}$.
We find the range of convergence of the method expressed in terms of
$p,q\in\mathbb{C}$ 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 (4.2) 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 (4.2) have a unit absolute value, and the other two
roots have an absolute value not exceeding one:

$(\lambda-e^{i\varphi_{1}})(\lambda-r_{2}e^{i\varphi_{2}})(\lambda-
r_{3}e^{i\varphi_{3}})=0$ (4.3)
$\lambda^{3}+\lambda^{2}(-e^{i\varphi_{1}}-r_{2}e^{i\varphi_{2}}-r_{3}e^{i\varphi_{3}})+\lambda(r_{2}e^{i\varphi_{1}}e^{i\varphi_{2}}+r_{3}e^{i\varphi_{1}}e^{i\varphi_{3}}+r_{2}r_{3}e^{i\varphi_{2}}e^{i\varphi_{3}})-r_{2}r_{3}e^{i\varphi_{1}}e^{i\varphi_{2}}e^{i\varphi_{3}}=0$
$r_{2},r_{3}\leq 1$

Comparing (4.2) and (4.3), we obtain a system of equations for the first
boundary of the convergence range:

$\begin{cases}e^{i\varphi_{1}}+r_{2}e^{i\varphi_{2}}+r_{3}e^{i\varphi_{3}}=0\\\
r_{2}e^{i\varphi_{1}}e^{i\varphi_{2}}+r_{3}e^{i\varphi_{1}}e^{i\varphi_{3}}+r_{2}r_{3}e^{i\varphi_{2}}e^{i\varphi_{3}}=p\\\
-r_{2}r_{3}e^{i\varphi_{1}}e^{i\varphi_{2}}e^{i\varphi_{3}}=q\end{cases}$
(4.4)

from which follow:

$|q|\leq 1$ (4.5) $\arg(q)=\pi+\varphi_{1}+\varphi_{2}+\varphi_{3}$ (4.6)
$r_{3}e^{i\varphi_{3}}=-r_{2}e^{i\varphi_{2}}-e^{i\varphi_{1}}$ (4.7)

The first equation of the system (4.4) has a geometric interpretation (fig.
1):

Figure 1. Geometric interpretation of the first equation of the system (4.4)

Substituting (4.7) into the second and third equations of the system (4.4), we
have expressions:

$p=-e^{i\varphi_{1}}e^{i\varphi_{1}}-r_{2}e^{i\varphi_{1}}e^{i\varphi_{2}}-r_{2}^{2}e^{i\varphi_{2}}e^{i\varphi_{2}}$
$q=r_{2}e^{i\varphi_{1}}e^{i\varphi_{1}}e^{i\varphi_{2}}+r_{2}^{2}e^{i\varphi_{1}}e^{i\varphi_{2}}e^{i\varphi_{2}}$

comparing which, we obtain the equation of the first boundary of the
convergence range of the Jacobi method in the general case:

$p=-qe^{-i\varphi_{1}}-e^{2i\varphi_{1}}$ (4.8)

from which follow the relationships of absolute value and argument for $p$ and
$q$:

$|p|=\sqrt{r_{q}^{2}+1+2r_{q}\cos(\varphi_{q}-3\varphi_{1})}$ (4.9)
$\arg(p)=\arctan(\frac{r_{q}\sin(\varphi_{q}-\varphi_{1})+\sin(2\varphi_{1})}{r_{q}\cos(\varphi_{q}-\varphi_{1})+\cos(2\varphi_{1})})$
(4.10)

The relationships (4.9) and (4.10) show that the absolute value and argument
for $p$ depend on three parameters - the absolute value $r_{q}$, the argument
$\varphi_{1}$, and the argument $\varphi_{q}$, the last of which depends not
only on $\varphi_{1}$ (4.6), so we will take the argument $\varphi_{1}$ as a
parameter to visualize the absolute value and argument for $p$.
Let us take for example $\varphi_{1}=0$ and $\varphi_{1}=\pi$ (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 (4.2) is real): at
$\varphi_{1}=0$, according to geometrical considerations (fig. 1) and (4.6)

$\varphi_{q}=\varphi_{2}+\varphi_{3}+\pi\in[-\pi,-\frac{\pi}{2}]\cup[\frac{\pi}{2},\pi]$
(4.11) $p=-q-1$ (4.12)

Similarly, when $\varphi_{1}=\pi$

$\varphi_{q}=\varphi_{2}+\varphi_{3}+\pi+\pi\in[-\frac{\pi}{2},\frac{\pi}{2}]$
(4.13) $p=q-1$ (4.14)

The picture of the absolute value (4.9) in this case is as follows (fig. 2):

Figure 2. Dependence of the absolute value $|p|$ at the boundary (4.8) at
$\varphi_{1}=0$ ($\operatorname{Re}(q)\leq 0$) and $\varphi_{1}=\pi$
($\operatorname{Re}(q)\geq 0$)

To find the second boundary of the convergence range, consider the case when
two roots of the equation (4.2) 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
$\varphi_{1}$, which is the argument of the first root: taking into account
that before the rotation by the angle $\varphi_{1}$ one root was real, and the
other two roots were complex-conjugate with arguments $\varphi_{2}$ and
$-\varphi_{2}$ respectively, after the rotation the picture on the complex
plane will be as follows (fig. 4, fig. 4):

Figure 3. Location of roots of the equation (4.2) on the complex plane before
rotation

Figure 4. Location of roots of the equation (4.2) on the complex plane after
rotation

$(\lambda-
r_{1}e^{i\varphi_{1}})(\lambda-e^{i(\varphi_{2}+\varphi_{1})})(\lambda-e^{i(-\varphi_{2}+\varphi_{1})})=0$
(4.15) $r_{1}\leq 1$

Opening the brackets and comparing (4.15) with equation (4.2), we obtain the
system for the second boundary of the convergence range:

$\begin{cases}e^{i(\varphi_{2}+\varphi_{1})}+e^{i(-\varphi_{2}+\varphi_{1})}+r_{1}e^{i\varphi_{1}}=0\\\
e^{2i\varphi_{1}}+r_{1}e^{(\varphi_{1}+\varphi_{2})}e^{i\varphi_{1}}+r_{1}e^{(\varphi_{1}-\varphi_{2})}e^{i\varphi_{1}}=p\\\
-r_{1}e^{3i\varphi_{1}}=q\end{cases}$ (4.16)

from which follow:

$|q|\leq 1$ (4.17)
$r_{1}e^{i\varphi_{1}}=-e^{i(\varphi_{2}+\varphi_{1})}-e^{i(-\varphi_{2}+\varphi_{1})}$
(4.18)

Let’s substitute (4.18) into the second and third equations of the system
(4.16):

$p=-e^{2i(\varphi_{1}+\varphi_{2})}-e^{2i\varphi_{1}}-e^{2i(\varphi_{1}-\varphi_{2})}$
(4.19)
$q=e^{2i\varphi_{1}}(e^{i(\varphi_{1}+\varphi_{2})}+e^{i(\varphi_{1}-\varphi_{2})})$
$qe^{-2i\varphi_{1}}=e^{i(\varphi_{1}+\varphi_{2})}+e^{i(\varphi_{1}-\varphi_{2})}$
(4.20)

Comparing the square of the expression (4.20) and the expression (4.19), we
obtain the equation of the second boundary of the convergence range of the
Jacobi method in general case:

$p=-q^{2}e^{-4i\varphi_{1}}+e^{2i\varphi_{1}}$ (4.21)
$\arg(q)=3\varphi_{1}+\pi$ (4.22)

from which we find the relationships of absolute value and argument for $p$
and $q$ subject to the condition (4.22):

$|p|=1-r_{1}^{2}=1-|q|^{2}$ (4.23)
$\arg(p)=2\varphi_{1}=\frac{2}{3}\arg(q)-\frac{2}{3}\pi$ (4.24)

When $\varphi_{1}\in[-\pi,\pi]$, the picture for the absolute value (4.23) is
as follows (fig. 5):

Figure 5. Dependence of the absolute value $|p|$ on $q$

It follows from (4.22), (4.24) that when $\varphi_{1}=\pm\pi$, arguments

$\arg(q)=\arg(p)=0$

are in the same phase corresponding to the real case, and when $\varphi_{1}=0$

$\arg(q)=\pi$ $\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 (4.8), it follows from (4.10) that for
$\varphi_{q}\in\\{0,\pm\pi\\}$ and $\varphi_{1}\in\\{0,\pm\pi\\}$ , as in the
case of the second bound:

$\arg(p)\in\\{0,\pm\pi\\}$

Hence, when at least one of the roots of the equation (4.2) is real, the cross
section of the convergence range boundaries shown in the figures (fig. 2, fig.
5), at $\operatorname{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 (4.12)
and (4.14), moves to the area of negative values of $p$):

Figure 6. Boundary of the convergence range of the Jacobi method for the real
case

In particular, when $\varphi_{1}\in\\{0,\pm\pi\\}$, the equation (4.21) takes
the form of a parabola:

$p=q^{2}+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 (4.8) and (4.21) and the condition:

$|q|\leq 1$

In the real case, when $a_{ij}\in A$; $p,q$ are real numbers, and the roots of
equation (4.2) are either all real or two of them are complex-conjugate,
equations (4.8) and (4.21) form the following boundary of the convergence
range (fig. 6):

$\left[\begin{gathered}p=-q-1\\\ p=q-1\\\ p=-q^{2}+1\\\ -1\leq q\leq
1\end{gathered}\right.$ (4.25)

The boundary, due to the above assumptions that it contains at least one root
of the equation (4.2) 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 (4.25), which does not belong to this region, so for
the convergence range, given the conditions (4.5) and (4.17), it follows that:

$-1<q<1$

Let’s show that the area in the figure (fig. 6) can be filled completely:

$p=-\frac{a_{13}a_{31}}{a_{11}a_{33}}-\frac{a_{23}a_{32}}{a_{22}a_{33}}-\frac{a_{12}a_{21}}{a_{11}a_{22}}=-\frac{a_{13}a_{32}a_{21}a_{22}a_{31}}{a_{11}a_{22}a_{33}a_{21}a_{32}}-\frac{a_{12}a_{23}a_{31}a_{11}a_{32}}{a_{11}a_{22}a_{33}a_{12}a_{31}}-\frac{a_{12}a_{21}}{a_{11}a_{22}}$

Let’s denote:

$x=\frac{a_{22}}{a_{21}}$ $y=\frac{a_{11}}{a_{12}}$ $t=\frac{a_{31}}{a_{32}}$
$p=-\frac{a_{13}a_{32}a_{21}}{a_{11}a_{22}a_{33}}xt-\frac{a_{12}a_{23}a_{31}}{a_{11}a_{22}a_{33}}\frac{y}{t}-\frac{1}{xy}$

Assuming that:

$xt=\frac{y}{t}=a$ (4.26) $t\neq 0$ $xy=a^{2}$ $p=-aq-\frac{1}{a^{2}}$ (4.27)

$x,y,t$ are independent of each other, so it is always possible to choose the
coefficients of the matrix $A$ of the system (2.1) such that the condition
(4.26) is satisfied and the line (4.27) 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 (4.27), 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.

###  4.2. Gauss-Seidel method

The equation (2.3) has the form:

$\lambda^{3}a_{11}a_{22}a_{33}+\lambda^{2}(a_{21}a_{13}a_{32}-a_{13}a_{22}a_{31}-a_{32}a_{11}a_{23}-a_{21}a_{33}a_{12})+\lambda
a_{12}a_{23}a_{31}=0$ (4.28)

and the method converges if all its roots lie inside the unit circle.
One of the roots of the equation (4.28) is zero, and the other two roots are
found from the quadratic equation:

$\lambda^{2}a+\lambda d+b=0$ (4.29)

where

$d=a_{21}a_{13}a_{32}-a_{13}a_{22}a_{31}-a_{32}a_{11}a_{23}-a_{21}a_{33}a_{12}$
$a=a_{11}a_{22}a_{33}$ $b=a_{12}a_{23}a_{31}$

For the Gauss-Seidel method, we find the convergence range, given that within
it all roots of the equation (4.29) have an absolute value not exceeding one.
In the general case $a,b,d,\lambda_{1,2}\in\mathbb{C}$ and the equation (4.29)
is equivalent to equation:

$a\lambda^{2}-\lambda
a(r_{1}e^{i\varphi_{1}}+r_{2}e^{i\varphi_{2}})+ar_{1}r_{2}e^{i\varphi_{1}}e^{i\varphi_{2}}=0$
$0<r_{1},r_{2}<1$

comparing it with (4.29), obtain the system:

$\begin{cases}-ar_{1}e^{i\varphi_{1}}-ar_{2}e^{i\varphi_{2}}=d\\\
ar_{1}r_{2}e^{i\varphi_{1}}e^{i\varphi_{2}}=b\end{cases}$ (4.30)

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\neq 0$,
since $a_{ii}$ are the diagonal elements of the triangular matrix $L+D$,
$i\in\\{1,2,3\\}$):

$|\frac{b}{a}|<1$ (4.31)

and the boundary of the convergence range (4.30), assuming that at least one
of the roots of the equation (4.29) has a unit absolute value on it (let
$r_{1}=1$), is given by the conditions:

$\begin{cases}d=-ae^{i\varphi_{1}}-be^{-i\varphi_{1}}\\\
ar_{2}e^{i\varphi_{1}}e^{i\varphi_{2}}=b\\\ 0\leq r_{2}\leq 1\end{cases}$
(4.32)

For $d_{1}=\frac{d}{a}$ and $b_{1}=\frac{b}{a}$ on the boundary (4.32), we can
also find the relationships between absolute values and arguments:

$d_{1}=-e^{i\varphi_{1}}-b_{1}e^{-i\varphi_{1}}\\\ $
$d_{1}=-\cos\varphi_{1}-r_{b_{1}}\cos(\varphi_{b_{1}}-\varphi_{1})-i(\sin\varphi_{1}+r_{b_{1}}\sin(\varphi_{b_{1}}-\varphi_{1}))$
$\varphi_{b_{1}}=\varphi_{1}+\varphi_{2}$
$|d_{1}|=\sqrt{1+r_{b_{1}}^{2}+2r_{b_{1}}\cos(\varphi_{1}-\varphi_{2})}$
(4.33)
$\arg(d_{1})=\arctan(\frac{\sin\varphi_{1}+r_{b_{1}}\sin\varphi_{2}}{\cos\varphi_{1}+r_{b_{1}}\cos\varphi_{2}})$
(4.34)

In particular, when the matrix $A$ of the system (2.1) contains real
coefficients, $a,b,d\in\mathbb{R}$, solving directly the quadratic equation
(4.29) and applying the convergence criterion of the Gauss-Seidel method,
taking into account the condition (4.31), we obtain the convergence range of
the Gauss-Seidel method in the case of real roots of the equation (4.29):

$\begin{cases}|d|<|a+b|\\\ |\frac{b}{a}|<1\end{cases}$ (4.35)

and in the case of complex-conjugate roots of the equation (4.29).:

$0<\frac{b}{a}<1$ (4.36)

Note that in the latter case, the condition

$|d|<|a+b|$

follows directly from the condition (4.36) and the negativity of the
discriminant of the equation (4.29): $d^{2}<4ab$. Therefore, the system (4.35)
is a single range of convergence of the Gauss-Seidel method in the case of
real matrix elements of the system (2.1).
The first condition of the system (4.35) is interpreted as a segment $d$ on an
infinite line.
Note that the conditions (4.35) are consistent with the boundary (4.32).
Unlike the Jacobi method, the convergence range of the Gauss-Seidel method in
the case of real coefficients of the equation (4.29) 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 (2.1). For this purpose, we construct the
convergence range bounded by the boundary (4.25) and the range (4.35) 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 (4.35), we obtain:

$\begin{cases}|(p+q)a-b|<|a+b|\\\ |\frac{b}{a}|<1\end{cases}$ (4.37)

Expanding the absolute values in the first inequality of the system (4.37), we
find that one of the boundaries of the convergence range of the Gauss-Seidel
method is always the line

$p=-q-1$ (4.38)

which is also one of the boundaries (4.25) of the convergence range of the
Jacobi method, and the second one is also a straight line, which has the
following form:

$p=-q+\frac{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 (4.37) follows:

$-1<\frac{a+2b}{a}<3,$

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 (4.38),
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):

Figure 7. Convergence ranges of methods on the plane $qOp$. The band $\alpha$
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 figure 7, the advantages of the Gauss-Seidel method over the
Jacobi method when the system (2.1) has real matrix elements are obvious (in
the case in Figure 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:

$\begin{cases}|2(p+q)-1|<3\\\ |\frac{b}{a}|=\frac{1}{2}<1\end{cases}$

thus

$\frac{a+2b}{a}=2$

Then, by analogy with fig. 7, the convergence ranges for each method on the
plane $qOp$ look as follows (fig. 8):

Figure 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:

$A=\begin{pmatrix}-8&6&-4\\\ -9&8&6\\\ 4&-5&3\end{pmatrix}$

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:

$\begin{cases}|-192(p+q)-144|<48\\\
|\frac{b}{a}|=|\frac{144}{192}|<1\end{cases}$ $\frac{a+2b}{a}=-\frac{1}{2}$

By analogy with fig. 7, we obtain the following picture of convergence ranges
for both methods on the plane $qOp$ (fig. 9):

Figure 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:

$p=-\frac{50}{192}$ $q=\frac{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$.

##  5\. 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.2) or (2.3):

$f(\lambda)=a_{0}\lambda^{n}+a_{1}\lambda^{n-1}+...+a_{n}=0,a_{0}\neq 0$ (5.1)

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)=a_{0}(z+1)^{n}+a_{1}(z+1)^{n-1}(z-1)+a_{2}(z+1)^{n-2}(z-1)^{2}+...+a_{n}(z-1)^{n}=0$

obtained from (5.1) 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)=d_{0}z^{n}+d_{1}z^{n-1}+...+d_{n},$

it’s equivalent, under the assumption that $d_{0}\neq 0$, to the polynomial
whose first coefficient is equal to one:

$\widetilde{f}(z)=z^{n}+\frac{d_{1}}{d_{0}}z^{n-1}+...+\frac{d_{n}}{d_{0}}$
(5.2)

Replacing in (5.2) the variable $z$ by a purely imaginary number
$i\omega,\omega\in\mathbb{R}$, we have the polynomial:

$\widetilde{f}(i\omega)=(i\omega)^{n}+\frac{d_{1}}{d_{0}}(i\omega)^{n-1}+...+\frac{d_{n}}{d_{0}}$

which, by raising the multiplier $i\omega$ 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.

$\widetilde{f}(i\omega)=g(\omega)+ih(\omega),$

for which, according to [14], if the degree of the polynomial (5.2) is $n=2m$:

$\widetilde{g}=(-1)^{m}g$ $\widetilde{h}=(-1)^{m-1}h,$

if $n=2m+1$:

$\widetilde{g}=(-1)^{m}h$ $\widetilde{h}=(-1)^{m}g.$

Let

$B=b_{0}x^{n}+b_{1}x^{n-1}+...+b_{n}$

be an arbitrary polynomial of degree $n$ with real coefficients with positive
prime factor $b_{0}$, and let

$C=c_{0}x^{n-1}+c_{1}x^{n-2}+...+c_{n-1}$

be an arbitrary polynomial of degree at most $n-1$ with real coefficients.
Definition. Square matrix

$\begin{pmatrix}b_{0}&b_{1}&b_{2}&b_{3}&...&0\\\ 0&c_{0}&c_{1}&c_{2}&...&0\\\
0&b_{0}&b_{1}&b_{2}&...&0\\\ 0&0&c_{0}&c_{1}&...&0\\\
...&...&...&...&...&...\\\ 0&...&b_{0}&b_{1}&...&b_{n}\\\
0&...&0&c_{0}&...&c_{n-1}\end{pmatrix}$

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$

$\widetilde{f}(z)=z^{n}+\frac{d_{1}}{d_{0}}z^{n-1}+...+\frac{d_{n}}{d_{0}}=0$

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 $\widetilde{g}$ and $\widetilde{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 [15] 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 (4.25) 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(\lambda)=\lambda^{3}+p\lambda+q$
$k(z)=(z+1)^{3}+p(z+1)(z-1)^{2}+q(z-1)^{3}$
$k(z)=(1+p+q)z^{3}+(3-p-3q)z^{2}+(3-p+3q)z+(1+p-q)$

Assuming that $1+p+q\neq 0$, divide the last polynomial by this sum

$\widetilde{k}(z)=z^{3}+\frac{3-p-3q}{1+p+q}z^{2}+\frac{3-p+3q}{1+p+q}z+\frac{1+p-q}{1+p+q}$
$\widetilde{k}(iw)=-iw^{3}-\frac{3-p-3q}{1+p+q}w^{2}+\frac{3-p+3q}{1+p+q}iw+\frac{1+p-q}{1+p+q}$

Let’s separate the real and imaginary parts

$\widetilde{k}(iw)=(-\operatorname{Re}\frac{3-p-3q}{1+p+q}w^{2}-\operatorname{Im}\frac{3-p+3q}{1+p+q}w+\operatorname{Re}\frac{1+p-q}{1+p+q})+$
$+i(-w^{3}-\operatorname{Im}\frac{3-p-3q}{1+p+q}w^{2}+\operatorname{Re}\frac{3-p+3q}{1+p+q}w+\operatorname{Im}\frac{1+p-q}{1+p+q})$
$g(w)=-\operatorname{Re}\frac{3-p-3q}{1+p+q}w^{2}-\operatorname{Im}\frac{3-p+3q}{1+p+q}w+\operatorname{Re}\frac{1+p-q}{1+p+q}$
$h(w)=-w^{3}-\operatorname{Im}\frac{3-p-3q}{1+p+q}w^{2}+\operatorname{Re}\frac{3-p+3q}{1+p+q}w+\operatorname{Im}\frac{1+p-q}{1+p+q}$

The degree of the polynomial $\widetilde{k}(z)$ is odd, so

$\tilde{g}=-h(w)=w^{3}+\operatorname{Im}\frac{3-p-3q}{1+p+q}w^{2}-\operatorname{Re}\frac{3-p+3q}{1+p+q}w-\operatorname{Im}\frac{1+p-q}{1+p+q}$
$\tilde{h}=-g(w)=\operatorname{Re}\frac{3-p-3q}{1+p+q}w^{2}+\operatorname{Im}\frac{3-p+3q}{1+p+q}w-\operatorname{Re}\frac{1+p-q}{1+p+q}$

The Hurwitz matrix for polynomials $\tilde{g},\tilde{h}$ has the form:

$\begin{pmatrix}1&\operatorname{Im}\frac{3-p-3q}{1+p+q}&-\operatorname{Re}\frac{3-p+3q}{1+p+q}&-\operatorname{Im}\frac{1+p-q}{1+p+q}&0&0\\\
0&\operatorname{Re}\frac{3-p-3q}{1+p+q}&\operatorname{Im}\frac{3-p+3q}{1+p+q}&-\operatorname{Re}\frac{1+p-q}{1+p+q}&0&0\\\
0&1&\operatorname{Im}\frac{3-p-3q}{1+p+q}&-\operatorname{Re}\frac{3-p+3q}{1+p+q}&-\operatorname{Im}\frac{1+p-q}{1+p+q}&0\\\
0&0&\operatorname{Re}\frac{3-p-3q}{1+p+q}&\operatorname{Im}\frac{3-p+3q}{1+p+q}&-\operatorname{Re}\frac{1+p-q}{1+p+q}&0\\\
0&0&1&\operatorname{Im}\frac{3-p-3q}{1+p+q}&-\operatorname{Re}\frac{3-p+3q}{1+p+q}&-\operatorname{Im}\frac{1+p-q}{1+p+q}\\\
0&0&0&\operatorname{Re}\frac{3-p-3q}{1+p+q}&\operatorname{Im}\frac{3-p+3q}{1+p+q}&-\operatorname{Re}\frac{1+p-q}{1+p+q}\end{pmatrix}$

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 (4.25).
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:

In general, the following conclusion can be made about the comparison of
convergence of the two methods: in the equation (2.3) for the Gauss-Seidel
method, it is always possible to take $\lambda$ 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.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.

##  6\. Statistical comparison of convergence of Jacobi and Gauss-Seidel
methods

100000 random matrices of SLAEs (2.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.

##  7\. 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-Gurwitz 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] Nützi, G., Schweizer, A., Möller, 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 Jacobi-embedded 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élyi, 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 Gauss-Seidel 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 many-core 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,(IJMTT), 56(4). DOI: 10.14445/22315373/IJMTT-V56P531
  * [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, È. 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