Question

Write a program that solves the problem described in Exercise 10 , with a right-hand-side vector defined so that the solution is a vector of all \(1 \mathrm{~s}\). For example, for \(\omega=0\) the matrix and right-hand side can be generated by $$ \begin{aligned} &\mathrm{A}=\text { delsq }\left(\text { numgrid }\left(r \mathrm{~S}^{\prime}, \mathrm{N}+2\right)\right) ; \\ &\mathrm{b}=\mathrm{A} \text { *ones }\left(\mathrm{N}^{\wedge} 2,1\right) \end{aligned} $$ Your program will find the numerical solution using the Jacobi, Gauss-Seidel, SOR, and CG methods. For CG use the MATLAB command pcg and run it once without preconditioning and once preconditioned with incomplete Cholesky IC(0). For SOR, use the formula for finding the optimal \(\omega^{*}\) and apply the scheme only for this value. As a stopping criterion use \(\left\|\mathbf{r}_{k}\right\| /\left\|\mathbf{r}_{0}\right\|<10^{-6}\). Also, impose an upper bound of 2000 iterations. That is, if a scheme fails to satisfy the accuracy requirement on the relative norm of the residual after 2000 iterations, it should be stopped. For each of the methods, start with a zero initial guess. Your program should print out the following: \- Iteration counts for the five cases (Jacobi, Gauss-Seidel, SOR, CG, PCG). \- Plots of the relative residual norms \(\frac{\left\|\mathbf{r}_{k}\right\|}{\mathfrak{b} \|}\) vs. iterations. Use the MATLAB command semilogy for plotting. Use two grids with \(N=31\) and \(N=63\), and repeat your experiments for three values of \(\omega: \omega=0, \omega^{2}=10\), and \(\omega^{2}=1000 .\) Use your conclusions from Exercise 16 to explain the observed differences in speed of convergence.

Short answer

Question: Implement the iterative methods (Jacobi, Gauss-Seidel, SOR, CG, and preconditioned CG) to find the numerical solution to a given problem using the matrix A and the right-hand-side vector b, and the provided stopping criterion. Output the iteration counts and the relative residual norms for each method. Answer: To solve the given problem, we must first construct the matrix A and the vector b using prescribed formulas. Next, we will implement the Jacobi, Gauss-Seidel, SOR, CG, and preconditioned CG methods. We'll then apply the given stopping criterion and upper bound of 2000 iterations for each method. Finally, we'll output the iteration counts and plot the relative residual norms for each method.

Step-by-step solution

Construct the matrix A and the vector b

First, we need to generate the matrix A and the right-hand-side vector b using the given formula for different values of N and omega: $ \begin{aligned} &\mathrm{A}=\text { delsq }\left(\text { numgrid }\left(r \mathrm{~S}^{\prime}, \mathrm{N}+2\right)\right) ; \\\ &\mathrm{b}=\mathrm{A} \text { *ones }\left(\mathrm{N}^{\wedge} 2,1\right) \end{aligned} $ We will be using two values of N (31 and 63), and three values of omega (0, 10, and 1000).

Implement iterative methods

Next, we need to implement the Jacobi, Gauss-Seidel, SOR, and CG methods to find the numerical solution to the given problem. We will also apply the preconditioned CG method using the MATLAB command pcg and the incomplete Cholesky preconditioner IC(0).

Apply stopping criterion and upper bound of iterations

For each method, we need to apply the given stopping criterion and the upper bound of iterations. We must stop the iteration if the relative norm of the residual is less than \(10^{-6}\) or if the iteration count exceeds 2000.

Output iteration counts and residual norm plots

