Project: Iteration Methods

6. Project: Iteration Methods#

6.1. Introduction#

In this project, you will practice how to use Jacobi method and Gauss-Seidel to solve linear systems in MATLAB. To get started with this project, you will need to download the starter code and unzip its contents to the directory where you wish to complete the project.

Files included in this project

  • project6.m - MATLAB script that steps you through the project.

  • checkMyAnswer.m - MATLAB script that checks your code.

  • checkMyAnswer.mat - Data set.

  • data.mat - Data set.

  • * Jacobi.m - Jacobi iteration method.

  • * GS.m - Gauss-Seidel iteration method.

You need to complete the files with *.

Throughout the project, you will be using the script project6.m. This script set up all dataset for the problem and make calls to functions that you will write. Do not modify it. You are only required to modify functions in * files, following the instructions in this project.

6.2. Assignments#

Part 1 - Jacobi Iteration Method#

Recall in the Jacobi method to solve the linear system \(A\mathbf{x}=\mathbf{b}\), denote \(D\) the diagonal part of \(A\), then the iteration formula is

\[ \mathbf{x}_{k+1}=\mathbf{x}_{k}+D^{-1}\left(b-A\mathbf{x}_{k}\right) \]

You will implement the Jacobi method based on this form.

You should use the stop criteria as follows: stop the iteration when the latest iterates \(\mathbf{x}_{k+1}\) and \(\mathbf{x}_{k}\) satisfy

\[ \frac{||\mathbf{x}_{k+1}-\mathbf{x}_{k}||_\infty}{||\mathbf{x}_{k+1}||_\infty+eps} < tol \]

with a given error tolerance tol, or when the maximal number of iteration is reached. Here, eps is the machine epsilon, and it is used to avoid possible division by 0 when the relative error is computed. The quantity \(||\cdot||_\infty\) represents the infinity norm of a vector, which is defined as \(||\mathbf{x}||_\infty=\max_{1\le i\le n}{|x_i|}\). You can use the following MATLAB built-in function to calculate the infinity norm of the vector x.

norm(x,inf)

Hint

Do not forget to specify the inf argument when you use norm function in this project, as you need to use infinity norm in this project. If you omit inf argument, the function norm will use the default \(l_2\) norm, which is different.

Now you will find the outlines of MATLAB function in the file Jacobi.m. Modify it to return the solution of \(A\mathbf{x}=\mathbf{b}\) computed by Jacobi method.

Next you can run project6.m, and you should see the following output:

Part 1: Jacobi Iteration Method
Case 1: Diagonally Dominant
Matrix A =
     3     1    -1
     2     5    -1
    -3     2    10

Vector b =
     6
     6
    -8

Initial guess x0=
     0
     0
     0

Error tolerance delta = 1.00e-05.
Maximal number of iteration is 100.
Running Jacobi.m ...
A numerical solution satisfying the error tolerance is
found within 100 iterations.
The number of iterations used to compute x is 16.
Solution vector x =
    1.7333
    0.4333
   -0.3667

Case 2: Not Diagonally Dominant
Matrix A =
     3     4    -1
     2     1    -5
    -3     2     1

Vector b =
     1
     7
    -8

Initial guess x0=
     0
     0
     0

Error tolerance delta = 1.00e-05.
Maximal number of iteration is 100.
Running Jacobi.m ...
The program fails to produce a numerical solution in 100 iterations.
Solution vector x =
 -6.58e+49
 1.28e+50
 -8.99e+49

Part 2 - Gauss-Seidel Iteration Method#

For the Gauss-Seidel method, the iteration formula is

\[ \mathbf{x}_{k+1}=\mathbf{x}_{k}+E^{-1}\left(b-A\mathbf{x}_{k}\right) \]

where \(E\) is the lower triangular matrix consisting of corresponding elements of \(A\) and zeros elsewhere. Now in the file GS.m, you will find the outlines of MATLAB function. Modify it to return the solution of \(A\mathbf{x}=\mathbf{b}\) computed by Gauss-Seidel method.

Next you can continue running project6.m, and you should see the following output:

Part 2: Gauss-Seidel Iteration Method
Case 1: Diagonally Dominant
Matrix A =
     3     1    -1
     2     5    -1
    -3     2    10

Vector b =
     6
     6
    -8

Initial guess x0=
     0
     0
     0

Error tolerance delta = 1.00e-05.
Maximal number of iteration is 100.
Running Jacobi.m ...
A numerical solution satisfying the error tolerance
is found within 100 iterations.
The number of iterations used to compute x is 7.
Solution vector x =
    1.7333
    0.4333
   -0.3667

Case 2: Not Diagonally Dominant
Matrix A =
     3     4    -1
     2     1    -5
    -3     2     1

Vector b =
     1
     7
    -8

Initial guess x0=
     0
     0
     0

Error tolerance delta = 1.00e-05.
Maximal number of iteration is 100.
Running Jacobi.m ...
The program fails to produce a numerical solution in 100 iterations.
Solution vector x =
 -9.98e+64
 -8.00e+63
 -2.83e+65

6.3. Submit your code#

This is the last step. You should submit your code files to WTClass. After clicking Browse My Computer, you should select all files listed in Introduction to upload. Failure to do so will affect your grade.

Contents