Lab04: Systems of Equations (Contd.)

In this lab, you must implement three different methods of solving systems of equations.

Gauss elimination with back substitution
LU decomposition
Iterative Jacobi method

There is algorith and skeleton code provided for you to complete your exercise. A problem set is given at the end to make use of the programs that you wrote in the lab.

1. Gauss elimination with back substitution

Gauss elimination is a systematic way of eliminating unknowns, one at a time, from the given equations until we are left with an equation containing only one unknown. At the end of the elimination procedure, we end up with $N$ equations where the first contains all $N$ unknowns, the second contains $N-1$ unknowns, the third, $N-2$ unknowns and the last, only one unknown. Back substitution comes into play at that point. The value of the last unknown is found from the last equation with only one unknown, and then its value is plugged back into the second-last equation containing two unknonwns, and so on. At the end of back substitution, all the $N$ unknowns are found and the system of equations is solved.

Preparation step

Let the $N$ unknowns be $x_1, x_2, \ldots, x_N$ and the equations be $E_1, E_2, \ldots, E_N$.

Form a coefficient matrix $A$ from the given set of $N$ equations in $N$ unknowns so that $A_{i,j}$ element contains the coefficient of $x_j$ in $E_i$. $A$ is therefore an $N\times{N}$ matrix.
Form a matrix $B$ containing the values on the right hand sides of the $N$ equations. $B$ is a $N\times1$ vector.
Form a new $N\times(N+1)$-sized augmented matrix $C$ with the last column being the vector $B$ and the first $N\times{N}$ entries coming from $A$.

The code in the cells below reads the $A$ and $B$ matrices and generates the augmented matrix $C$.

import numpy as np
np.set_printoptions(precision=2, suppress=True)   # only for pretty printing outputs

# Read the matrices A and B from user
A = np.array([[-2, -4, -1, 1],
              [-4, -4, -4, 1],                           # Enter coefficients here for A
              [-5, 0, 4, -4],
              [3, 1, -2, 3]], dtype='float')
B = np.array([[-7, -32, -14, 17]])                          # Enter RHS values here for B
C = np.concatenate((A, B.T), axis=1)
print(C)

First Step: Gaussian elimination

The aim of this step is to remove one variable at a time from the augmented matrix $C$. Eliminating a variable means making its coefficient $0$ in the equations. Remember that in $C$, the first row is the first equation, the second is the second equation and so on. No variables are eliminated from the first equation; therefore the first row in $C$ is unaltered. We eliminate $x_1$ from the second equation, i.e. the first element in the second row of $C$ must be made $0$. We eliminate $x_1, x_2$ from the third equation; thus, the first two elements in the third row are made $0$s. Continuing the argument, we make the first $N-2$ elements as 0 in the $N-1$ row and the first $N-1$ elements as $0$ in the last row.

The result at the end of the Gaussian elimination step is a matrix that contains no $0$s in the first row, the first element as $0$ in the second row, the first two elements as $0$s in the third row, $\ldots$, the first $N-2$ elements as $0$ in the $N-1$ or the last but one row, and the first $N-1$ elements as $0$s in the last row. In fact, the last row contains $0$s for all the coefficients except for that of $x_N$ in the last but one column and some value in the last column. Mathematically, we say that all the elements $$ C_{i,j} = 0, ~~~~~~~\mbox{for}~ j < i $$ Such a matrix is said to be in Row Echelon form.

Gauss elimination is done using standard row operations on a matrix. We use the standard operation $$ R_j = kR_j - R_i $$ in our procedure. That is, we multiply all the elements in Row $j$ by a factor $k$ and then subtract the corresponding elements from Row $i$. For example, if $R_j = [ 4, 1, 2, 7]$ and $R_i = [1, 2, 3, -2]$, then $$ R_j = 3R_j - R_i $$ makes $$\begin{eqnarray*} R_j & = & [3\times4-1, 3\times1-2, 3\times2-3, 3\times7-(-2)] \\ & = & [11, 1, 3, 23] \end{eqnarray*}$$

