Iterative methods exercises

7.5. Iterative methods exercises#

Exercise 7.1

Using a pen and calculator, calculate the first 2 iterations of the Jacobi method for solving the system of linear equations below. Use starting values of \(x_i^{(0)} = 0 \) and work to 4 decimal places.

\[\begin{split} \begin{align*} 4x_1 +x_2 -x_3 +x_4 &=14,\\ x_1 +4x_2 -x_3 -x_4 &=10,\\ -x_1 -x_2 +5x_3 +x_4 &=-15,\\ x_1 -x_2 +x_3 +3x_4 &=3. \end{align*} \end{split}\]

Exercise 7.2

Repeat Exercise 7.1 using the Gauss-Seidel method.

Exercise 7.3

Repeat Exercise 7.1 using the SOR method using the optimum value for the relaxation parameter.

Exercise 7.4

Compute the spectral radii for the three methods used to compute the solution to the system from Exercise 7.1. Which of the three methods would you expect to converge to a solution the fastest?

Solution

The iteration matrices for the three methods applied to this system are

\[\begin{split} \begin{align*} T_J &= -D^{-1}(L + U) = \begin{pmatrix} 0 & \frac14 & -\frac14 & \frac14 \\ \frac14 & 0 & -\frac14 & -\frac14 \\ -\frac15 & -\frac15 & 0 & \frac15 \\ \frac13 & -\frac33 & \frac13 & 0 \end{pmatrix}, \\ T_{GS} &= -(L + D)^{-1}U = \begin{pmatrix} 0 & -\frac14 & \frac14 & -\frac14 \\ 0 & \frac{1}{16} & \frac{3}{16} & \frac{5}{16} \\ 0 & -\frac{3}{80} & \frac{7}{80} & -\frac{3}{16} \\ 0 & \frac{7}{60} & -\frac{1}{20} & \frac14 \end{pmatrix}, \\ T_{SOR} &= (D + \omega L)^{-1} ((1 - \omega)D - \omega U) = \begin{pmatrix} -0.1136 & -0.2784 & 0.2784 & -0.2784 \\ 0.0316 & -0.0361 & 0.2009 & 0.3559 \\ -0.0183 & -0.0700 & -0.0069 & -0.2055 \\ 0.0607 & 0.1159 & -0.0262 & 0.1981 \\ \end{pmatrix} \end{align*} \end{split}\]

The spectral radius of the iteration matrices are

\[\begin{split} \begin{align*} \rho(T_J) &= 0.6054, \\ \rho(T_{GS}) &= 0.3553, \\ \rho(T_{SOR}) &= 0.1968. \end{align*} \end{split}\]

So the SOR method should converge the fastest, followed by the Gauss-Seidel method (as expected).

Python

# Define linear system
A = np.array([[ 4, 1, -1, 1 ], 
              [1, 4, -1, -1], 
              [-1, -1, 5, 1], 
              [1, -1, 1, 3 ]])
omega = 1.1136

# Extract L, D and U matrices from A
L = np.tril(A, -1)
D = np.diag(np.diag(A))
U = np.triu(A, 1)

# Compute the iteration matrices
TJ = -np.linalg.inv(D) @ (L + U)
TGS = -np.linalg.inv(L + D) @ U
TSOR = np.linalg.inv(D + omega * L) @ ((1 - omega) * D - omega * U)

# Compute the spectral radius
rho_TJ = max(abs(np.linalg.eigvals(TJ)))
rho_TGS = max(abs(np.linalg.eigvals(TGS)))
rho_TSOR = max(abs(np.linalg.eigvals(TSOR)))

# Output results
print(f"TJ = \n", TJ)
print(f"TGS = \n", TGS)
print(f"TSOR = \n", TSOR)

print(f"\nrho_TJ = {rho_TJ:0.4f}")
print(f"rho_TGS = {rho_TGS:0.4f}")
print(f"rho_TSOR = {rho_TSOR:0.4f}")

MATLAB

% Define linear system
A = [ 4, 1, -1, 1 ;
    1, 4, -1, -1 ;
    -1, -1, 5, 1 ;
    1, -1, 1, 3 ];
omega = 1.1136;

% Extract L, D and U matrices from A
L = tril(A, -1);
D = diag(diag(A));
U = triu(A, 1);

% Compute the iteration matrices
TJ = -inv(D) * (L + U);
TGS = -inv(L + D) * U;
TSOR = inv(D + omega * L) * ((1 - omega) * D - omega * U);

% Compute the spectral radius
rho_TJ = max(abs(eig(TJ)));
rho_TGS = max(abs(eig(TGS)));
rho_TSOR = max(abs(eig(TSOR)));

% Output results
if true
    TJ
    TGS
    TSOR
    fprintf("rho_TJ = %0.4f", rho_TJ)
    fprintf("rho_TGS = %0.4f", rho_TGS)
    fprintf("rho_TSOR = %0.4f", rho_TSOR)
