WS03: Triangular systems and Gaussian elimination#

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.

Partial solutions are available to check your understanding. Please create Issues and Pull requests with your solutions or share your answers for peer feedback in the class team Worksheet solutions channel

Part a (pen and paper warm up)#

1. Triangular systems#

Solve the upper triangular linear system given by

\[\begin{split} \begin{aligned} 2 x_1 &+& x_2 &+& 4 x_3 &=& 12 \\ && 1.5 x_2 && &=& 3 \\ && && 2 x_3 &=& 4 \end{aligned}. \end{split}\]

2. Elementary row operations#

Consider the system

\[\begin{split} \begin{aligned} x_1 + 2 x_2 & = 1 && (3) \\ 4 x_1 + x_2 & = -3 && (4). \end{aligned} \end{split}\]

Find:

  • \(2 \times (3)\) \(\rightarrow\)

  • \(0.25 \times (4)\) \(\rightarrow\)

  • \((4) + (-1) \times (3)\) \(\rightarrow\)

  • \((4) + (-4) \times (3)\) \(\rightarrow\)

3. Gaussian elimination#

Use Gaussian elimination followed by backward substitution 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}\]
\[\begin{split} \begin{pmatrix} 4 & 3 & 2 & 1 \\ 1 & 2 & 2 & 2 \\ 1 & 1 & 3 & 0 \\ 2 & 1 & 2 & 3 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \\ x_4 \end{pmatrix} = \begin{pmatrix} 10 \\ 7 \\ 5 \\ 8 \end{pmatrix}. \end{split}\]

The solution is \(\vec{x} = (1, 1, 1, 1)^T\).

Part b (code implementations and testing)#

import numpy as np

4. Implementations#

Implement the following functions with doc-strings given.

def lower_triangular_solve(A, b):
    """
    Solve the system  A x = b  where A is assumed to be lower triangular,
    i.e. A(i,j) = 0 for j > i, and the diagonal is assumed to be nonzero,
    i.e. A(i,i) != 0.

    The code checks that A is lower triangular and converts A and b to
    double precision before computing.

    ARGUMENTS:  A   lower triangular n x n array
                b   right hand side column n-vector

    RETURNS:    x   column n-vector solution
    """

    # TODO implement me

def upper_triangular_solve(A, b):
    """
    Solve the system  A x = b  where A is assumed to be upper triangular,
    i.e. A(i,j) = 0 for j < i, and the diagonal is assumed to be nonzero,
    i.e. A(i,i) != 0.

    The code checks that A is upper triangular and converts A and b to
    double precision before computing.

    ARGUMENTS:  A   upper triangular n x n array
                b   right hand side column n-vector

    RETURNS:    x   column n-vector solution
    """

    # TODO implement me

def gauss_elimination(A, b, verbose=False):
    """
    Reduce the system  A x = b  to upper triangular form, assuming that
    the diagonal is nonzero, i.e. A(i,i) != 0.

    Before computing A and b are converted to double precision.

    ARGUMENTS:  A   n x n matrix
                b   right hand side column n-vector

                verbose  (optional) if true print elimination steps

    RETURNS:    A   upper triangular n x n matrix
                b   modified column n-vector
    """

    # TODO implement me

5. Testing#

Test your solutions using your answers to part a. The first test has been implemented for you.

U = np.array([[2, 1, 4], [0, 1.5, 0], [0, 0, 2]])
b = np.array([[12], [3], [4]])

x0 = # TODO fill with precomputed answer

x = upper_triangular_solve(U, b)
np.testing.assert_almost_equal(x, x0)

# 3x3 Gaussian elimination test

# 4x4 Gaussian elimination test

Part c: Extension#

Next we want to know how good the method is. We do this by testing robustness and efficiency in practical settings

6. Robustness#

Consider the system of linear equations given by:

\[\begin{split} \begin{aligned} x_1 + x_2 & = 3 \\ \varepsilon x_1 + x_2 & = 2 + \varepsilon. \end{aligned} \end{split}\]

  • First verify that the true solution is \((x_1, x_2)^T = (1, 2)^T\).

  • Write the problem in matrix form. What does your function give for small values \(\varepsilon = 10^{-4}, 10^{-8}, 10^{-12}\)?

  • Rewrite the matrix form by considering the equations the other way round. What does your function given now?

7. Efficiency#

Consider the family of systems of linear equations given by

\[ A \vec{x} = \vec{b}, \]

where \(A\) is given by \(\varepsilon I_n + B B^T\) for a random matrix \(B\).

def generate_test_set(n):
    """
    Generate a random solvable system of linear equations.

    ARGUMENTS:  n   size of the problem

    RETURNS:    A   n x n matrix for system of linear equations
                x   n vector of solutions
                b   n vector for right hand side
    """

    B = np.random.rand(n, n)
    eps = 0.1

    A = eps * np.eye(n) + B * B.T
    x = np.ones((n, 1))
    b = np.matmul(A, x)

    return A, x, b

Test the run time of your Gaussian elimination implementation for a range of values n based on the data from generate_test_set.