ENGI 2211 Engineering Mathematics 2: Numerical Methods

1 Introduction

In this assignment,  you will learn more  about iterative solvers for linear systems. You will also learn about data structures for sparse matrices. Sparse matrices are matrices with only a small number of non-zero entries.  They are very common in Engineering since all finite difference and finite element methods generate sparse matrices.  The default matrix data type in Matlab is dense, however, sparse matrices are supported. Have a look at the Matlab documentation!

Iterative methods use much less memory than direct methods like backslash in Matlab, allowing then to solve much larger problems on the same hardware. In this assignment, you will learn three iterative solvers:  Jacobi, Gauss-Seidel and  SOR. They  are not covered in class,  nor  sparse matrices in  Matlab, which means that in order to be able to do this assignment, you have to learn something new on your own.  Part of the assignment is programming the iterative methods in Matlab.  This assignment consists of five MatLab codes that you have to implement and a report that you have to write.

2 Programs

This assignment requires that you implement five MatLab programs. In the next sections, the programs are described.

2.1 Jacobi method for dense matrices

For this part of the assignment, you have to write a MatLab function jacobi_dense. m to apply the Jacobi method to dense linear systems.

The MatLab function jacobi_dense.m must have the following input-output format

function  x  =  jacobi_dense(A,b,x0)

where:

❼ x - Computed approximated solution

❼ A - Linear system matrix

❼ b - Linear system right-hand-side

❼ x0 - Initial guess

Any change in the declaration of the function will result in zero marks being awarded for the corresponding part.

To help you to understand better, have a look to the provided example for jacobi_dense.m.

The code that you implement should be written in such a way that for ex- ample to solve the following linear system:

we use the following lines in Matlab:

A  =  [1,-2,2;-1,1,-1;-2,-2,1];

b  =  [1,2,3]’;

x0  =  zeros(3,1);

x  =  jacobi_dense(A,b,x0);

2.2 Gauss-Seidel method for dense matrices

For this part of the assignment, you have to write a MatLab function gauss_seidel_dense. m to apply the Gauss-Seidel method to dense linear systems.

The MatLab function gauss_seidel_dense.m must have the following input- output format

function  x  =  gauss_seidel_dense(A,b,x0)

where:

❼ x - Computed approximated solution

❼ A - Linear system matrix

❼ b - Linear system right-hand-side

❼ x0 - Initial guess

Any change in the declaration of the function will result in zero marks being awarded for the corresponding part.

To help you to understand better, have a look to the provided example for gauss_seidel_dense.m.

The code that you implement should be written in such a way that for ex- ample to solve the following linear system:

we use the following lines in Matlab:

A  =  [1,-2,2;-1,1,-1;-2,-2,1];

b  =  [1,2,3]’;

x0  =  zeros(3,1);

x  =  gauss_seidel_dense(A,b,x0);

2.3 SOR method for dense matrices

For this part of the assignment, you have to write a MatLab function sor_dense. m to apply the SOR method to dense linear systems.

The MatLab function sor_dense.m must have the following input-output format

function  x  =  sor_dense(A,b,x0,omega)

where:

❼ x - Computed approximated solution

❼ A - Linear system matrix

❼ b - Linear system right-hand-side

❼ x0 - Initial guess

❼ omega - Relaxation parameter

Any change in the declaration of the function will result in zero marks being awarded for the corresponding part.

To help you to understand better, have a look to the provided example for sor_dense.m.

The code that you implement should be written in such a way that for ex- ample to solve the following linear system:

we use the following lines in Matlab:

A  =  [1,-2,2;-1,1,-1;-2,-2,1];

b  =  [1,2,3]’;

x0  =  zeros(3,1);

x  =  sor_dense(A,b,x0,0.5);

2.4 Jacobi method for sparse matrices

For this part of the assignment, you have to write a MatLab function jacobi_sparse. m to apply the Jacobi method to sparse linear systems.

