Thomas Smotzer, John Buoni
Department of Mathematics and Statistics, Youngstown State University, Youngstown, Ohio, USA.
DOI: 10.4236/ajcm.2024.142009   PDF    HTML   XML   128 Downloads   593 Views   Citations

Abstract

Kellogg gave a version of the Peaceman-Radford method. In this paper, we introduce a SSOR iteration method which uses Kellogg’s method. The new algorithm has some advantages over the traditional SSOR algorithm. A Cyclic Reduction algorithm is introduced via a decoupling in Kellogg’s method.

Keywords

Matrix Splitting, SSOR Iteration, KSSOR Iteration Method, Kellogg-Type SSOR Iteration, Cyclic Reduction

Share and Cite:

1. Introduction

Throughout this paper, one studies a solution method for the system of linear equations which occurs throughout the literature [1] - [11] .

$Ax=b$ (1)

where $A=\left(a_{ij}\right)\in R^{n\times n}$ is non-singular matrix with unit diagonal b and $x\in R^n$ with x being an unknown vector in $R^n$ .

Splitting the matrix A into

$A=I-L-U$ (2)

where I is the identity matrix, L and U are the strictly lower and strictly upper triangular parts of the matrix A, respectively; see [8] and [9] , the successive overrelaxation iteration formula for the SOR iteration method for an optimal relaxation factor $\omega$ is as follows:

$x^{k+1}={\left(I-\omega L\right)}^{-1}\left(\left(1-\omega \right)I+\omega U\right)x^k+\omega b\mathrm{.}$ (3)

While the following 2-step scheme:

$\begin{array}{l}\left(I-\omega L\right)y^k=\left(\left(1-\omega \right)I+\omega U\right)x^k+\omega b\mathrm{,}\\ \left(I-\omega U\right)x^{k+1}=\left(\left(I-\omega \right)I+\omega L\right)y^k+\omega b\end{array}$ (4)

for any initial guess $x^0$ defines the iteration formula for the SSOR iteration method with relaxation factor $\omega >0$ , [12] and [10] . The study of allowable $\omega$ for convergence of (4) appears throughout the literature [2] and [5] .

Define the iteration matrix for the SOR iteration method as

$\Delta \left(\omega \right)\mathrm{<mo>\; :\; </mo>}={\left(I-\omega L\right)}^{-1}\left(\left(1-\omega \right)I+\omega U\right)\mathrm{,}$ (5)

and the iteration matrix for the SSOR iteration method as

$\Gamma \left(\omega \right)\mathrm{<mo>\; :\; </mo>}={\left(I-\omega U\right)}^{-1}\left(\left(1-\omega \right)I+\omega L\right){\left(I-\omega L\right)}^{-1}\left(\left(1-\omega \right)I+\omega U\right)\mathrm{.}$ (6)

The required two step iteration is a drawback of the SSOR iteration method, when compared to SOR iteration method. However, for a symmetric matrix A, the iteration matrix has real spectrum, and the spectral radius of $\Gamma \left(\omega \right)$ i.e. $\rho \left(\Gamma \left(\omega \right)\right)$ , is equal to the square of a norm of $\Delta \left(\omega \right)$ ( [9] : pg. 237). This paper introduces a variation of the SSOR iteration method based on the work of Kellogg [13] and is similar to work that appears in [7] for the HSS algorithm. This new algorithm requires approximately the amount of work per iteration as the SOR iteration method, and uses an iteration matrix which is similar to the SSOR iteration method matrix with equal spectral radius. This iteration scheme will be entitled the Kellogg-type SSOR iteration method, and its iteration matrix will be denoted by the KSSOR iteration matrix.

2. Different Iterations

Letting the terms in $\Gamma \left(\omega \right)$ be somewhat permuted, one defines

$\begin{array}{l}\Psi \left(\omega \right):={\left(I-\omega L\right)}^{-1}\left(\left(1-\omega \right)I+\omega L\right){\left(I-\omega U\right)}^{-1}\left(\left(1-\omega \right)I+\omega U\right),\\ \Theta \left(\omega \right):={\left(I-\omega U\right)}^{-1}\left(\left(1-\omega \right)I+\omega U\right){\left(I-\omega L\right)}^{-1}\left(\left(1-\omega \right)I+\omega L\right).\end{array}$ (7)