In practice, we introduce a small computational extra into the above elimination procedure. When operating with decimal numbers on a computer, it is better to deal with large numbers. So we introduce a pivot, which is the largest element (by magnitude) in each column of $C$. We use that element to compute $k$.

Finally, here are the steps in Gaussian elimination of matrix $C$.

Repeat the following for each row r_piv beginning with the first and up to the last but one row
Find the pivot row (the row containing the pivot) in the column r_piv. Call the pivot row so found as pivot.
Swap the r_piv row with the pivot row.
After swapping, for each row row below the r_piv row, i.e. from row r_piv + 1 to the last row, do the following two steps.
1. Find the multiplicative factor k as k = C[r_piv,r_piv] / C[row,r_piv]
2. Perform the row operation C[row] = k*C[row] - C[r_piv]
We now eliminated the r_piv unknown from all the equations below the r_piv row.

Use the skeleton code below and complete the Gaussian elimination method.

#############################################################################
# Code for gauss_eliminate step. It takes the augmented matrix as input:
# named as eqns in the code. The debug parameter, if True prints a lot
# of intermediate values so that the code can be debugged.
# The function returns the augmented matrix in row echelon form so that
# it can be given to a back substitution function.
#############################################################################
def gauss_eliminate(eqns, debug=False) :
    nrows = eqns.shape[0]
    ncols = eqns.shape[1]
    
    for r_piv in range(nrows - 1) :               # Step 1
        
        # The below step finds the pivot row. Do not try to understand it
        # on the first attempt. You must read the numpy tutorial to
        # find out about np.ix_() function.         Step 2
        pivot = np.ix_(abs(eqns[r_piv:,r_piv]) == 
                       abs(eqns[r_piv:,r_piv]).max())[0][0] + r_piv
        # In rare cases, it is possible that the unknown 
        # already is not present in the input equations. 
        # If the pivot is already 0, skip remaining steps
        
        if eqns[pivot,r_piv] != 0 :
            
            # We swap rows if necessary! Python has an easy operation to swap
            # rows of a matrix. We use it. Search the Internet to understand
            # how it is done!                       Step 3
            if pivot != r_piv :
                eqns[[r_piv, pivot]] = eqns[[pivot, r_piv]]
                
                if debug :
                    print('Pivot row: ', pivot)
                    print('Swapping Rows: ', r_piv, pivot)
                    print('Swapped Equations Array: ')
                    print(eqns)
            
            # For all the rows below the pivot, 
            # row = mult_fact * row - pivot_row
            for row in range(r_piv+1, nrows) :    # Step 4
                ############################################
                # Insert your code here for Step 4A and 4B 
                ############################################
                
                if debug :
                    print('Multiplicative Factor: ', k)
                    print('After reduction of row: ', row+1)
                    print(eqns)
    # We are done with eliminating the r_piv unknown. Repeat for others.

    return eqns

Second Step: Back substitution

From the row echelon matrix, we can compute all the unknowns starting with $x_N$ from Equation $E_N$. We then substitute the value of $X_N$ in Equation $E_{N-1}$ to get the value of $x_{N-1}$. We substitute the values of $x_{N-1}$ and $x_{N-2}$ in Equation $E_{N-2}$ to get the value of $x_{N-2}$. We continue this process of working backwards all the way to the first equation $E_1$ to get $x_1$.

The code skeleton given below keeps track of the unknown whose value will be computed in the current step as var_to_compute and the current equation as eqn_num. Remember that we could have used only one variable but using two is not all that expensive (memory-wise!) and makes the code a lot clearer to understand.

Here are the steps on the row echelon matrix $R$ obtained from Gaussian elimination step.

