Assignment 2 - 2010 - Solution/Question 4

<rst> <rst-options: 'toc' = False/> <rst-options: 'reset-figures' = False/> Question 4 [3]

===

Write a MATLAB or Python function that implements the Jacobi (simultaneous update) method for solving a general system of equations, :math:`Ax = b`. The function should accept the following inputs:

• :math:`A`: any :math:`N \times N` matrix.
• :math:`b`: any :math:`N \times 1` vector.
• :math:`x_0`: any :math:`N \times 1` vector that is the initial guess for :math:`x`.
• :math:`\text{tol}`: a small number, so that we will stop the algorithm when the relative difference between successive iterations is small: i.e. stop when :math:`\displaystyle \frac{||x^{(k)} - x^{(k-1)}|| }{||x^{(k)}||} < \text{tol}`, where the norms are assumed to be the Euclidean norm.
• :math:`\text{maxiter}`: an integer, for example, 200, representing an upper bound on the number of iterations we are willing to perform.

The function does not have to perform pivoting to achieve diagonal dominance. The function should calculate and print out, at every iteration, the relative difference between the :math:`k^\text{th}` and :math:`k-1^\text{th}` iteration.

Consider this system of equations:

.. math::

y_F &= 0.5\\ G/L y_F - (\delta+1)x_1 + x_2 &= 0 \\ \delta x_1 - (\delta+1)x_2 + x_3 &= 0 \\ \delta x_2 - (\delta+1)x_3 + x_F &= 0 \\ x_F &= 0.1

where :math:`G/L = 1.6`, and :math:`\delta = 0.8`

Your answer should include:

1. . A print out of the function
2. . The relative difference between the estimates for the first 5 iterations of solving the system of **five** equations, start from :math:`y^{(0)}_F = x^{(0)}_1 = x^{(0)}_2 = x^{(0)}_3 = x^{(0)}_F = 0`.
3. . The number of iterations required, and the final solution, using a termination tolerance of :math:`\epsilon = 1 \times 10^{-9}`.
4. . What is the condition number of :math:`A`?

Solution

---

Given :math:`Ax = b`, the Jacobi method can be derived as shown in class, or an alternative derivation is given here, which leads to a slightly cleaner implementation. Let :math:`A = D + R` where D is a diagonal matrix containing diagonal elements of :math:`A`. Then :math:`x^{k+1}=D^{-1}(b-Rx^k)`. This procedure is generically implemented as shown in the MATLAB code below, and gives exactly the same solution as the Python code (the Python code implements the version shown in the course notes).

Note the termination criteria is for either exceeding the maximum iterations **or** achieving the desired tolerance. Now for the given system of 5 equations, we can simply call the function we have written, and provide appropriate input values as given in the problem:

.. twocolumncode:: :code1: ../che3e4/Assignments/Assignment-2/code/jacobi.m :language1: matlab :header1: MATLAB code :code2: ../che3e4/Assignments/Assignment-2/code/jacobi.py :language2: python :header2: Python code

The relative difference :math:`\displaystyle \frac{||x^{(k)} - x^{(k-1)}|| }{||x^{(k)}||}` and the first five iterations is given as::

            X-vector updates                                     Tolerance
   Iter 0 : [ 0.5,  0.        , -0.        , -0.        ,  0.1], 1.0
   Iter 1 : [ 0.5,  0.44444444, -0.        ,  0.05555556,  0.1], 0.65995
   Iter 2 : [ 0.5,  0.44444444,  0.22839506,  0.05555556,  0.1], 0.31895
   Iter 3 : [ 0.5,  0.57133059,  0.22839506,  0.15706447,  0.1], 0.19952
   Iter 4 : [ 0.5,  0.57133059,  0.34118275,  0.15706447,  0.1], 0.13224

It is seen that the relative difference is decreasing. The number of iterations required is 58 to achieve a relative error of :math:`8.604 \times 10^{-10}`, less than the required :math:`1 \times 10^{-9}`.

The solution at the 58th iteration is :math:`x=[0.500000, 0.695122, 0.451220, 0.256098, 0.100000]`, which agrees with the solution found using built-in functions.

The condition number of :math:`A` can be obtained by computing :math:`||A||||A^{-1}||` or using MATLAB's built-in function as cond(A). The result is :math:`||A||||A^{-1}||=8.7259`. Note that it is the 2-norm condition number.

</rst>