end

Exercise 7.5

Write a program to calculate the solution to the system of linear equations from Exercise 7.1 using the Jacobi, Gauss-Seidel and SOR methods using a convergence tolerance of \(tol=10^{-6}\). How many iterations did each of the three methods take to converge to the solution?

Solution

\(x_1 = 1.931507, x_2 = 1.821918, x_3 = -2.616438, x_4 = 1.835616\)

Jacobi: 28 iterations
Gauss-Seidel: 15 iterations
SOR: 10 iterations

Python

import numpy as np

def jacobi(A, b, tol=1e-6):
    n = len(b)
    x = np.zeros(n)
    maxiter = 100

    for k in range(maxiter):
        xold = np.copy(x)
        for i in range(n):
            sum_ = 0
            for j in range(n):
                if i != j:
                    sum_ += A[i,j] * xold[j]
        
            x[i] = (b[i] - sum_) / A[i,i]
            
        r = b - np.dot(A, x)

        if np.linalg.norm(r, np.inf) < tol:
            break

    print(f"The Jacobi method took {k} iterations to converge")
    
    return x


def gauss_seidel(A, b, tol=1e-6):
    n = len(b)
    x = np.zeros(n)
    maxiter = 100

    for k in range(maxiter):
        for i in range(n):
            sum_ = 0
            for j in range(n):
                if i != j:
                    sum_ += A[i,j] * x[j]
        
            x[i] = (b[i] - sum_) / A[i,i]
            
        r = b - np.dot(A, x)

        if np.linalg.norm(r, np.inf) < tol:
            break

    print(f"The Gauss-Seidel method took {k} iterations to converge")

    return x


def sor(A, b, omega, tol=1e-6):
    n = len(b)
    x = np.zeros(n)
    maxiter = 100

    for k in range(maxiter):
        for i in range(n):
            sum_ = 0
            for j in range(n):
                if i != j:
                    sum_ += A[i,j] * x[j]
        
            x[i] = (1 - omega) * x[i] + omega * (b[i] - sum_) / A[i,i]
            
        r = b - np.dot(A, x)

        if max(abs(r)) < tol:
            break

    print(f"The SOR method took {k} iterations to converge")

    return x


 # Define system
A = np.array([[4, 1, -1, 1], 
              [1, 4, -1, -1], 
              [-1, -1, 5, 1], 
              [1, -1, 1, 3]])
b = np.array([14, 10, -15, 3])

# Solve system
x = jacobi(A, b)
x = gauss_seidel(A, b)
x = sor(A, b, omega)

print(f"\nx1 = {x[0]:0.6f}, x2 = {x[1]:0.6f}, x3 = {x[2]:0.6f}, x4 = {x[3]:0.6f}\n")

MATLAB

function x = jacobi(A, b, tol)

n = length(b);
x = zeros(n, 1);
maxiter = 100;

for k = 0 : 100
    xold = x;
    for i = 1 : n
        sum_ = 0;
        for j = 1 : n
            if j ~= i
                sum_ = sum_ + A(i,j) * xold(j);
            end
        end
        x(i) = (b(i) - sum_) / A(i,i);
    end
    r = b - A * x;
    if norm(r, inf) < tol
        break
    end

end

fprintf("The Jacobi method took %i iterations to converge", k)

end

function x = gauss_seidel(A, b, tol)

n = length(b);
x = zeros(n, 1);
maxiter = 100;
for k = 0 : 100
    for i = 1 : n
        sum_ = 0;
        for j = 1 : n
            if j ~= i
                sum_ = sum_ + A(i,j) * x(j);
            end
        end
        x(i) = (b(i) - sum_) / A(i,i);
    end
    r = b - A * x;
    if norm(r, inf) < tol
        break
    end
end

fprintf("The Gauss-Seidel method took %i iterations to converge", k)

end

function x = sor(A, b, omega, tol)

n = length(b);
x = zeros(n, 1);
for k = 0 : 100
    for i = 1 : n
        sum_ = 0;
        for j = 1 : n
            if j ~= i
                sum_ = sum_ + A(i,j) * x(j);
            end
        end
        x(i) = (1 - omega) * x(i) + omega * (b(i) - sum_) / A(i,i);
    end
    r = b - A * x;
    if norm(r, inf) < tol
        break
    end
end

fprintf("The SOR method took %i iterations to converge", k)

end

% Define system
A = [ 4, 1, -1, 1 ;
      1, 4, -1, -1 ;
     -1, -1, 5, 1 ;
      1, -1, 1, 3 ];
b = [ 14 ; 10 ; -15 ; 3];

% Solve system
tol = 1e-6;
omega = 1.1136;
x = jacobi(A, b, tol);
x = gauss_seidel(A, b, tol);
x = sor(A, b, omega, tol)