Create a $N\times1$ vector sol to store the computed values of $x_1, x_2, \ldots, x_N$.
Start from the variable $x_N$ and Equation $E_N$ in Row $N$. Repeat the following steps for all the variables until $x_1$. Current variable to compute is sol[var_to_compute] and equation is Row eqns[var_to_compute].
Get the RHS value rhs of the current equation from the last column of the current row.
Find the sum of the products of the unknowns already computed and their coefficients in the current equation. This is $$ psum = \sum_{i={j+1}}^N R_{j,i}x_{j} $$ where $j$ is the current variable (var_to_compute) from the row eqns[eqn_num].
Adjust rhs = rhs - sum.
Compute the current variable as rhs / eqns[eqn_num,var_to_compute]

We are ready to repeat the steps to compute another unknown.

def back_substitute(eqns) :
    nrows = eqns.shape[0]
    ncols = eqns.shape[1]
    
    sol =                                           # Step 1: Insert code

    for var_to_compute in range(nrows-1, -1, -1) :  # Step 2
        eqn_num = var_to_compute
        
        # Check to see if the coefficient is 0. We cannot
        # solve the system of equation in that case.
        if eqns[eqn_num, var_to_compute] == 0 :
            print('Unique solution does not exist!')
            return False
        
        rhs =                                       # Step 3: Insert code
        
        psum = 0.0                                  # Step 4: Insert code
        # Compute psum value                        #

        rhs =                                       # Step 5: Insert code
        
        sol[var_to_compute] =                       # Step 6: Insert code

    return sol

1P. Problem Set for Gauss Elimination Method

Given the following A and B matrices, write down the actual equations that they represent and solve them.
When asked for the pincode, a mathematician gave the following system of equations. What is the pincode? $$ \begin{eqnarray*} 6x_1 + 5x_2 + 6x_3 + 4x_4 + 4x_5 + 7x_6 & = & 92 \\ 5x_1 + 8x_2 + 4x_3 + 2x_4 + 2x_5 + 4x_6 & = & 91 \\ 7x_1 + 7x_2 + x_3 + 4x_4 + 3x_5 + 2x_6 & = & 88 \\ 5x_1 + 4x_2 + 7x_3 + 3x_4 + 3x_5 + x_6 & = & 56 \\ 4x_1 + 7x_2 + 2x_3 + 4x_4 + 3x_5 + 6x_6 & = & 89 \\ 5x_1 + 6x_2 + x_3 + 8x_4 + 2x_5 + 7x_6 & = & 91 \end{eqnarray*} $$

2. LU decomposition

This method uses the information from Gaussian elimination to simplify solving systems of equations when only the $B$ matrix changes. The idea is to split the $A$ matrix (that remains the same between problems) into two parts $L$ and $U$ such that $A = LU$. $L$ and $U$ are special: $L$ is a lower triangular matrix and $U$ is an upper triangular matrix. In a lower triangular matrix, all the elements above the diagonal are 0, and in an upper triangular matrix, all the elements below the digaonal are 0. For example, if we take only the $A$ part of the row echelon form of the augmented matrix $R$ (previous section), it forms an upper triangular matrix.

We now see how to get these two matrices using the intermediate steps of Gaussian elimination. Remember, the we eliminate the unknowns one at a time. To eliminate $x_1$, we find multiplication factors that make all the elements of column 1 in rows 2 to $N$ $= 0$. Thus, there are $N-1$ multiplication factors in the first column after eliminating $x_1$. Similarly, we make $0$ all the elements of column 2 in rows 3 to $N$ by finding $N-2$ multiplication factors while eliminating $x_2$. The idea is to keep track of these multiplicative factors.

While we do the Gaussian elimination, let us keep track of the multiplication factors in each column in a matrix called MFM for multiplication factors matrix. In addition, we can make all its diagonal elements as 1.

The elimination must be done using the following formula instead of the one used in Gauss elimination - this helps compute the row-echelon matrix as the desired U.

