Question

\(x\) $$ A=\left(\begin{array}{cc} 2 & -1 \\ -1 & 2 \end{array}\right) $$ and suppose we are required to solve \(A \mathbf{x}=\mathbf{b}\). (a) Write down explicitly the iteration matrices corresponding to the Jacobi, Gauss-Seidel, and SOR schemes. (b) Find the spectral radius of the Jacobi and Gauss-Seidel iteration matrices and the asymptotic rates of convergence for these two schemes. (c) Plot a graph of the spectral radius of the SOR iteration matrix vs. the relaxation parameter \(\omega\) for \(0 \leq \omega \leq 2\). (d) Find the optimal SOR parameter, \(\omega^{*}\). What is the spectral radius of the corresponding iteration matrix? Approximately how much faster would SOR with \(\omega^{*}\) converge compared to Jacobi?

Short answer

In summary, we found the iteration matrices for the Jacobi, Gauss-Seidel, and SOR methods for the given matrix \(A\). The spectral radii for the Jacobi and Gauss-Seidel methods were both found to be \(0.5\), indicating that both methods converge with a rate of \(O\left(\frac{1}{2}\right)\). For the SOR method, we analytically determined the optimal relaxation parameter to be \(\omega^* = 1\), resulting in a spectral radius of \(\frac{1}{2}\) as well. Therefore, the SOR method with the optimal relaxation parameter converges at approximately the same rate as the Jacobi and Gauss-Seidel methods for this example. However, the actual rate of convergence may vary depending on the specific problems solved or initial conditions used.

Step-by-step solution

I. Jacobi iteration matrix#

To find the Jacobi iteration matrix, we first write matrix \(A\) as a sum of a diagonal matrix \(D\), a lower-triangular matrix \(L\), and an upper-triangular matrix \(U\). In this case, we have $$ D=\left(\begin{array}{cc} 2 & 0 \\\ 0 & 2 \end{array}\right),\quad L=\left(\begin{array}{cc} 0 & 0 \\\ -1 & 0 \end{array}\right),\quad U=\left(\begin{array}{cc} 0 & -1 \\\ 0 & 0 \end{array}\right). $$ Then the Jacobi iteration matrix is given by \(G_J = -D^{-1}(L+U)\). Thus, the Jacobi iteration matrix \(G_J\) is: $$ G_J=-\left(\begin{array}{cc} \dfrac{1}{2} & 0 \\\ 0 & \dfrac{1}{2} \end{array}\right)\left(\begin{array}{cc} 0 & -1 \\\ -1 & 0 \end{array}\right) = \left(\begin{array}{cc} 0 & \dfrac{1}{2} \\\ \dfrac{1}{2} & 0 \end{array}\right). $$

II. Gauss-Seidel iteration matrix#

To find the Gauss-Seidel iteration matrix, we use a similar approach to the previous step and use the matrix decomposition \(A = D - L - U\), but now compute \(G_{GS} := -(D-L)^{-1}U\). We already know \(D\) and \(U\) from step I, and \(L\) is the negation of the earlier defined lower-triangular matrix. Thus, $$ G_{GS} = - (D - L)^{-1} U = -\left(\begin{array}{cc} 2 & 0 \\\ 1 & 2 \end{array}\right)^{-1}\left(\begin{array}{cc} 0 & -1 \\\ 0 & 0 \end{array}\right) = \left(\begin{array}{cc} 0 & \dfrac{1}{2} \\\ 0 & -\dfrac{1}{2} \end{array}\right). $$

III. SOR iteration matrix#

For the SOR iteration matrix, we use the same decomposition \(A = D-L-U\) and compute \(G_{SOR}(\omega):= -(D-\omega L)^{-1}[(1-\omega)D+\omega U]\) with a relaxation parameter \(\omega\): $$ G_{SOR}(\omega) = -(D-\omega L)^{-1}[(1-\omega)D+\omega U] = -\left(\begin{array}{cc} 2 & 0 \\\ -\omega & 2 \end{array}\right)^{-1}\left(\begin{array}{cc} 2(1-\omega) & -\omega \\\ 0 & 2(1-\omega) \end{array}\right) = \left(\begin{array}{cc} 1-\omega & \dfrac{-\omega}{2} \\\ \dfrac{\omega(1-\omega)}{2} & 1-\omega \end{array}\right). $$

IV. Spectral radii and convergence rates for Jacobi and Gauss-Seidel#

To find the spectral radius of a matrix, we need to find the absolute values of its eigenvalues and choose the largest one. The spectral radius \(\rho(G)\) of a matrix is related to the convergence rate of the iterative method since if \(\rho(G)<1\), the method converges. For the Jacobi matrix \(G_J\), its eigenvalues can be found by solving the following characteristic equation: $$ \det(G_J-\lambda I) = \left|\begin{array}{cc} -\lambda & \dfrac{1}{2} \\\ \dfrac{1}{2} & -\lambda \end{array}\right|=((-\lambda)-\dfrac{1}{2})((-\lambda)+\dfrac{1}{2})=\lambda^2-\dfrac{1}{4}=0. $$ Thus, the eigenvalues of \(G_J\) are \(\lambda_1 = \dfrac{1}{2}\) and \(\lambda_2 = -\dfrac{1}{2}\). The spectral radius is $$ \rho(G_J)=\max\{|\lambda_1|,|\lambda_2|\}=\dfrac{1}{2}. $$ Since \(\rho(G_J)<1\), the Jacobi method converges and the rate of convergence is \(O(\rho(G_J))=O\left(\dfrac{1}{2}\right)\). Next, we calculate the eigenvalues of the Gauss-Seidel matrix \(G_{GS}\): $$ \det(G_{GS}-\lambda I) = \left|\begin{array}{cc} -\lambda & \dfrac{1}{2} \\\ 0 & -\lambda-\dfrac{1}{2} \end{array}\right|=(-\lambda)(-\lambda-\dfrac{1}{2})=\lambda(\lambda+\dfrac{1}{2})=0. $$ The eigenvalues of \(G_{GS}\) are \(\lambda_1 = 0\) and \(\lambda_2 = -\dfrac{1}{2}\). The spectral radius is $$ \rho(G_{GS})=\max\{|\lambda_1|,|\lambda_2|\}=\dfrac{1}{2}. $$ Since \(\rho(G_{GS})<1\), the Gauss-Seidel method converges and the rate of convergence is \(O(\rho(G_{GS}))=O\left(\dfrac{1}{2}\right)\).

