Assignment 2 - 2010 - Solution/Question 4

From Process Model Formulation and Solution: 3E4
Jump to: navigation, search
← Question 3 (previous step) Back to all questions Next step: Question 5 →

Question 4 [3]

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

  • \(A\): any \(N \times N\) matrix.
  • \(b\): any \(N \times 1\) vector.
  • \(x_0\): any \(N \times 1\) vector that is the initial guess for \(x\).
  • \(\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 \(\displaystyle \frac{||x^{(k)} - x^{(k-1)}|| }{||x^{(k)}||} < \text{tol}\), where the norms are assumed to be the Euclidean norm.
  • \(\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 \(k^\text{th}\) and \(k-1^\text{th}\) iteration.

Consider this system of equations:

\[\begin{split}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\end{split}\]

where \(G/L = 1.6\), and \(\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 \(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 \(\epsilon = 1 \times 10^{-9}\).
  4. What is the condition number of \(A\)?

Solution

Given \(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 \(A = D + R\) where D is a diagonal matrix containing diagonal elements of \(A\). Then \(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:

MATLAB code Python code
% This code goes in a file called: jacobi.m

function [x, rel_error] = jacobi(A, b, x0, tol, maxiter)

% Implements the Jacobi method to solve Ax = b
% A = D + R where D is a diagonal matrix containing 
% diagonal elements of A.
%
% NOTE: this an alternative algorithm to what was presented
% ----  in the course notes.  However, the course notes
%       algorithm gives the same values of x at each 
%       iteration.

D = diag(diag(A));      % construct the diagonal matrix
R = A - D;              % Construct R
iter = 1;               % iteration counter
%x = inv(D)*(b - R*x0);  % Perform the first iteration
x = x0;
rel_error = tol * 2;% norm(x - x0)/norm(x);
%fprintf('\n Iteration %i: relative error =%d ',iter, rel_error)
%iter = iter + 1;
while (rel_error>tol && iter <= maxiter)
    
    xprev = x;    
    x = inv(D)*(b - R*xprev);
    rel_error = norm(x - xprev)/norm(x);
    
    fprintf('\n Iteration %i: Relative error =%d ',iter, rel_error)    
    iter = iter + 1;        
end

% The rest of this code goes in a file called anything else,
% for example, question4_assign2.m - it uses the jacobi.m file
GL = 1.6;
d = 0.8;

A = [1.0,      0,      0,      0,   0; ...
      GL, -(d+1),    1.0,      0,   0; ...
       0,      d, -(d+1),    1.0,   0; ...
       0,      0,      d, -(d+1), 1.0; ...
       0,      0,      0,      0, 1.0];
b = [0.5,      0,      0,      0, 0.1]';

% Solver parameters
tol = 1e-9;
maxiter = 200;
x0 = zeros(5,1);

fprintf('\nCondition number of A=%d\n',cond(A))

% Solve the system, using our jacobi.m function file
[x, rel_error] = jacobi(A, b, x0, tol, maxiter);

fprintf('\nJacobi function: x=[%f %f %f %f %f]\n rel error=%d',...
    x, rel_error)
fprintf('norm(A*x-b)=%d',norm(A*x-b))   % 1.081579e-09
x=A\b;                                  % MATLAB's solver
fprintf('\nMATLAB built-in solver: x=[%f %f %f %f %f]\n',x)
import numpy as np
from numpy.linalg import *

def jacobi(A, b, x0, tol, maxiter=200):
    """
    Performs Jacobi iterations to solve the line system of
    equations, Ax=b, starting from an initial guess, ``x0``.

    Terminates when the change in x is less than ``tol``, or
    if ``maxiter`` [default=200] iterations have been exceeded.

    Returns 3 variables:
        1.  x, the estimated solution
        2.  rel_diff, the relative difference between last 2
            iterations for x
        3.  k, the number of iterations used.  If k=maxiter,
            then the required tolerance was not met.
    """
    n = A.shape[0]
    x = x0.copy()
    x_prev = x0.copy()
    k = 0
    rel_diff = tol * 2

    while (rel_diff > tol) and (k < maxiter):

        for i in range(0, n):
            subs = 0.0
            for j in range(0, n):
                if i != j:
                    subs += A[i,j] * x_prev[j]

            x[i] = (b[i] - subs ) / A[i,i]
        k += 1

        rel_diff = norm(x - x_prev) / norm(x)
        print(x, rel_diff)
        x_prev = x.copy()

    return x, rel_diff, k

# Main code starts here
# ---------------------
GL = 1.6
d = 0.8
A = np.array([
    [1.0,      0,      0,      0,   0],
    [ GL, -(d+1),    1.0,      0,   0],
    [  0,      d, -(d+1),    1.0,   0],
    [  0,      0,      d, -(d+1), 1.0],
    [  0,      0,      0,      0, 1.0]])
b = [0.5,      0,      0,      0, 0.1]
x0 = np.zeros(5);

tol = 1E-9
maxiter = 200
x, rel_diff, k = jacobi(A, b, x0, tol, maxiter)
if k == maxiter:
    print(('WARNING: the Jacobi iterations did not '
           'converge within the required tolerance.'))
print(('The solution is %s; within a tolerance of %g, '
        'using %d iterations.' % (x, rel_diff, k)))
print('Solution error = norm(Ax-b) = %g' % \
            norm(np.dot(A,x)-b)) 
print('Condition number of A = %0.5f' % cond(A))
print('Solution from built-in functions = %s' % solve(A, b))

The relative difference \(\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 \(8.604 \times 10^{-10}\), less than the required \(1 \times 10^{-9}\).

The solution at the 58th iteration is \(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 \(A\) can be obtained by computing \(||A||||A^{-1}||\) or using MATLAB's built-in function as cond(A). The result is \(||A||||A^{-1}||=8.7259\). Note that it is the 2-norm condition number.

← Question 3 (previous step) Back to all questions Next step: Question 5 →