If A is a symmetric matrix then $\Psi \left(\omega \right)$ and $\Theta \left(\omega \right)$ are symmetric matrices and in this case, the SSOR iteration matrix has real spectrum. A relationship between the spectral radii of the three matrices $\Gamma \left(\omega \right)$ , $\Psi \left(\omega \right)$ and $\Theta \left(\omega \right)$ follows.

Lemma 1. $\rho \left(\Gamma \left(\omega \right)\right)=\rho \left(\Psi \left(\omega \right)\right)=\rho \left(\Theta \left(\omega \right)\right)$ , for all $\omega >0$ .

Proof. Since $\left(I-\omega U\right)\Gamma \left(\omega \right){\left(I-\omega U\right)}^{-1}=\left(\left(1-\omega \right)I+\omega L\right){\left(I-\omega L\right)}^{-1}\left(\left(1-\omega \right)I+\omega U\right){\left(I-\omega U\right)}^{-1}$ $=\Psi \left(\omega \right)$ , then $\Gamma \left(\omega \right)$ is similar to $\Psi \left(\omega \right)$ . Now $\Theta \left(\omega \right)$ and $\Psi \left(\omega \right)$ have the same spectral radius from the well known product theorem [14] , i.e. $\rho \left(AB\right)=\rho \left(BA\right)$ for any A and B. Hence, $\rho \left(\Gamma \left(\omega \right)\right)$ , $\rho \left(\Psi \left(\omega \right)\right)$ and $\rho \left(\Theta \left(\omega \right)\right)$ are equal.

Lemma 1 suggests using a different ordering for the SSOR iteration scheme.

Assuming that the SSOR iterative method converges for $\omega >0$ , i.e. $\rho \left(\Gamma \left(\omega \right)\right)<1$ with $\omega$ a positive acceleration factor, then the Kellogg-type SSOR iteration method is written as:

$\begin{array}{l}\left(I-\omega L\right)y^k=\left(\left(1-\omega \right)I+\omega L\right)x^k+\omega b_1\mathrm{,}\\ \left(I-\omega U\right)x^{k+1}=\left(\left(1-\omega \right)I+\omega U\right)y^k+\omega b_2\end{array}$ (8)

where $b=b_1+b_2$ . Now for $z=y^k+x^{k+1}$ in our trial solution, one tests $\Vert Az-b\Vert$ where $\Vert \cdot \Vert$ is the 2-norm. If this quantity fails to satisfy our convergence criteria, one assigns $x^{k+1}$ to $x^k$ , and reiterates.

Upon eliminating $y^k$ and the $x^k$ , from (8) one has

$\begin{array}{l}{x}^{k+1}=\Theta \left(\omega \right)x^k+{\left(I-\omega L\right)}^{-1}\left(\left(1-\omega \right)I+\omega L\right){\left(I-\omega U\right)}^{-1}\omega b_1+{\left(I-\omega U\right)}^{-1}\omega b_2\\ \text{and}\\ y^{k+1}=\Psi \left(\omega \right)y^k+{\left(I-\omega U\right)}^{-1}\left(\left(1-\omega \right)I+\omega U\right){\left(I-\omega L\right)}^{-1}\omega b_2+{\left(I-\omega L\right)}^{-1}\omega b_1\mathrm{.}\end{array}$ (9)

This decouples the original system (4).

The following lemma will be used to show the convergence of the Kellogg-type SSOR iteration method to a unique solution to the system of equations $Ax=b$ .

Lemma 2. For any $c_1\mathrm{,}{c}_2$ then there exists unique $y\mathrm{,}x$ satisfying

$\begin{array}{l}x={\left(I-\omega L\right)}^{-1}\left(\left(1-\omega \right)I+\omega L\right)y+c_1\mathrm{,}\\ y={\left(I-\omega U\right)}^{-1}\left(\left(1-\omega \right)I+\omega U\right)x+c_2\mathrm{.}\end{array}$ (10)

Furthemore,

$\begin{array}{l}y=\Psi \left(\omega \right)y+{\left(I-\omega L\right)}^{-1}\left(\left(I-\omega \right)I+\omega L\right)c_2+c_1,\\ x=\Theta \left(\omega \right)x+{\left(I-\omega U\right)}^{-1}\left(\left(I-\omega \right)I+\omega U\right)c_1+c_2.\end{array}$ (11)

