WS04: LU Factorisation and iterative methods - partial solutions

WS04: LU Factorisation and iterative methods - partial solutions#

These exercises are indented to give you practice at using the material on numerical approximation and are intended to reinforce the material that was covered in lectures.

Please attempt the worksheet before your tutorial. Support is available in your tutorial or in the Class Team.

These are partial solutions. Please create Issues and Pull requests with your solutions.

Part a (pen and paper warm up)#

1. LU factorisation#

Use LU factorisation to solve the linear system of equations given by

\[\begin{split} \begin{pmatrix} 4 & -1 & -1 \\ 2 & 4 & 2 \\ 1 & 2 & 4 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} = \begin{pmatrix} 9 \\ -6 \\ 3 \end{pmatrix}. \end{split}\]

2. Jacobi iteration#

Rewrite the Jacobi iteration as:

\[ \vec{x}^{(k+1)} = D^{-1} \vec{b} - D^{-1}(L+U) \vec{x}^{(k)}, \]

or

\[ \vec{x}^{(k+1)} = \vec{p} - P \vec{x}^{(k)}, \]

with \(P=D^{-1}(L+U)\) and \(\vec{p}=D^{-1} \vec{b}\), which facilitates computing and coding. Apply the above Jacobi iteration to the following linear system of equations. Take three iterations starting from an initial guess \((0, 0)^T\).

\[\begin{split} \begin{pmatrix} 2 & 1 \\ 1 & 2 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} = \begin{pmatrix} 3 \\ 3 \end{pmatrix}. \end{split}\]

Note that the exact solution is \((1,1)^T\). If we want to reduce the error – \( \|\text{exact solution }- \text{approximated solution}\|\), less than 0.1, how many iterations do we need? what is the corresponding residual?

3. Gauss-Seidel iteration#

Apply Gauss-Seidel iteration to the above equations, and take three iterations and check whether it gives you more accurate approximation than Jacobi iteration – all you need to do is to update \((x_1, x_2)\) in terms of components, i.e.: update \(x_1\) first, then update \(x_2\) using your new \(x_1\).

Part b (code implementations and testing)#

4. Implement Gauss-Seidel iteration and Jacobi iteration#

Notice the only difference between Gauss-Seidel iteration and Jacobi iteration is that the former updates vector \(\vec{x}\) component by component, while latter updates vector \(\vec{x}\) as a whole.

import numpy as np

# import triangular solvers from previous workbook
%run ws03_implemented.ipynb

def Jacobi_iteration(A, b, max_iteration, x0=None):
    # we should take care to ensure that arrays are stored with the correct type - float!
    A = A.astype(np.float64)
    b = b.astype(np.float64)

    # check sizes of A and b match appropriately
    nb = len(b)
    n, m = A.shape
    if n != m:
        raise ValueError(f"A is not a square matrix! {A.shape=}")
    if n != nb:
        raise ValueError(f"shapes of A and b do not match! {A.shape=} {b.shape=}")

    # check diagonal is non zero
    for i in range(n):
        if np.isclose(A[i, i], 0):
            raise ValueError(f"A[{i}, {i}] is zero")

    # construct iteration matrices
    P = np.zeros([n, n])  # matrix P = D^{-1}(L+U)
    p = np.zeros(n)  # vector p = D^{-1} b
    for i in range(n):
        p[i] = b[i] / A[i, i]
        for j in range(n):
            P[i, j] = A[i, j] / A[i, i]
        P[i, i] = 0

    # create a new array to store the results, initialised as zero
    if x0 is None:
        x = np.zeros_like(b)
    else:
        x = x0.copy()

    # perform iteration x <- p - P * x
    for it in range(max_iteration):
        xnew = np.empty_like(x)
        for i in range(n):
            xnew[i] = p[i]
            for j in range(n):
                xnew[i] -= P[i, j] * x[j]
        x = xnew.copy()

    return x