k = C[row,r_piv] / C[r_piv,r_piv]
C[row] = C[row] - k * C[r_piv]

As the multiplication factors are used to eliminate only entries below the diagonal in the $A$ portion of the augmented matrix, the matrix containing all the multiplication factors with the diagonal elements as 1 is in a lower triangular matrix form. This is our matrix $L$.

Where do we get the $U$ matrix from? Well, we already have it! It is the row echelon matrix from the Gaussian elimination method without the last column. Remember the last column is from the matrix $B$ and has nothing to do with $A$. Thus, we have both the $L$ and $U$ matrices from one single run of the Gaussian elimination method.

You can verify that these are indeed the required $L$ and $U$ matrices by multiplying them. The product is the $A$ matrix.

The code skeleton is given in the cells below. Use the following steps to complete the code. In the code, matrix $C$ is given by the variable eqns.

Take the augmented matrix and copy the code from Gauss Eliminate until you come to the pivot step. In this simple version, we don't find any pivot nor do we swap rows. We simply assume that the diagonal elements of matrix $A$ are not 0. If they are, rearrange the matrix yourself by hand!
Create a new $N\times{N}$ matrix MFM and initialise it to an identity matrix. An identity matrix has $1$s on the diagonal and $0$s everywhere else.
Make the first column as the current column and the pivot as the diagonal element in the current column.
Compute the multiplication factors for the current column and store them in the corresponding elements of the MFM matrix.
Repeat the above step now by incrememting the current column by $1$ until the last column.

def LU_decomposition(eqns, debug=False) :
    nrows = eqns.shape[0]
    ncols = eqns.shape[1]
    temp = eqns.copy()
    
    MFM = np.eye(eqns.shape[0])            # Step 2: Done for you

    for r_piv in range(ncols - 1) :        # Repeat step
        # The diagonal element cannot be 0 in this
        # simple algorithm. Check for that.
        if temp[r_piv,r_piv] != 0 :        # Part of Step 1
            # Insert code to compute multiplication factors
            # and store in MFM matrix      # Step 4
    
    return MFM, temp[:,:-1]

Verifying that your code works!

Run your code by executing the following cells. You should get $L$ in lower triangular form, $U$ in upper triangular form and also that $A = LU$.

A = np.array([[5, 1, -4, -3], 
    [-5, -3., 1., -3],
    [-3. , 1., -2., 0],
    [5., -2., 3., -4.]
    ], dtype='float')
B = np.array([[-41., -30., -19., -3.]])

C = np.concatenate((A, B.T), axis=1)
L, U = LU_decomposition(C, debug=True)
print('############# U MATRIX ###############')
print(U)
print('############# L MATRIX ###############')
print(L)
print('############# VERIFY A = LU ###############')
print(L.dot(U))

One more problem

Verify that your code works on the following rather large system of 10 equations in 10 unknowns. Try both Gaussian elimination and LU decomposition methods. In LU decomposition, of course, you can only verify that $A = LU$.

A = np.array([[12, -1, -1, -1, -1, -1, -1, -3, -1, -1],
              [-1, 10, -1, -1, -1, -3, -1, -1, -1, -1],
              [-1, -1, 7, 0, 0, 0, -2, 1, 1, 1],
              [-2, 1, -1, 15, -5, -1, 0, -2, 3, -1],
              [1, 1, 1, 1, -22, 1, 1, 1, -4, 1],
              [-3, -3, 0, -1, 0, 9, 0, 0, 0, 0],
              [5, -3, -2, -2, 0, 3, 15, -1, -3, -1],
              [1, -2, 1, -2, 1, -2, 1, 25, -1, 1],
              [1, 1, 1, 1, 1, -1, -1, -1, 18, -1],
              [-1, -1, -1, -1, 1, 1, 1, 1, 1, 10]], dtype='float')
B = np.array([70, 0, 32, 5, 22, 0, 106, 30, 5, 69], dtype='float')