Proof. The system of equations (10) may be represented by

$\tilde{A}\left[\begin{array}{c}x\\ y\end{array}\right]=\left[\begin{array}{c}{c}_1\\ c_2\end{array}\right]$ ,

where

$\begin{array}{c}\tilde{A}=\left[\begin{array}{cc}I& -{\left(I-\omega L\right)}^{-1}\left(\left(I-\omega \right)I+\omega L\right)\\ -{\left(I-\omega U\right)}^{-1}\left(\left(I-\omega \right)I+\omega U\right)& I\end{array}\right]\\ =\left[\begin{array}{cc}I& 0\\ -{\left(I-\omega U\right)}^{-1}\left(\left(I-\omega \right)I+\omega U\right)& I\end{array}\right]\left[\begin{array}{cc}I& 0\\ 0& I-\Theta \left(\omega \right)\end{array}\right]\\ \ast \left[\begin{array}{cc}I& -{\left(I-\omega L\right)}^{-1}\left(\left(I-\omega \right)I+\omega L\right)\\ 0& I\end{array}\right]\mathrm{.}\end{array}$

The det of $\tilde{A}$ is non-zero, completing the proof of the first part. The second part follows immediately.

The following theorem proves that the sequences $\left\{x^k\right\}$ and $\left\{y^k\right\}$ converge to x and y, respectively, and $x+y$ is the solution to $Ax=b$ .

Theorem 1. Let $x^k$ and $y^k$ be defined by (9) Then (i) ${\lim }_{k\to \infty }{x}^k=x$ , ${\lim }_{k\to \infty }{y}^k=y$ , and (ii) $A\left(x+y\right)=b$ .

Proof. Setting $e^n=x-x^n$ , $c_1={\left(I-\omega L\right)}^{-1}\omega b_1$ and $c_2={\left(I-\omega U\right)}^{-1}\omega b_2$ . Then $\begin{array}{c}y-y^n={\left(I-\omega L\right)}^{-1}\left(\left(I-\omega \right)I+\omega L\right)x+c_1-{\left(I-\omega L\right)}^{-1}\left(\left(I-\omega \right)I+\omega L\right)x^n-c_1\\ ={\left(I-\omega L\right)}^{-1}\left(\left(I-\omega \right)I+\omega L\right)x-{\left(I-\omega L\right)}^{-1}\left(\left(I-\omega \right)I+\omega L\right)x^n\end{array}$ . Since $x={\left(I-\omega U\right)}^{-1}\left(\left(I-\omega \right)I+\omega U\right)y+c_2$ and $x^n={\left(I-\omega U\right)}^{-1}\left(\left(I-\omega \right)I+\omega U\right)y^{n-1}+c_2$ , then ${\left(I-\omega L\right)}^{-1}\left(\left(I-\omega \right)I+\omega L\right)x-{\left(I-\omega L\right)}^{-1}\left(\left(I-\omega \right)I+\omega L\right)x^n$ becomes

$\Theta \left(\omega \right)y+{\left(I-\omega L\right)}^{-1}\left(\left(I-\omega \right)I+\omega L\right)c_2-\Theta \left(\omega \right)x^n-{\left(I-\omega L\right)}^{-1}\left(\left(I-\omega \right)I+\omega L\right)c_2$

In which case, $e^{n+1}=\Theta \left(\omega \right)e^n$ and since $\rho \left(\Theta \left(\omega \right)\right)<1$ , then ${\lim }_{k\to \infty }{x}^k=x$ . In a similar fashion, ${\lim }_{k\to \infty }{y}^k=y$ .

To show that $A\left(x+y\right)=b$ , set $c_1={\left(I-\omega L\right)}^{-1}\omega b_1$ and $c_2={\left(I-\omega U\right)}^{-1}\omega b_2$ by (i) and Lemma 2, the vector $\left[\begin{array}{c}x\\ y\end{array}\right]$ satisfies

$\begin{array}{l}y={\left(I-\omega L\right)}^{-1}\left(\left(1-\omega \right)I+\omega L\right)x+c_1\\ x={\left(I-\omega U\right)}^{-1}\left(\left(1-\omega \right)I+\omega U\right)y+c_2.\end{array}$ (12)