The MatLab function jacobi_sparse.m must have the following input-output format

function  x  =  jacobi_sparse(A,b,x0)

where:

❼ x - Computed approximated solution

❼ A - Linear system sparse matrix

❼ b - Linear system right-hand-side

❼ x0 - Initial guess

Any change in the declaration of the function will result in zero marks being awarded for the corresponding part.

2.5 Gauss-Seidel method for sparse matrices

For this part of the assignment, you have to write a MatLab function gauss_seidel_sparse. m to apply the Gauss-Seidel method to sparse linear systems.

The MatLab function gauss_seidel_sparse.m must have the following input- output format

function  x  =  gauss_seidel_sparse(A,b,x0)

where:

❼ x - Computed approximated solution

❼ A - Linear system sparse matrix

❼ b - Linear system right-hand-side

❼ x0 - Initial guess

Any change in the declaration of the function will result in zero marks being awarded for the corresponding part.

3 Report

The report should not be longer than five sides of A4 and the minimum allowed font size is 12pt.  Margins should be no smaller than 2 cm. Everything including tables and references has to be within the five sides of A4. The report should be in pdf format and named using the following convention:  SURNAME-Firstname_ENGI2211_NM.pdf. All the files have to be submitted in one Zip file.  The Zip file should contain the report in pdf format and the five codes.  The Zip file should be named using the following convention:  SURNAME-Firstname_ENGI2211_NM.zip.  In the report, you have to include the tables described below that you can fill in using the results from your programs. For the tables, report your answers to 3 significant figures.

You can use the report to explain the methods you have implemented.  In particular, you must include the following points in your report:

1.  Present the sparse matrix format in Matlab.

2.  Explain your implementation of the codes.

3.  Explain any technique used to make the code run faster.

4. Using the code jacobi_dense.m fill in Table 1 with the entries in the

solution vectors x computed using the matrix:

A  =  [1,-2,2;-1,1,-1;-2,-2,1];

and the three right-hand-sides:  r1  =  [1,2,3]’;

r2  =  [4,  -9,  7]’;

r3  =  [5,  -2,  4]’;

5. Using the code gauss_seidel_dense.m fill in Table 2 with the entries

in the solution vectors x computed using the matrix:

A  =  [2,-1,1;2,2,2;-1,-1,2];

and the three right-hand-sides:  r1  =  [1,2,3]’;

r2  =  [4,  -9,  7]’;

Table 1:  Computed the solutions using the three right-hand-sides  (answers to 3sf).

r3  =  [5,  -2,  4]’;

Table 2:  Computed the solutions using the three right-hand-sides  (answers to 3sf).

6.  Using the code  sor_dense.m fill in Table  3 with the entries in the

solution vectors x computed using the matrix:

A  =  [2,-1,1;2,2,2;-1,-1,2];

and the three right-hand-sides:  r1  =  [1,2,3]’;

r2  =  [4,  -9,  7]’;

r3  =  [5,  -2,  4]’;

and with ω = 0.5.

Table 3:  Computed the solutions using the three right-hand-sides  (answers to 3sf).

7. Using the codes jacobi_sparse.m and gauss_seidel_sparse.m fill in Table 4 with the norm the solution vectors x computed solving the linear system solved inside  fdelliptic. m using the following code: n  =  20;

f=@(x,y)(-1+0.*x);

[p,u]  =  fdelliptic(0,1,n,f);

Make sure to take into account how boundary conditions are imposed!

The linear system solved inside fdelliptic. m issparse but the routine treat it as dense.  Before passing the matrix to jacobi_sparse.m and gauss_seidel_sparse.m, convert it to sparse format.

Table 4: Computed the solutions using the right-hand-side from fdelliptic.m (answers to 3sf).

8.  Discuss engineering applications of iterative solvers.

Feel free to discuss any other aspect of your work that you consider interest- ing.