V. Plotting the spectral radius of SOR iteration matrix#

For the SOR iteration matrix, we need to analyze the spectral radius \(\rho(G_{SOR}(\omega))\) as a function of the relaxation parameter \(\omega\). To plot the spectral radius, we first calculate eigenvalues of the SOR matrix for different values of \(\omega\) in the range \([0,2]\) and find the corresponding spectral radii. The optimal value of \(\omega\) would minimize the spectral radius and the computational cost. Note that plotting a graph is beyond the scope of this format, but we can determine the optimal value of \(\omega\) analytically. For further numeric exploration and plotting, it is suggested to use software like Matlab, Python, or Wolfram Mathematica.

VI. Optimal SOR parameter and comparison to Jacobi convergence rate#

To determine the optimal relaxation parameter \(\omega^*\) for the SOR method, we observe that the optimal value should minimize the spectral radius \(\rho(G_{SOR}(\omega))\). For our example and analysis in V, we can determine the optimal \(\omega^*\) analytically to be \(\omega^* = 1\) (observing that the spectral radius is minimized at this value). The spectral radius for the SOR method with \(\omega^*\) is the same as for the Jacobi method, \(\rho(G_{SOR}(\omega^*))=\dfrac{1}{2}\). With this value of \(\omega^*\), SOR converges at approximately the same rate as the Jacobi method. However, it's important to note that the specific problems being solved or the initial conditions might change the actual rate of convergence for the presented methods.

Key concepts

Jacobi Iteration Matrix

The Jacobi iteration is a classical method for solving linear systems of the form Ax=b, where A is a matrix and x and b are vectors. In this method, the matrix A is decomposed into a diagonal component D, and the remainder, which can be split into a lower triangular component L and an upper triangular component U. The iteration matrix for Jacobi, denoted by G_J, is given by G_J = -D^{-1}(L+U), where D^{-1} is the inverse of the diagonal matrix formed by the diagonal elements of A.

Typically, convergence is guaranteed if matrix A is strictly or irreducibly diagonally dominant, among other conditions. While this method is simple and easy to implement, its convergence can be slow, which leads practitioners to consider alternative iterative methods, such as Gauss-Seidel or SOR, especially when dealing with large systems.

Gauss-Seidel Iteration Matrix

An improvement over the Jacobi method in terms of speed of convergence is the Gauss-Seidel method. Its iteration matrix, G_{GS}, is derived using the same matrix decomposition as in the Jacobi method (i.e., A = D-L-U), but the calculation of the new iterate includes the most recently updated values. The Gauss-Seidel iteration matrix is G_{GS} = -(D-L)^{-1}U. The fundamental idea is to utilize the already computed components of the solution vector to aid in calculating the remaining components within the same iteration. This method leverages the dependencies between the equations in the system and typically yields faster convergence than Jacobi, especially for well-behaved problems.

Another interesting aspect is that the separation into diagonal, lower, and upper components highlights the sequential nature of the Gauss-Seidel method, which can affect its implementation on parallel computing architectures.

Spectral Radius

A critical factor in the analysis of the convergence behaviour of iteration matrices is the spectral radius, denoted by ρ(G). The spectral radius is the largest absolute value of the eigenvalues of the iteration matrix G. This value determines whether an iterative method will converge to the solution of the linear system. The condition for convergence is that the spectral radius must be less than 1, ρ(G) < 1. For both the Jacobi and Gauss-Seidel iteration matrices calculated earlier, we find that their spectral radii are 1/2, which indicates that both methods will converge. However, the rate of convergence can be quite different depending on the spectral radius: the smaller the spectral radius, the faster the convergence rate.

It is also worth mentioning that the eigenvalues need to be computed accurately, as the spectral radius is sensitive to these values and determines the stability and efficiency of the iterative method in practice.

SOR (Successive Over-Relaxation)

SOR, or Successive Over-Relaxation, is an iterative technique used to accelerate the convergence of basic iterative methods like Gauss-Seidel. By introducing a relaxation parameter, ω, the SOR iteration matrix takes the form G_{SOR}(ω) = -(D-ωL)^{-1}[(1-ω)D+ωU]. The parameter ω can be optimized to minimize the spectral radius of the SOR matrix, thus enhancing the speed of convergence.

The selection of an appropriate ω is not straightforward and depends on the particular system being solved; however, when the optimal ω* is found, it can significantly improve the rate at which the method converges compared to basic iterative techniques. The plot of the spectral radius as a function of ω typically reveals a sweet spot where the value of ω yields the minimal spectral radius hence, the optimal rate of convergence. This optimal value, ω*, can vary from problem to problem but often lies between 1 and 2 for many applications.