Rewriting this as

$\begin{array}{l}\left(I-\omega L\right)y=\left(\left(1-\omega \right)I+\omega L\right)x+\omega b_1\\ \left(I-\omega U\right)x=\left(\left(1-\omega \right)I+\omega U\right)y+\omega b_2\mathrm{,}\end{array}$ (13)

and then adding them together, canceling, and dividing by $\omega$ . This simplifies to

$\left(x+y\right)-L\left(x+y\right)-U\left(x+y\right)=\left(I-L-U\right)\left(x+y\right)=b_1+b_2$ as promised.

A flowchart for the Kellog-type SSOR iteration method follows, see Figure 1.

3. A Closer Look

For strictly lower and strictly upper triangular parts of an $n\times n$ matrix, the following holds.

Lemma 3. For L and U the strictly lower and strictly upper triangular parts of the $n\times n$ matrix A then

${\left(I-\omega L\right)}^{-1}\left(\left(1-\omega \right)I+\omega L\right)\text{and}{\left(I-\omega U\right)}^{-1}\left(\left(1-\omega \right)I+\omega U\right)$ (14)

are lower triangular and upper triangular matrices respectively.

Proof. Since $L^n$ is equal to the zero matrix and $\rho \left(\omega L\right)<1$ then ${\left(I-\omega L\right)}^{-1}$ exists and is given by ${\displaystyle \underset{i=0}{\overset{n-1}{\sum }}}{\left(\omega L\right)}^i$ .

Upon multiplying by $\left(\left(1-\omega \right)I+\omega L\right)$ gives the result that ${\left(I-\omega L\right)}^{-1}\left(\left(1-\omega \right)I+\omega L\right)$ is a lower triangular matrix.

Similarly, ${\left(I-\omega U\right)}^{-1}\left(\left(1-\omega \right)I+\omega U\right)$ is upper triangular.

One of the main drawbacks of the SSOR iteration method when compared to SOR iteration method is its two step iteration. Some benefits of SSOR iteration method for a symmetric matrix A, is $\Gamma \left(\omega \right)$ will have positive real spectrum and that $\rho \left(\Gamma \left(\omega \right)\right)$ is equal to the square of a norm of $\Delta \left(\omega \right)$ . See [9] pg. 237. In the KSSOR iteration method, if one takes $b_1=b$ and $b_2=0$ , then each step is a triangular matrix (one lower, in the first step of the iteration, while upper in the second step of the iteration) times a vector operation. Since the amount of

work done by KSSOR iteration may be considered less than the SSOR iteration method. Furthermore, since the spectral radius of the SSOR iteration matrix and the KSSOR iteration matrix are equal, their iteration counts should be fairly equal. However, in most instances the inverse would not be taken.

4. Cyclic Reduction

A Cyclic Reduction Algorithm, similar to [3] , and [8] pg. 170, is as follows.

Set $\tilde{A}=\left[\begin{array}{cc}I& -{\left(I-\omega L\right)}^{-1}\left(\left(I-\omega \right)I+\omega L\right)\\ -{\left(I-\omega U\right)}^{-1}\left(\left(I-\omega \right)I+\omega U\right)& I\end{array}\right]$ one sets $K=I-\tilde{A}=\left[\begin{array}{cc}0& {\left(I-\omega L\right)}^{-1}\left(\left(I-\omega \right)I+\omega L\right)\\ {\left(I-\omega U\right)}^{-1}\left(\left(I-\omega \right)I+\omega U\right)& 0\end{array}\right]$ and then $K^2=\left[\begin{array}{cc}\Theta \left(\omega \right)& 0\\ 0& \Psi \left(\omega \right)\end{array}\right]$ .

Lemma 4. If $\rho \left(\Gamma \left(\omega \right)\right)<1$ then $\rho \left(K^2\right)<1$ and $\rho \left(K\right)<1$ .

