A.Kaw矩阵代数初步学习笔记 8. Gauss-Seidel Method

“矩阵代数初步”(Introduction to MATRIX ALGEBRA)课程由Prof. A.K.Kaw(University of South Florida)设计并讲授。
PDF格式学习笔记下载（Academia.edu）
第8章课程讲义下载（PDF）

Summary

• Algorithm
  Given a general set of $n$ equations and $n$ unknowns $$\begin{cases}a_{11}x_1 + a_{12}x_2 +\cdots + a_{1n}x_n = c_1 \\ a_{21}x_1 + a_{22}x_2 +\cdots + a_{2n}x_n = c_2 \\ \vdots\\ a_{n1}x_1 + a_{n2}x_2 +\cdots + a_{nn}x_n = c_n \\\end{cases}$$ If the diagonal elements are non-zero, each equation is rewritten for the corresponding unknown, that is, $$\begin{cases}x_1 = \displaystyle{c_1-a_{12}x_2 - a_{13}x_3-\cdots -a_{1n}x_n\over a_{11}}\\ x_2 = \displaystyle{c_2-a_{21}x_1 - a_{23}x_3-\cdots -a_{2n}x_n\over a_{22}}\\ \vdots\\ x_n = \displaystyle{c_n-a_{n1}x_1 - a_{n2}x_2-\cdots -a_{n,n-1}x_{n-1}\over a_{11}}\\ \end{cases}$$ $$\Rightarrow \begin{cases}x_1 = \displaystyle{c_1-\displaystyle\sum_{j=1,j\neq1}^{n}a_{1j}x_j\over a_{11}}\\ x_2 = \displaystyle{c_1-\displaystyle\sum_{j=1,j\neq2}^{n}a_{2j}x_j\over a_{22}}\\ \vdots\\ x_n = \displaystyle{c_n-\displaystyle\sum_{j=1,j\neq n}^{n}a_{nj}x_j\over a_{nn}}\\ \end{cases}$$ Hence for any row $i$, $$x_i = {c_i-\displaystyle\sum_{j=1,j\neq i}^{n}a_{ij}x_j\over a_{ii}}$$ where $i=1$, $2$, $\cdots$, $n$.
• Iteration
  To find $x_i$, we assume an initial guess for the $x_i$ and then use the rewritten equations to calculate the new estimates. We always use the most recent estimates to calculate the next estimates, $x_i$. At the end of each iteration, we calculate the absolute relative approximate error for each $x_i$ as $$\varepsilon_i = \left|{x_i^{\text{new}} - x_i^{\text{old}}\over x_i^{\text{new}}}\right|$$ where $x_i^{\text{new}}$ is the recently obtained value of $x_i$, and $x_i^{\text{old}}$ is the previous value of $x_i$. When the absolute relative approximate error for each $x_i$ is less than the pre-specified tolerance, the iterations are stopped.
• Convergent
  The coefficient matrix $[A]$ in $[A][X]=[B]$ must be diagonally dominant, that is, $$\begin{cases}|a_{ii}| \geq \displaystyle\sum_{j=1, j\neq i}^{n}a_{ij}&\text{for all}\ i\\ |a_{ii}| > \displaystyle\sum_{j=1, j\neq i}^{n}a_{ij} & \text{for at least one}\ i\end{cases}$$