def Gauss_Seidel_iteration(A, b, max_iteration, x0=None):
    # we should take care to ensure that arrays are stored with the correct type - float!
    A = A.astype(np.float64)
    b = b.astype(np.float64)
    if x0 is not None:
        x0 = x0.astype(np.float64)

    # check sizes of A and b match appropriately
    nb = len(b)
    n, m = A.shape
    if n != m:
        raise ValueError(f"A is not a square matrix! {A.shape=}")
    if n != nb:
        raise ValueError(f"shapes of A and b do not match! {A.shape=} {b.shape=}")

    for i in range(n):
        if np.isclose(A[i, i], 0):
            raise ValueError(f"A[{i}, {i}] is zero")

    # do not construct iteration matrices explicitly
    LD = np.zeros_like(A)
    U = np.zeros_like(A)
    for i in range(n):
        for j in range(n):
            if i < j:
                U[i, j] = A[i, j]
            else:
                LD[i, j] = A[i, j]

    # p = (L + D)^{-1} b --> found by solving triangular system
    # (L + D) p = b
    p = lower_triangular_solve(LD, b)

    # create a new array to store the results, initialised as zero
    if x0 is None:
        x = np.zeros_like(b)
    else:
        x = x0.copy()

    # perform iteration x <- p - ( P * (U * x) )
    # (L+D)(xnew - p) = -U*x
    Ux = np.empty_like(x)
    for it in range(max_iteration):
        for i in range(n):
            Ux[i] = 0.0
            for j in range(i + 1, n):
                Ux[i] += U[i, j] * x[j]

        Px = lower_triangular_solve(LD, Ux)
        x = p - Px

    return x

5. Test your iterative methods#

# Test different linear solvers starting from the above two-dimensional linear system
A = np.array([[2, 1], [1, 2]])
b = np.array([3, 3])
x_exact = np.array([1, 1])

# numpy linear solver
x0 = np.linalg.solve(A, b)
print("Solution by numpy solver:", x0)

x = Jacobi_iteration(A, b, 4)
print("Solution by Jacobi iteration: ", x)
print("Error: ", x - x_exact)
print("Residual: ", np.matmul(A, x) - b)

x = Gauss_Seidel_iteration(A, b, 4)
print("Solution by Gauss Seidel iteration: ", x)
print("Error: ", x - x_exact)
print("Residual: ", np.matmul(A, x) - b)

# Note we do not expect these values to be very small for small numbers of iterations!

A = np.array(
    [[-10, 2, 0, 67], [-2, 50, -77, 1.0e-5], [1, 7, 30, 8], [-10, -7, 0.001, 80]]
)
b = np.array([1, 2, 9, 0])

# numpy linear solvers
x0 = np.linalg.solve(A, b)
# x0 = np.linalg.inv(A).dot(b)
print("Solution by numpy solver:", x0)

x = Jacobi_iteration(A, b, 100)
print("Solution by Jacobi iteration: ", x)
print("Residual: ", np.matmul(A, x) - b)
# larger tolerance needed since the convergence is slower!
np.testing.assert_allclose(x, x0, rtol=1.0e-4)

x = Gauss_Seidel_iteration(A, b, 100)
print("Solution by Gauss Seidel iteration: ", x)
print("Residual: ", np.matmul(A, x) - b)
np.testing.assert_allclose(x, x0)

n = 20
B = np.random.rand(n, n)
eps = 10
A = eps * np.eye(n) + B * B.T
b = np.random.rand(n)

# numpy linear solvers
x0 = np.linalg.solve(A, b)
# x0 = np.linalg.inv(A).dot(b)
print("Solution by numpy solver:", x0)

x = Jacobi_iteration(A, b, 100)
print("Solution by Jacobi iteration: ", x)
print("Residual: ", np.matmul(A, x) - b)
np.testing.assert_allclose(x, x0)

x = Gauss_Seidel_iteration(A, b, 100)
print("Solution by Gauss Seidel iteration: ", x)
print("Residual: ", np.matmul(A, x) - b)
np.testing.assert_allclose(x, x0)

Part c: Extension#

6. Convergence#

For the following iteration formulas:

\[\begin{split} \begin{aligned} \vec{x}^{(k+1)} & = \begin{pmatrix} 3 \\ 1 \end{pmatrix} - \begin{pmatrix} 2 & 0 \\ 0 & 2 \end{pmatrix} \vec{x}^{(k)}, \\ \vec{x}^{(k+1)} & = \begin{pmatrix} 3 \\ 1 \end{pmatrix} - \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} \vec{x}^{(k)}, \\ \vec{x}^{(k+1)} & = \begin{pmatrix} 3 \\ 1 \end{pmatrix} - \begin{pmatrix} 1/2 & 0 \\ 0 & 1/2 \end{pmatrix} \vec{x}^{(k)}. \end{aligned} \end{split}\]

Take two or three iterations from an arbitrary initial guess. Do you think which iteration can converge?
Do you think the convergence of Jacobi iteration is completely determined by the properties of matrix \(P\)?
Wikipedia: Eigenvaluesdecomposition is the right tool to decompose a matrix to diagonal matrix, which allows us to analyse the convergence of an iteration method but is beyond the scope of this module.