Proof. Since $\sigma \left(K^2\right)=\sigma \left(\Theta \left(\omega \right)\right)\cup \sigma \left(\Psi \left(\omega \right)\right)=\sigma \left(\Theta \left(\omega \right)\right)$ , where the last equality follows from the similarity of $\Theta \left(\omega \right)$ and $\Psi \left(\omega \right)$ then $\rho \left(K^2\right)<1$ . By the spectral mapping theorem $\pm \sqrt{\lambda }\in \sigma \left(K\right)$ if and only if $\lambda \in \sigma \left(K^2\right)$ ; i.e. $\text{determinant}\left(A^2-\lambda I\right)=\text{determinant}\left(A+\sqrt{\lambda }I\right)\text{determinant}\left(A-\sqrt{\lambda }I\right)$ . hence, $\rho \left(K\right)<1$ , see [14] . This completes the proof of the lemma.

$K^2$ defines the following uncoupled iteration scheme.

$\begin{array}{l}{x}^{m+2}=\Theta \left(\omega \right)x^m+k_1\\ y^{m+2}=\Psi \left(\omega \right)y^m+k_2\end{array}$ (15)

where $k_1={\left(I-\omega L\right)}^{-1}\left(\left(I-\omega \right)I+\omega L\right){\left(I-\omega U\right)}^{-1}{b}_2+{\left(I-\omega L\right)}^{-1}{b}_1$ and $k_2={\left(I-\omega U\right)}^{-1}\left(\left(I-\omega \right)I+\omega U\right){\left(I-\omega L\right)}^{-1}{b}_1+{\left(I-\omega U\right)}^{-1}{b}_2$ .

We now focus our attention on the solution of this reduced equation $x^{m+2}=\Theta \left(\omega \right)x^m+k_1$ .

Since $\rho \left(\Theta \left(\omega \right)\right)<1$ , then the Cyclic Reduction Scheme $x^m$ converges to x and having found the vector x in (10), form the vector $y={\left(I-\omega U\right)}^{-1}\left(\left(I-\omega \right)I+\omega U\right)x+{\left(I-\omega U\right)}^{-1}{b}_2$ and from which $x+y$ is a solution.

A flowchart for the Cyclic Reduction iteration method follows, see Figure 2.

5. Examples

In the following examples, a solution vector x is constructed as

$x\left(i\right)=\frac{i}{N}\sin \left(\frac{i\ast \pi }{6}\right)$ for $i=1,\cdots ,N$ and an initial guess is $x_0={\left[\mathrm{1,1,}\cdots \mathrm{,1}\right]}^{\text{T}}$ .

The error in each of our schemes will be measured by the 2-norm $\Vert Ay-b\Vert$ . where y is the resulting approximate answer and our desired error tolerance will be 10⁻⁶. Our table includes the $\Vert x-y\Vert$ i.e. the actual 2-norm difference between the true solution, x, and our approximate solution y, the iteration count, and the elapsed time in seconds for the work of the loop doing the iteration. ω is found through an iterative search, all decimals are rounded to two places and in the KSSOR iteration method and the Cyclic Reduction Method, $b_1=b$ and $b_2=0$ .

Example 1. Consider the system of linear equations with the n × n symmetric matrix which appears in [4] Set P to be the m × m matrix

$P=\left[\begin{array}{cccccccc}2& -1& 0& 0& 0& 0& 0& 0\\ -1& 2& -1& 0& 0& 0& 0& 0\\ 0& -1& 2& -1& 0& 0& \mathrm{...}& 0\\ \mathrm{.}& \mathrm{.}& \mathrm{.}& \mathrm{.}& \mathrm{.}& \mathrm{.}& & \\ 0& 0& 0& \mathrm{..}& & -1& 2& -1\\ 0& 0& 0& 0& \mathrm{..}& 0& -1& 2\end{array}\right]$

and $T=\left(I\otimes P+P\otimes I\right)-\left(P\otimes P\right)/6$ , and $A=D^{-1}\ast T$ where D is the diagonal of T. The resulting matrix A, is then n × n where $n=m^2$ . For $m=32$ , the n × n matrix is 1024 × 1024 and Table 1 summarizes our results.

Example 2. Consider the system of linear equations with 1024 × 1024 coefficient matrix which appears in [11] .

$T=\left[\begin{array}{cccccccc}1& 1& 0& 0& 0& 0& 0& 0\\ -1& 3& 2& 0& 0& 0& 0& 0\\ 0& -1& 5& 3& 0& 0& \mathrm{...}& 0\\ \mathrm{.}& \mathrm{.}& \mathrm{.}& \mathrm{.}& \mathrm{.}& \mathrm{.}& & \\ 0& 0& 0& \mathrm{..}& & -1& 2\ast N-3& N-1\\ 0& 0& 0& 0& \mathrm{..}& 0& -1& 2\ast N-1\end{array}\right]$