Finally, we need to output the iteration counts for the five cases (Jacobi, Gauss-Seidel, SOR, CG, PCG) and plot the relative residual norms (\(\frac{\left\|\mathbf{r}_{k}\right\|}{\mathfrak{b} \|\)) for each method using the MATLAB command semilogy. After completing these steps, we can analyze the results and observe the differences in the speed of convergence for each method. We can then use the conclusions from Exercise 16 to explain these differences.

Key concepts

Jacobi Method

The Jacobi Method is one of the simpler yet fundamental iterative techniques used to solve systems of linear equations. The key principle of this method involves decomposing the matrix into its diagonal and off-diagonal components. The iteration formula is derived by isolating each diagonal entry and progressively updating the solution based on values from the previous iteration. This method is straightforward to implement and requires minimal computational effort per iteration. Each variable in the linear system is updated simultaneously, relying solely on the approximate values from the preceding step.

The Jacobi Method is ideal for systems where the diagonal dominance condition is met. It ensures convergence, meaning the calculated solutions progressively approach the actual solution. Nevertheless, its convergence speed might not be the fastest, primarily if the system is large or poorly conditioned. Iterations continue until the difference between consecutive approximations falls below a specified tolerance level, ensuring the solution's accuracy.

To illustrate, for a matrix \( A \) and a right-hand-side vector \( b \), the Jacobi iteration can be expressed as:
\[x^{(k+1)} = D^{-1}(b - Rx^{(k)})\]
where \( D \) is the diagonal component, and \( R \) is the remainder of the matrix \( A \). Each update depends entirely on the values from the previous round, which is the hallmark of the Jacobi method.

Gauss-Seidel Method

The Gauss-Seidel Method builds upon the foundation of the Jacobi Method, offering an improvement in convergence speed. While the Jacobi Method updates all variables simultaneously using the previous values, the Gauss-Seidel Method makes updates in-place at each step in the iteration. As soon as a new value of a variable is computed, it is immediately used for calculating the subsequent variables. This in-place update helps accelerate convergence, ultimately requiring fewer iterations.

A significant advantage of the Gauss-Seidel Method is its efficiency for diagonally dominant matrices. In a linear system where this condition is met, the Gauss-Seidel Method often converges more quickly than the Jacobi Method. The Gauss-Seidel iteration for a matrix \( A \) and vector \( b \) can be written as:
\[x^{(k+1)}_i = \frac{1}{a_{ii}} \left( b_i - \sum_{j=1}^{i-1} a_{ij}x^{(k+1)}_j - \sum_{j=i+1}^{n} a_{ij}x^{(k)}_j \right)\]
Here, \( a_{ij} \) represents the entries of the matrix \( A \), and \( x_i \) represents the components of the solution vector.

Although an improvement over the Jacobi Method, the Gauss-Seidel approach still has challenges with ill-conditioned systems, requiring careful tuning of the stopping criteria to ensure reliable solutions.

SOR Method

The Successive Over-Relaxation (SOR) Method enhances the Gauss-Seidel Method by introducing a relaxation factor \( \omega \). This factor is used to accelerate the convergence by adjusting the update of the solution vector. The essence of the SOR Method lies in refining the iterative approach by fine-tuning the increment pace with which the new values approximate the true solution.

The relaxation factor \( \omega \) significantly impacts the method's effectiveness. When chosen optimally, \( \omega \) can lead to faster convergence than either the Jacobi or Gauss-Seidel methods. The iteration formula resembles:
\[x^{(k+1)}_i = (1-\omega)x^{(k)}_i + \frac{\omega}{a_{ii}} \left( b_i - \sum_{j=1}^{i-1} a_{ij}x^{(k+1)}_j - \sum_{j=i+1}^{n} a_{ij}x^{(k)}_j \right)\]
This formula shows how the new approximation \( x^{(k+1)}_i \) is computed based on a blend of the historical value and a partly updated state, controlled via \( \omega \). The optimal relaxation factor depends on the matrix characteristics and requires numerical experiments or specific techniques to determine.

The SOR Method is particularly useful in handling large sparse systems but requires careful calibration, as a poorly chosen \( \omega \) might slow convergence or even lead to divergence.

Conjugate Gradient Method

The Conjugate Gradient Method is a sophisticated iterative approach designed for solving sparse, symmetric positive-definite linear systems. Its strength lies in combining several mathematical ideas, including conjugacy and orthogonality, to achieve a quick convergence.

Whereas methods like Jacobi or Gauss-Seidel are more direct and straightforward, the Conjugate Gradient Method follows a strategic path that seeks to minimize the error in the solution space systematically. This strategy involves generating a series of search directions that are mutually orthogonal (conjugate) concerning the matrix. Starting from an initial guess, the method iterates to find the optimal path and reduce the residual error from one step to the next efficiently.

In practice, the Conjugate Gradient Method is very effective for large-scale problems where other methods struggle, both in terms of computation time and memory requirements. Its adaptability is further enhanced when combined with preconditioning techniques such as the incomplete Cholesky preconditioner to improve convergence with ill-conditioned systems. The preconditioned version can be expressed through the MATLAB function `pcg`, ensuring robustness and scalability for demanding problems.

This method is defined by its elegant approach and ability to adopt the underlying geometry of problem space, which helps it find solutions faster by reducing redundant calculations across the iterative process.