In this tutorial, we want you to:
import numpy as np
def eliminate(A, b, k):
"""Eliminate leading coefficients in the rows below the k-th row of A
in the system np.matmul(A, x) == b. The elimination is done in place."""
n = A.shape[0]
for i in range(k + 1, n):
m = A[i, k] / A[k, k]
A[i, :] = A[i, :] - m * A[k, :]
b[i] = b[i] - m * b[k]
def gauss_elimination(A, b):
"""Convert the system np.matmul(A, x) == b into the equivalent system np.matmul(U, x) == c
via Gaussian elimination. The elimination is done in place."""
if A.shape[0] == 1:
return
eliminate(A, b, 0)
gauss_elimination(A[1:, 1:], b[1:])
We want you to implement back-substitution because it follows quite naturally from Gaussian elimination (a non-singular row-reduced matrix is upper triangular and hence a candidate for back-substitution). Let's say we have the linear system $U \vec{x} = \vec{b}$ where $U$ is square and upper-triangular, the reason the method is called back-substitution is because we solve for the "back" elements of $\vec{x}$ first and work our way up the vector.
The algorithm (from the lecture slides) for back-substitution is:
$x_i = (b_i - \sum_{j = i + 1}^{n}u_{i,j}x_j)/u_{i, i}$
where $U$ is an $n \times n$ upper-triangular non-singular matrix. This one may require some tinkering on your end so we are providing you with the following function skeleton to complete:
def back_substitution(U, b):
"""Return a vector x with np.matmul(U, x) == b, where
* U is an nxn numpy matrix that is upper-triangular and non-singular
* b is an nx1 numpy vector
"""
n = U.shape[0] # get the height of matrix U (square matrix)
x = np.zeros_like(b, dtype=np.float) # x has the same shape as b, initialize it with zeros
for i in range(n-1, -1, -1): # looping through rows backwards (n-1 is the last row)
pass
# TODO: complete this function
return x
Now that we have Gaussian elimination and back-substitution implemented, complete the following function that solves an arbitrary non-singular system of linear equations (you can use any code from this tutorial):
def solve(A, b):
"""Return a vector x with np.matmul(A, x) == b, where
* A is an nxn numpy matrix that is non-singular
* b is an nx1 numpy vector
"""
# TODO: complete this function
Forward-substitution follows the same idea as back-substitution but is used
when a matrix is lower triangular. What changes would have to be made to
your back_substitution function in order for it to implement forward-substitution
instead of back-substitution? You do not need to implement these changes,
just describe what you would do.