Then A becomes $D^{-1}\ast T$ where D is the diagonal of T. Table 2 summarizes our results.

Example 3 Consider the system of linear equations with n × n coefficient matrix

$T=\left[\begin{array}{cc}B& E\\ -E^{\prime }& 0.5\ast I\end{array}\right]$

where $B=\left[\begin{array}{cc}I\otimes C+C\otimes I& 0\\ 0& I\otimes C+C\otimes I\end{array}\right]$ , $E=\left[\begin{array}{c}I\otimes F\\ F\otimes I\end{array}\right]$ .

Where $C=\text{tridiag}\left(-\mathrm{1,2,}-1\right)$ , $F=h\ast \text{tridiag}\left(-\mathrm{1,1,0}\right)$ are m × m matrices with $h=\frac{1}{m+1}$ resulting in $n=3m^2$ . Then A becomes $D^{-1}\ast T$ where D is the diagonal of T and for $m=20$ is a 1200 × 1200 matrix. Table 3 summarizes our results.

Since the spectral radius of all iteration matrices is the same in all the iteration schemes, the convergence to a solution of the system of linear equations has about the same number of iterations for all the iteration convergence schemes considered in this paper.

6. Conclusion

The KSSOR algorithm requires an upper and lower triangular solution to a linear system which may be more accurate than the SSOR algorithm which on the surface requires full matrix solves at each step.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Aljuhani, A., Maturi, D. and Alshehri, H. (2022) Variational Iteration Method for Solving Boussinesq Equations Using Maple. Applied Mathematics, 13, 960-967.[CrossRef]
[2]	Hadjidimos, A. and Neumann, M. (1989) Convergence Domains of SSOR Method Associated with p-Cyclic Matrices. BIT Numerical Mathematics, 29, 311-320.[CrossRef]
[3]	Hageman, L. and Varga, R. (1964) Block Iterative Methods for Cyclically Reduced Matrix Equations. Numerische Mathematik, 6, 106-119.[CrossRef]
[4]	Lynch, R.E., Rice, J.R. and Thomas, D.R. (1965) Tensor Product Analysus of Alternating Direction Implicit Methods. Journal of the Society for Industrial and Applied Mathematics, 13, 995-1006.[CrossRef]
[5]	Neumaier, G. and Varga, R.S. (1984) Exact Convergence and Divergence Domains for the Symmetric Successive Overrelaxation Iterative (SSOR) Method Applied to H-Matrices. Linear Algebra and Its Applications, 58, 261-274.[CrossRef]
[6]	Ojarikre, H. and Mamadu, E. (2022) Time Discretized Variational Iteration Method for the Stochastic Volatility Process with Jumps. Advances in Pure Mathematics, 12, 693-700.[CrossRef]
[7]	Smotzer, T. and Buoni, J. (2023) Another HSS Iteration Method. Numerical Algorithms.[CrossRef]
[8]	Varga, R.S. (2000) Matrix Iterative Analysis. 2nd Edition, Springer.[CrossRef]
[9]	Young, D. (1971) Iterative Solutions of Large Linear Systems. Academic Press.
[10]	Young, D. (1971) On the Accelerated SSOR Method for Solving Large Linear Systems. Advances in Mathematics, 23, 215-271.[CrossRef]
[11]	Zhang, C.-Y., Xue, Z. and Luo, S. (2016) A Convergence Analysis of SOR Iterative Methods for Linear Systems with Weak H-Matrices. Open Mathematics, 14, 747-760.[CrossRef]
[12]	Sheldon, J. (1955) On the Numerical Solution of Elliptic Difference Equations. Mathematical Tables and Other Aids to Computation, 9, 101-112.[CrossRef]
[13]	Kellogg, R.B. (1963) Another Alternating-Direction Method. Journal of the Society for Industrial and Applied Mathematics, 11, 976-981.[CrossRef]
[14]	Halmos, P.R. (1974) Finite-Dimensional Vector Spaces. Springer-Verlag.[CrossRef]