• An example
  Suppose the following system of equations $$\begin{cases}12x_1 + 3x_2 -5x_3 =1\\ x_1 + 5x_2 +3x_3 =28\\ 3x_1 + 7x_2 +13x_3 =76\end{cases}$$ Use $$\begin{bmatrix}x_1\\ x_2\\ x_3 \end{bmatrix} = \begin{bmatrix}1\\ 0\\ 1\end{bmatrix}$$ as the initial guess and conduct two iterations.
  • Diagonally dominant test: $$\begin{cases}|a_{11}|=12 > |a_{12}| + |a_{13}| = 3+5=8\\ |a_{22}|=5 > |a_{21}| + |a_{23}|=1+3=4\\ |a_{33}|=13 > |a_{31}|+|a_{32}| = 3+7=10 \end{cases}$$ Hence the solution should converge using Gauss Seidel method.
  • Rewriting the equations: $$\begin{cases}x_1 = \displaystyle{1-3x_2 +5x_3 \over 12}\\ x_2 = \displaystyle{28-x_1-3x_3\over 5}\\ x_3= \displaystyle{76-3x_1-7x_2\over 13}\end{cases}$$ And the initial value is $$\begin{bmatrix}x_1\\ x_2\\ x_3 \end{bmatrix} = \begin{bmatrix}1\\ 0\\ 1\end{bmatrix}$$
  • Iteration 1: $$\begin{cases}x_1 = \displaystyle{1 - 3\times0 + 5\times1 \over 12} = 0.5\\ x_2 = \displaystyle{28 - 0.5 - 3\times1 \over 5} = 4.9\\ x_3= \displaystyle{76-3\times 0.5-7\times 4.9\over 13}=3.0923\end{cases}$$ Notice that the second and the third equations above, $x_1$ and $x_2$ are updated immediately. And the absolute relative approximate error is $$\begin{cases}\varepsilon_1 = \displaystyle{|0.5-1|\over0.5} = 1\\ \varepsilon_2 = \displaystyle{|4.9 - 0|\over 4.9} = 1\\ \varepsilon_3 = \displaystyle{|3.0923 - 1|\over3.0923}=0.67662 \end{cases}$$
  • Iteration 2: $$\begin{cases}x_1 = \displaystyle{1 - 3\times4.9 + 5\times3.0923 \over 12} = 0.14679\\ x_2 = \displaystyle{28 - 0.14679 - 3\times3.0923 \over 5} = 3.7153\\ x_3= \displaystyle{76-3\times 0.14679-7\times 3.7153\over 13} = 3.8118 \end{cases}$$ And the absolute relative approximate error is $$\begin{cases}\varepsilon_1 = \displaystyle{|0.14679-0.5|\over0.14679} = 2.4\\ \varepsilon_2 = \displaystyle{|3.7153 - 4.9|\over 3.7153} = 0.31889\\ \varepsilon_3 = \displaystyle{|3.8118 - 3.0923|\over3.8118}=0.18874 \end{cases}$$
  • Final result:
    After 6 iterations, we have the solution $$\begin{bmatrix}x_1\\ x_2\\ x_3 \end{bmatrix} = \begin{bmatrix}0.99919\\ 3.0001\\ 4.0001 \end{bmatrix}$$ which is very close to the exact solution $$\begin{bmatrix}x_1\\ x_2\\ x_3 \end{bmatrix} = \begin{bmatrix}1\\ 3\\ 4 \end{bmatrix}$$
• R code:
  Some comments:
  • In the second function PrepA, we use the elementary row operation $R_i + mR_j$ if the diagonal element in $R_i$ equals to zero.
  • In the third function IterSolve, we use $$\varepsilon=\sum\left|x^{\text{new}}-x^{\text{old}}\right|$$ instead of the absolute relative approximate error.
  • x0 is a vector in the main function, which is the initial guess of the system. And eps is the tolerance of the error, which can be smaller or bigger in different cases. The last parameter is maxit is the number of iterations, it does not need to be too much in most cases.
  • Using this code to calculate the previous example:

    A = matrix(c(12, 1, 3, 3, 5, 7, -5, 3, 13), ncol = 3)
    b = matrix(c(1, 28, 76), ncol = 1)
    IterSolve(A, b, c(1, 0, 1))$x
    # Result
    # Converged after  11 iterations
    # [1] 1 3 4

Selected Problems

1. Given the system of equations $$\begin{cases}3x_1 + 7x_2 + 13x_3 =76\\ x_1 +5x_2 + 3x_3 =28\\ 12x_1 +3x_2 -5x_3 =1\end{cases}$$ find the solutions using the Gauss-Seidel method. Use $$\begin{bmatrix}x_1\\ x_2\\ x_3 \end{bmatrix} = \begin{bmatrix}1\\ 0\\ 1\end{bmatrix}$$ as the initial guess.

Solution: Note that the coefficient matrix is not diagonal dominant: $$\begin{cases}|a_{11}| = 3 < |a_{12}| + |a_{13}| = 7+13 =20\\ |a_{33}| = 5 < |a_{31}| + |a_{32}| = 12 + 3 = 15\end{cases}$$ Hence it may diverge. Moreover, we can use our R code to test it:

A = matrix(c(3, 1, 2, 7, 5, 3, 13, 3, -5), ncol = 3)
b = matrix(c(76, 28, 1), ncol = 1)
IterSolve(A, b, c(1, 0, 1))$x
# Result
# [1] -2.496896e+172  1.261843e+171 -9.230477e+171
# Warning message:
# In IterSolve(A, b, c(1, 0, 1)) : Maxit reached

2. Solve the following system equations using Gauss-Seidel method. $$\begin{cases}12x_1 + 7x_2 + 3x_3 = 17\\ 3x_1 + 6x_2 +2x_3 =9\\ 2x_1 + 7x_2 -11x_3 =49\end{cases}$$ Choose the initial guess as $$\begin{bmatrix}x_1\\ x_2\\ x_3 \end{bmatrix} = \begin{bmatrix}1\\ 3\\ 5\end{bmatrix}$$

Solution:

Firstly, we test whether the coefficient matrix is diagonal dominant: $$\begin{cases}|a_{11}| = 12 > |a_{12}| + |a_{13}| =10\\ |a_{22}|=6 > |a_{21}| +|a_{23}|=5\\ |a_{33}| = 11 > |a_{31}| +|a_{32}|= 9\end{cases}$$ which means it is diagonal dominant. Then we will conduct two iterations: $$I_1=\begin{cases}x_1 = \displaystyle{17-7x_2-3x_3\over12} = {17-7\times3 -3\times5\over12} = -1.583333\\ x_2 = \displaystyle{9-3x_1-2x_3\over6} = {9-3\times(-1.583333) -2\times5\over6} = 0.625000\\ x_3=\displaystyle{49-2x_1-7x_2\over -11} = {49-2\times(-1.583333) -7\times 0.625000\over -11} = -4.344697\end{cases}$$ $$I_2=\begin{cases}x_1 = \displaystyle{17-7x_2-3x_3\over12} = \displaystyle{17-7\times0.625000 -3\times(-4.344697)\over12} = 2.138258\\ x_2 = \displaystyle{9-3x_1-2x_3\over6} = {9-3\times2.138258-2\times(-4.344697)\over6} = 1.879104\\ x_3=\displaystyle{49-2x_1-7x_2\over -11} = {49-2\times2.138258 -7\times 1.879104\over -11} = -2.869978\end{cases}$$ Alternatively, we can use R code to solve it directly:

A = matrix(c(12, 3, 2, 7, 6, 7, 3, 2, -11), ncol = 3)
b = matrix(c(17, 9, 49), ncol = 1)
IterSolve(A, b, c(1, 3, 5), eps = 1e-8)$x
# Result
# Converged after  16 iterations
# [1]  1  2 -3

That is, the solution is $$\begin{bmatrix}x_1\\ x_2\\ x_3 \end{bmatrix} = \begin{bmatrix}1\\ 2\\ -3\end{bmatrix}$$

3. Solve the following system equations using Gauss-Seidel method. $$\begin{cases}3x_1 + 6x_2 + 2x_3 =9\\ 12x_1 +7x_2 + 3x_3=17\\ 2x_1 +7x_2 -11x_3=49\end{cases}$$ Choose the initial guess as $$\begin{bmatrix}x_1\\ x_2\\ x_3 \end{bmatrix} = \begin{bmatrix} 1.1\\ 2.1\\ -2.9 \end{bmatrix}$$

Solution:

We will use the R code to solve it directly:

A = matrix(c(3, 12, 2, 6, 7, 7, 2, 3, -11), ncol = 3)
b = matrix(c(9, 17, 49), ncol = 1)
IterSolve(A, b, c(1, 0, 1), eps = 1e-8)$x
# Result
# Error in IterSolve(A, b, c(1, 3, 5)) : The algorithm diverges

Recall the R function, the result is divergent when the solution in the iterations goes to infinity. Moreover, we can read off its non-convergent according to it is not diagonal dominant since $$\begin{cases}|a_{11}|=3 < |a_{12}|+|a_{13}|=8\\ |a_{22}|=7 < |a_{21}|+|a_{23}|=15\end{cases}$$