Linear Algebra – Matrices – Numerical methods

This subject is no longer in Lay, however it was in the first editions.

We consider two methods for iterative estimates of solutions of a system of linear equations. We suppose that a system of linear equations of the form \(A\mathbf{x}=\mathbf{b}\) is uniquely solvable (so \(A\) is invertible). Both methods are based on a recursive algorithm of the form

\[M\mathbf{x}_{k+1}=N\mathbf{x}_k+\mathbf{b},\quad k=0,1,2,\ldots,\]

where \(A=M-N\). If the sequence \(\{\mathbf{x}_k\}_{k=0}^{\infty}\) converges to a vector \(\mathbf{x}\), then we have:

\[M\mathbf{x}=N\mathbf{x}+\mathbf{b}\quad\Longleftrightarrow\quad(M-N)\mathbf{x}=\mathbf{b}\quad\Longleftrightarrow\quad A\mathbf{x}=\mathbf{b}.\]

Hence: \(\mathbf{x}\) is the unique solution of \(A\mathbf{x}=\mathbf{b}\).

The split up \(A=M-N\) can be made in several different ways. In both methods that we study we assume that the diagonal entries of \(A\) are all nonzero. In the Jacobi method \(M\) is a diagonal matrix with the diagonal entries of \(A\) on its main diagonal and in the Gauss-Seidel method \(M\) is a lower triangular matrix where the lower triangular part is equal to that of \(A\).

Problems that may arise:

1. The methods will not always work.


2. Sometimes the processes converge extremely slow or not at all.


3. How do we have to choose the initial vector \(\mathbf{x}_0\) such that the process converges (as rapid as possible)?


4. If one can choose, which method works best?

We consider an example first and will shortly comment on these problems afterwards.

Example: Consider \(A\mathbf{x}=\mathbf{b}\) with \(A=\begin{pmatrix}5&2&-1\\-1&5&2\\2&1&5\end{pmatrix}\) and \(\mathbf{b}=\begin{pmatrix}20\\5\\3\end{pmatrix}\).

The Jacobi method: \(M=\begin{pmatrix}5&0&0\\0&5&0\\0&0&5\end{pmatrix}\) and \(N=\begin{pmatrix}0&-2&1\\1&0&-2\\-2&-1&0\end{pmatrix}\). Hence:

\[M\mathbf{y}=N\mathbf{x}+\mathbf{b}\quad\Longleftrightarrow\quad\left\{\begin{array}{rcrrrrr}5y_1&=&&-2x_2&+x_3&+20\\[2.5mm] 5y_2&=&x_1&&-2x_3&+5\\[2.5mm]5y_3&=&-2x_1&-x_2&&+3.\end{array}\right.\]

Starting with \(\mathbf{x}_0=\mathbf{0}\) we now obtain: \(\mathbf{x}_1=\begin{pmatrix}4\\1\\3/5\end{pmatrix}\), \(\mathbf{x}_2=\begin{pmatrix}93/25\\39/25\\-6/5\end{pmatrix}\) and so on. In decimals:

\[\mathbf{x}_1=\begin{pmatrix}4.0\\1.0\\0.6\end{pmatrix},\quad\mathbf{x}_2=\begin{pmatrix}3.72\\1.56\\-1.20\end{pmatrix},\quad \mathbf{x}_3=\begin{pmatrix}3.136\\2.224\\-1.200\end{pmatrix},\quad\mathbf{x}_4=\begin{pmatrix}2.8704\\2.1072\\-1.0992\end{pmatrix}, \quad\mathbf{x}_5=\begin{pmatrix}2.93728\\2.01376\\-0.96960\end{pmatrix},\quad\ldots.\]

The Gauss-Seidel method: \(M=\begin{pmatrix}5&0&0\\-1&5&0\\2&1&5\end{pmatrix}\) and \(N=\begin{pmatrix}0&-2&1\\0&0&-2\\0&0&0\end{pmatrix}\). Hence:

\[M\mathbf{y}=N\mathbf{x}+\mathbf{b}\quad\Longleftrightarrow\quad\left\{\begin{array}{rrrcrrrrr}5y_1&&&=&&-2x_2&+x_3&+20\\[2.5mm] -y_1&+5y_2&&=&&&-2x_3&+5\\[2.5mm]2y_1&+y_2&+5y_3&=&&&&3.\end{array}\right.\]

Starting with \(\mathbf{x}_0=\mathbf{0}\) we now obtain: \(\mathbf{x}_1=\begin{pmatrix}4\\9/5\\-34/25\end{pmatrix}\), \(\mathbf{x}_2=\begin{pmatrix}376/125\\1341/625\\-3226/3125\end{pmatrix}\) and so on. In decimals:

\[\mathbf{x}_1=\begin{pmatrix}4.00\\1.80\\-1.36\end{pmatrix},\quad\mathbf{x}_2=\begin{pmatrix}3.00800\\2.14560\\-1.03232\end{pmatrix}, \quad\mathbf{x}_3=\begin{pmatrix}2.93529600\\1.99998720\\-0.97411584\end{pmatrix}, \quad\mathbf{x}_4=\begin{pmatrix}3.00518195200\\1.99068272640\\-1.00020932608\end{pmatrix},\quad\ldots.\]

The exact solution is \(\mathbf{x}=\begin{pmatrix}3\\2\\-1\end{pmatrix}\). We see that the Gauss-Seidel methode converges (quite) faster than the Jacobi method. This is often the case, however the Gauss-Seidel method is somewhat more complicared than the Jacobi method. It is also possible to construct examples where the Jacobi method is the fastest. It is alos possible to construct examples where one of the methods or even both methods will not converges.

However there exists an simple condition that guarantees that both methods will converge. This condition is not necessary (also in other cases we can have convergence), however in practice it is often possible to satisfy this condition. In that case the convergence is guaranteed for both methods. For this condition it is necessary that the matrix \(A\) is square. A square matrix \(A\) is called strictly diagonally dominant if the absolute value of each diagonal entry of \(A\) is larger than the sum of the absolute values of the other entries in the same row. In that case it can be shown that the matrix \(A\) is invertible (so \(A\mathbf{x}=\mathbf{b}\) is uniquely solvable) and that both the Jacobi method and the Gauss-Seidel method converge to the unique solution for every initial vector \(\mathbf{x}_0\). The proof of this statement is left out here.

Note that the matrix \(\begin{pmatrix}5&2&-1\\-1&5&2\\2&1&5\end{pmatrix}\) is strictly diagonally dominant, since: \(|5| > |2|+|\pm1|=3\).

Example: Consider \(A\mathbf{x}=\mathbf{b}\) with \(A=\begin{pmatrix}10&3&2\\3&10&2\\3&2&10\end{pmatrix}\) and \(\mathbf{b}=\begin{pmatrix}15\\15\\15\end{pmatrix}\).

It is easily seen that \(A\) is strictly diagonally dominant (since: \(|10| > |3| + |2| = 5\)) and that the unique solution is \(\mathbf{x}=\begin{pmatrix}1\\1\\1\end{pmatrix}\) is.

The Jacobi method: \(M=\begin{pmatrix}10&0&0\\0&10&0\\0&0&10\end{pmatrix}\) and \(N=\begin{pmatrix}0&-3&-2\\-3&0&-2\\-3&-2&0\end{pmatrix}\). Hence:

\[M\mathbf{y}=N\mathbf{x}+\mathbf{b}\quad\Longleftrightarrow\quad\left\{\begin{array}{rcrrrrr}10y_1&=&&-3x_2&-2x_3&+15\\[2.5mm] 10y_2&=&-3x_1&&-2x_3&+15\\[2.5mm]10y_3&=&-3x_1&-2x_2&&+15.\end{array}\right.\]

Starting with \(\mathbf{x}_0=\mathbf{0}\) we now obtain:

\[\mathbf{x}_1=\begin{pmatrix}1.5\\1.5\\1.5\end{pmatrix},\quad\mathbf{x}_2=\begin{pmatrix}0.75\\0.75\\0.75\end{pmatrix},\quad \mathbf{x}_3=\begin{pmatrix}1.125\\1.125\\1.125\end{pmatrix},\quad\mathbf{x}_4=\begin{pmatrix}0.9375\\0.9375\\0.9735\end{pmatrix}, \quad\mathbf{x}_5=\begin{pmatrix}1.03125\\1.03125\\1.03125\end{pmatrix},\quad\ldots.\]

The Gauss-Seidel method: \(M=\begin{pmatrix}10&0&0\\3&10&0\\3&2&10\end{pmatrix}\) and \(N=\begin{pmatrix}0&-3&-2\\0&0&-2\\0&0&0\end{pmatrix}\). Hence:

\[M\mathbf{y}=N\mathbf{x}+\mathbf{b}\quad\Longleftrightarrow\quad\left\{\begin{array}{rrrcrrrrr}10y_1&&&=&&-3x_2&-2x_3&+15\\[2.5mm] 3y_1&+10y_2&&=&&&-2x_3&+15\\[2.5mm]3y_1&+2y_2&+10y_3&=&&&&15.\end{array}\right.\]

Starting with \(\mathbf{x}_0=\mathbf{0}\) we now obtain:

\[\mathbf{x}_1=\begin{pmatrix}1.50\\1.05\\0.84\end{pmatrix},\quad\mathbf{x}_2=\begin{pmatrix}1.0170\\1.0269\\0.98952\end{pmatrix}, \quad\mathbf{x}_3=\begin{pmatrix}0.99402600\\1.00388820\\1.00101456\end{pmatrix},\quad \mathbf{x}_4=\begin{pmatrix}0.998630628000\\1.00020789960\\1.00036923168\end{pmatrix},\quad\ldots.\]

---

Last modified on April 12, 2021