# Numerical Methods

# Lecture 5: Numerical Linear Algebra I


```python
import numpy as np
import scipy.linalg as sl
```

## Learning objectives:

* Manipulation of matrices and matrix equations in Python
* Reminder on properties of matrices (from MM1): determinants, singularity etc
* Algorithms for the solution of linear systems
* Gaussian elimination, including back substitution

## Contents
## I.    Introduction - Linear (Matrix Systems)
## II.   Matrices in Python
## III.  Properties of matrices: determinants, singularity, solvability of linear systems, etc
## IV.  Gaussian Elimination - Method
## V.  Gaussian Elimination - Algorithm and Code
## VI.   Gaussian Jordan Elimination


## I.    Introduction - Linear (Matrix) Systems

Recall from your Mathematical Methods I course that the we can re-write a system of simultaneous (linear) equations in matrix form.  For example, in week 4 of MM1 you considered the following example:

\begin{eqnarray*}
  2x + 3y &=& 7 \\
   x - 4y &=& 3
\end{eqnarray*} 

and it was noted that this can be written in matrix form as 

$$
\left(
  \begin{array}{rr}
    2 & 3 \\
    1 & -4  \\
  \end{array}
\right)\left(
  \begin{array}{c}
    x \\
    y \\
  \end{array}
\right) = \left(
  \begin{array}{c}
    7 \\
    3 \\
  \end{array}
\right)
$$

We understand that this system of simultaneous equation with 2 equations always has the form of 

$$
\left(
  \begin{array}{rr}
    a & b \\
    c & d  \\
  \end{array}
\right)\left(
  \begin{array}{c}
    x \\
    y \\
  \end{array}
\right) = \left(
  \begin{array}{c}
    e \\
    f \\
  \end{array}
\right)
$$

Where $a,b,c,d,e,f$ are arbitrary constants.  


Let's call the matrix which stores the coefficients of our system of linear equations to be **A**

$$
A=
\left(
  \begin{array}{rr}
    a & b \\
    c & d  \\
  \end{array}
\right)
$$

And the matrix that contains our variables to be **x**

$$
x=
\left(
  \begin{array}{c}
    x \\
    y  \\
  \end{array}
\right)
$$

Yeah, I know it's kinda confusing why this matrix is called **x** when it's got both $x$ and $y$ inside, because instinctively you would think that it should contain only $x$. However, this is just a name we give to the matrix, and it should be well known by now that those doing Engineering and Sciences are not the best people at naming various things. We could have equally named this **c**, or **z**, or **insert some other letter of the alphabet**, as long as we remember that this is a matrix used to store the variable. To help yourself in the future, you might want to write the letters used to represent the matrix in one color, for example, red, and the letters used to represent the variables or constants in the matrix in another color for example, blue, since you should have all those wonderful colored pens from your various field trips. This would be especially helpful during your Math Methods (or whatever they may have renamed this course to be in future) you are doing this year and next year to avoid confusion, and I mean this seriously, it really helps a lot. The matrices here were intentionally not colored to encourage you to write these things down in your favorite colored pens, and reinforce which are matrices and which are not, and totally not because I am too lazy to change the color,<span class= "irony"><span style="color:red"> TOTALLY NOT! BELIEVE ME!</span></span>

And the matrix that contains the results of our system of linear equation to be **b**

$$
b=
\left(
  \begin{array}{c}
    e \\
    f  \\
  \end{array}
\right)
$$

This system of equations can be represented as the matrix equation $A\pmb{x}=\pmb{b}$




We could do the same procedure for equations which have $3$ variables inside of them, for example

\begin{eqnarray*}
  2x + 3y +z &=& 7 \\
   x - 4y + 2z&=& 3 \\
  3x + 5y + 3z&=& 10 \\ 
\end{eqnarray*} 

and it was noted that this can be written in matrix form as 

$$
\left(
  \begin{array}{rrr}
    2 & 3 & 1 \\
    1 & -4  & 2\\
    3 & 5  &3 
  \end{array}
\right)\left(
  \begin{array}{c}
    x \\
    y \\
    z \\
  \end{array}
\right) = \left(
  \begin{array}{c}
    7 \\
    3 \\
    10 \\
  \end{array}
\right)
$$

We understand that this system of simultaneous equation with 2 equations always has the form of 

$$
\left(
  \begin{array}{rr}
    a & b & c \\
    d & e & f \\
    g  &h &i  \\
  \end{array}
\right)\left(
  \begin{array}{c}
    x \\
    y \\
    z \\
  \end{array}
\right) = \left(
  \begin{array}{c}
    j \\
    k \\
    l \\
  \end{array}
\right)
$$

Where $a,b,c,d,e,f$ are arbitrary constants.  


Let's call the matrix which stores the coefficients of our system of linear equations to be **A**

$$
A=
\left(
  \begin{array}{rrr}
    a & b & c \\
    d & e & f \\
    g  &h &i  \\
  \end{array}
\right)
$$

And the matrix that contains our variables to be **x**

$$
x=
\left(
  \begin{array}{c}
    x \\
    y \\
    z \\
  \end{array}
\right)
$$


And the matrix that contains the results of our system of linear equation to be **b**

$$
b=
\left(
  \begin{array}{c}
    j \\
    k \\
    l \\
  \end{array}
\right)
$$

This system of equations can again be represented as the matrix equation $A\pmb{x}=\pmb{b}$

Note here that although the letters I have used for the matrices are in alphabetical order, it is not neccesary to be in alphabetical order. Basically, you could put whatever letter you want instead of $a, b, c, d, etc. $ that I have placed here. Indeed, because many letters in Engineering are often used to represent stuff, like $c$ is often used for the speed of light, $u,v$ is used for velocity, $i,j,k$ is often used for indexing, it is often better to not use letters as placeholder for values in the matrix and instead let's use something else to avoid confusion. While you as a human might be able to distinguish which $c$ means the speed of light, and which $c$ is simply a placeholder in your matrix, your computer might not be smart enough to do so, and might get very confused. 

For example, taking our 
$$
\left(
  \begin{array}{rrr}
    2 & 3 & 1 \\
    1 & -4  & 2\\
    3 & 5  &3 
  \end{array}
\right)\left(
  \begin{array}{c}
    x \\
    y \\
    z \\
  \end{array}
\right) = \left(
  \begin{array}{c}
    7 \\
    3 \\
    10 \\
  \end{array}
\right)
$$

and the generalization 

$$
\left(
  \begin{array}{rr}
    a & b & c \\
    d & e & f \\
    g  &h &i  \\
  \end{array}
\right)\left(
  \begin{array}{c}
    x \\
    y \\
    z \\
  \end{array}
\right) = \left(
  \begin{array}{c}
    j \\
    k \\
    l \\
  \end{array}
\right)
$$

Here, $c$ is simply a placeholder for the number 1. If you used $c$ for light speed before, the computer might get very confused and instead try to solve 

$$
\left(
  \begin{array}{rrr}
    2 & 3 & 3000000 \\
    1 & -4  & 2\\
    3 & 5  &3 
  \end{array}
\right)\left(
  \begin{array}{c}
    x \\
    y \\
    z \\
  \end{array}
\right) = \left(
  \begin{array}{c}
    7 \\
    3 \\
    10 \\
  \end{array}
\right)
$$

while the you are actually trying to solve for 

$$
\left(
  \begin{array}{rrr}
    2 & 3 & 1 \\
    1 & -4  & 2\\
    3 & 5  &3 
  \end{array}
\right)\left(
  \begin{array}{c}
    x \\
    y \\
    z \\
  \end{array}
\right) = \left(
  \begin{array}{c}
    7 \\
    3 \\
    10 \\
  \end{array}
\right)
$$
.



Indeed, to avoid your computer getting confused at why there is strange 3 million in there, 
We could instead use the position in the matrix to name them 
For example:

$$
A=
\left(
  \begin{array}{rrr}
    A_{11} & A_{12} & A_{13} \\
    A_{21} & A_{22} & A_{23} \\
    A_{31} & A_{32} & A_{33}  \\
  \end{array}
\right)
$$

and 

$$
x=
\left(
  \begin{array}{c}
    x_1 \\
    x_2 \\
    x_3 \\
  \end{array}
\right)
$$

We call them $x_1$,$x_2$,$x_3$, because we are running out of letters!

$$
b=
\left(
  \begin{array}{c}
    b_1 \\
    b_2 \\
    b_3 \\
  \end{array}
\right)
$$

Thus we write

$$
\left(
  \begin{array}{rrr}
    A_{11} & A_{12} & A_{13} \\
    A_{21} & A_{22} & A_{23} \\
    A_{31} & A_{32} & A_{33}  \\
  \end{array}
\right)
\left(
  \begin{array}{c}
    x_1 \\
    x_2 \\
    x_3 \\
  \end{array}
\right)
=
\left(
  \begin{array}{c}
    b_1 \\
    b_2 \\
    b_3 \\
  \end{array}
\right)
$$

And still, our matrix form of $Ax = b$ remains

We understand that the above can easily be extended to 4 variables and we get 

$$
A=
\left(
  \begin{array}{rrrr}
    A_{11} & A_{12} & A_{13} & A_{14} \\
    A_{21} & A_{22} & A_{23} & A_{24}\\
    A_{31} & A_{32} & A_{33} & A_{34} \\
    A_{41} & A_{42} & A_{43} & A_{44} \\
  \end{array}
\right)
$$

and 

$$
x=
\left(
  \begin{array}{c}
    x_1 \\
    x_2 \\
    x_3 \\
    x_4 \\
  \end{array}
\right)
$$

We call them $x_1$,$x_2$,$x_3$, because we are running out of letters!

$$
b=
\left(
  \begin{array}{c}
    b_1 \\
    b_2 \\
    b_3 \\
    b_4 \\
  \end{array}
\right)
$$

Thus we write

$$
\left(
  \begin{array}{rrrr}
    A_{11} & A_{12} & A_{13} & A_{14} \\
    A_{21} & A_{22} & A_{23} & A_{24}\\
    A_{31} & A_{32} & A_{33} & A_{34} \\
    A_{41} & A_{42} & A_{43} & A_{44} \\
  \end{array}
\right)
\left(
  \begin{array}{c}
    x_1 \\
    x_2 \\
    x_3 \\
    x_4 \\
  \end{array}
\right)
=
\left(
  \begin{array}{c}
    b_1 \\
    b_2 \\
    b_3 \\
    b_4 \\
  \end{array}
\right)
$$



More generally, consider the arbitrary system of $n$ linear equations for $n$ unknowns

\begin{eqnarray*}
  A_{11}x_1 + A_{12}x_2 + \dots + A_{1n}x_n &=& b_1 \\ 
  A_{21}x_1 + A_{22}x_2 + \dots + A_{2n}x_n &=& b_2 \\ 
  \vdots &=& \vdots \\ 
  A_{n1}x_1 + A_{n2}x_2 + \dots + A_{nn}x_n &=& b_n
\end{eqnarray*}

where $A_{ij}$ are the constant coefficients of the linear system, $x_j$ are the unknown variables, and $b_i$
are the terms on the right hand side (RHS).  Here the index $i$ is referring to the equation number
(the row in the matrix below), with the index $j$ referring to the component of the unknown
vector $\pmb{x}$ (the column of the matrix).

This system of equations can be represented as the matrix equation $A\pmb{x}=\pmb{b}$

$$
\left(
  \begin{array}{cccc}
    A_{11} & A_{12} & \dots & A_{1n} \\
    A_{21} & A_{22} & \dots & A_{2n} \\
    \vdots & \vdots & \ddots & \vdots \\
    A_{n1} & A_{n2} & \dots & A_{nn} \\
  \end{array}
\right)\left(
         \begin{array}{c}
           x_1 \\
           x_2 \\
           \vdots \\
           x_n \\
         \end{array}
       \right)  = \left(
                   \begin{array}{c}
                     b_1 \\
                     b_2 \\
                     \vdots \\
                     b_n \\
                   \end{array}
                 \right)
$$


We can easily solve the above $2 \times 2$ example of two equations and two unknowns using substitution (e.g. multiply the second equation by 2 and subtract the first equation from the resulting equation to eliminate $x$ and hence allowing us to find $y$, then we could compute $x$ from the first equation). We find:

$$ x=37/11, \quad y=1/11.$$

In MM1 you also considered $3 \times 3$ examples which were a little more complicated but still doable.  This lecture considers the case of $n \times n$ where $n$ could easily be billions(think about AI and ML and think about how complex human decision process can be. What are the various different factors that compelled you to make such a choice in such a situation, there could be billions of different things; even very trivial things could be important in your decision making, don't worry, you will learn how to make the computer do these equations for you, instead of doing them yourselves. The method introduced in this lecture, the Gaussian elimination is a classic, and you probably have been doing it without knowing about it. Actual real world examples would use much more advanced and sophisticated algorithms than Gaussian elimination to deal with these massive matrices, but Gaussian elimination is a good bridge between your A Level Math and becoming an AI or ML expert!)! This case arises when you solve a differential equation numerically on a discrete mesh or grid. Here you would typically obtain one unknown and one (discrete, linear or nonlinear) equation at every grid point. You could generate an arbitrarily large matrix system simply by generating a finer mesh.

Note that you will solve differential equations numerically in the follow-up course Numerical Methods II.

Cases where the matrix is non-square, i.e. of shape $m \times n$ where $m\ne n$ correspond to the (over- or under-determined) system where you have more or less equations than unknowns - we won't consider these in this lecture. 


## II.    Matrices in Python

We have already used numpy arrays to store one-dimensional vectors of numbers.

The convention is that these are generally considered to be *column* vectors and have shape $n \times 1$.

We can extend to higher dimensions through the introduction of matrices as two-dimensional arrays (more generally vectors and matrices are just two examples of tensors). 

We use subscript indices to identify each component of the array or matrix, i.e. we can identify each component of the vector $\pmb{v}$ by $v_i$, and each component of the matrix $A$ by $A_{ij}$.  

Note that it is a convention that vectors are either underlined or bold, and generally lower case letters, whereas matrices are plain capital letters.

The *dimension* or *shape* of a vector/matrix is the number of rows and columns it posesses, i.e. $n \times 1$ and $m \times n$ for the examples above.

Here is an example of how we can extend our use of the numpy array object to two dimensions in order to define a matrix $A$ and some examples of some operations we can make on it.


```python
A = np.array([[10., 2., 1.],[6., 5., 4.],[1., 4., 7.]])

print(A)
```

    [[10.  2.  1.]
     [ 6.  5.  4.]
     [ 1.  4.  7.]]



```python
# the total size of the array storing A - here 9 for a 3x3 matrix
print(np.size(A))  
```

    9



```python
# the number of dimensions of the matrix A
print(np.ndim(A) )    
```

    2



```python
# the shape of the matrix A
print(np.shape(A))
```

    (3, 3)



```python
# the transpose of the matrix A
print(A.T)
```

    [[10.  6.  1.]
     [ 2.  5.  4.]
     [ 1.  4.  7.]]



```python
# the inverse of the matrix A - computed using a scipy algorithm
print(sl.inv(A))
```


```python
# the determinant of the matrix A - computed using a scipy algorithm
print(sl.det(A))
```


```python
# Multiply A with its inverse using the @ matrix multiplication operator. 
# Note that due to roundoff errors the off diagonal values are not exactly zero.
print(A @ sl.inv(A))
```

    [[ 1.00000000e+00 -1.11022302e-16  5.55111512e-17]
     [-2.22044605e-16  1.00000000e+00  2.22044605e-16]
     [-8.32667268e-17 -3.33066907e-16  1.00000000e+00]]



```python
# The @ operator is fairly new, you may also see use of numpy.dot:
print(np.dot(A,sl.inv(A)))
```

    [[ 1.00000000e+00 -1.11022302e-16  5.55111512e-17]
     [-2.22044605e-16  1.00000000e+00  2.22044605e-16]
     [-8.32667268e-17 -3.33066907e-16  1.00000000e+00]]



```python
# same way to achieve the same thing
print(A.dot(sl.inv(A)))
```

    [[ 1.00000000e+00 -1.11022302e-16  5.55111512e-17]
     [-2.22044605e-16  1.00000000e+00  2.22044605e-16]
     [-8.32667268e-17 -3.33066907e-16  1.00000000e+00]]



```python
# note that the * operator simply does operations element-wise - here this
# is not what we want!
print(A*sl.inv(A))
```

    [[ 1.42857143 -0.15037594  0.02255639]
     [-1.71428571  2.59398496 -1.02255639]
     [ 0.14285714 -1.14285714  2.        ]]



```python
# how to initialise a vector of zeros 
print(np.zeros(3))
```

    [0. 0. 0.]



```python
# how to initialise a matrix of zeros 
print(np.zeros((3,3)))
```

    [[0. 0. 0.]
     [0. 0. 0.]
     [0. 0. 0.]]



```python
# how to initialise a 3rd-order tensor of zeros
print(np.zeros((3,3,3)))
```

    [[[0. 0. 0.]
      [0. 0. 0.]
      [0. 0. 0.]]
    
     [[0. 0. 0.]
      [0. 0. 0.]
      [0. 0. 0.]]
    
     [[0. 0. 0.]
      [0. 0. 0.]
      [0. 0. 0.]]]



```python
# how to initialise the identity matrix, I or Id
print(np.eye(3))
```

    [[1. 0. 0.]
     [0. 1. 0.]
     [0. 0. 1.]]


Let's quickly consider the $2 \times 2$ case from MM1 recreated above where we claimed that $x=37/11$ and $y=1/11$.  

To solve the matrix equation 

$$ A\pmb{x}=\pmb{b} $$

we can simply multiply both sides by the inverse of the matrix $A$ (if $A$ is invertible and if we know what the inverse is of course!):

\begin{align}
A\pmb{x} & = \pmb{b}\\
\implies A^{-1}A\pmb{x} & = A^{-1}\pmb{b}\\
\implies I\pmb{x} & = A^{-1}\pmb{b}\\
\implies \pmb{x} & = A^{-1}\pmb{b}
\end{align}

so we can find the solution $\pmb{x}$ by multiplying the inverse of $A$ with the RHS vector $\pmb{b}$.

### <span style="color:blue">Exercise 5.1: Solving a linear system </span>
Formulate and solve the linear system and check you get the answer quoted above.


```python

```

### <span style="color:blue">Aside: matrix objects </span>

Note that numpy does possess a matrix object as a sub-class of the numpy array.  We can cast the above two-dimensional arrays into matrix objects and then the star operator does yield the expected matrix product:


```python
# A is an n-dimensional array (n-2 here)
A = np.array([[10., 2., 1.],[6., 5., 4.],[1., 4., 7.]])
print(type(A))
```

    <class 'numpy.ndarray'>



```python
# this casts the array A into the matrix class
print(type(np.mat(A)))
```

    <class 'numpy.matrix'>



```python
# for these objects * is standard matrix multuiplication 
# and we can check that A*A^{-1}=I as expected
print(np.mat(A)*np.mat(sl.inv(A)))
```

    [[ 1.00000000e+00 -1.11022302e-16  5.55111512e-17]
     [-2.22044605e-16  1.00000000e+00  2.22044605e-16]
     [-8.32667268e-17 -3.33066907e-16  1.00000000e+00]]


### <span style="color:blue">Slicing </span>

Remember from last year's Python module that just as for arrays or lists, we can use *slicing*  in order to extract components of matrices, for example:


```python
A=np.array([[10., 2., 1.],[6., 5., 4.],[1., 4., 7.]])
print(A)
```

    [[10.  2.  1.]
     [ 6.  5.  4.]
     [ 1.  4.  7.]]



```python
# single entry, first row, second column
print(A[0,1])
```

    2.0



```python
# first row
print(A[0,:])
```

    [10.  2.  1.]



```python
# last row
print(A[-1,:])
```

    [1. 4. 7.]



```python
# second column
print(A[:,1])
```

    [2. 5. 4.]



```python
# extract a 2x2 sub-matrix
print(A[1:3,1:3])
```

    [[5. 4.]
     [4. 7.]]


### <span style="color:blue">Exercise 5.2: Matrix manipulation in Python </span>

Let
$$
A = \left(
  \begin{array}{ccc}
    1 & 2 & 3 \\
    4 & 5 & 6 \\
    7 & 8 & 9 \\
  \end{array}
\right)
\mathrm{\quad\quad and \quad\quad}
b = \left(
  \begin{array}{c}
    2 \\
    4 \\
    6 \\
  \end{array}
\right)
$$

- Store $A$ and $b$ in NumPy array structures. Print them.
- Print their shape and size. What is the difference ?
- Create a NumPy array $I$ containing the identity matrix $I_3$. Do $A = A+I$.
- Substitute the 3rd column of $A$ with $b$.
- Solve $Ax=b$.


```python

```

## III.    Properties of matrices: determinants, singularity, solvability of linear systems, etc

Consider $N$ linear equations in $N$ unknowns, $A\pmb{x}=\pmb{b}$. 

From MM1 you learnt that this system has a *unique solution* provided that the determinant of A, $\det(A)$, is non-zero. In this case the matrix is said to be *non-singular*.

If $\det(A)=0$ (with $A$ then termed a *singular matrix*), then the linear system does *not* have a unique solution, it may have either infinite *or* no solutions.

For example, consider

$$
\left(
  \begin{array}{rr}
    2 & 3 \\
    4 & 6  \\
  \end{array}
\right)\left(
  \begin{array}{c}
    x \\
    y \\
  \end{array}
\right) = \left(
  \begin{array}{c}
    4 \\
    8 \\
  \end{array}
\right)
$$

The second equation is simply twice the first, and hence a solution to the first equation is also automatically a solution to the second equation. 

We hence only have one *linearly-independent* equation, and our problem is under-constrained: we effectively only have one eqution for two unknowns with infinitely many possibly solutions. 

If we replaced the RHS vector with $(4,7)^T$, then the two equations would be contradictory: in this case we have no solutions.

Note that a set of vectors where one can be written as a linear sum of the others are termed *linearly-dependent*. When this is not the case the vectors are termed *linearly-independent*.

The following properties of a square $n\times n$ matrix are equivalent:


* $\det(A)\ne 0$ - A is non-singular


* The columns of $A$ are linearly independent


* The rows of $A$ are linearly independent


* The columns of A *span* $n$-dimensional space (recall MM1 - we can reach any point in $\mathbb{R}^N$ through a linear combination of these vectors - note that this is simply what the operation $A\pmb{x}$ is doing of course if you write it out)


* $A$ is invertible, i.e. there exists a matrix $A^{-1}$ such that $A^{-1}A = A A^{-1}=I$


* the matrix system $A\pmb{x}=\pmb{b}$ has a unique solution for every vector $b$


## IV.    Gaussian elimination - Method


The Gaussian elimination algorithm is simply a systematic implementation of the method of equation substitution we used above to solve the $2\times 2$ system (i.e. where we "multiply the second equation by 2 and subtract the first equation from the resulting equation to *eliminate* $x$ and hence allowing us to find $y$, we can then compute $x$ from the first equation").

So *Gaussian elimination*  as the method is atributed to the mathematician *Gauss* (although it was certainly known before his time) and *elimination* as we seek to eliminate unknowns.

To perform this method for arbitrarily large systems (on paper) we form the so-called *augmented matrix*

$$
[A|\pmb{b}] = 
\left[
  \begin{array}{rr|r}
    2 & 3 & 7 \\
    1 & -4 & 3  \\
  \end{array}
\right]
$$

Note that this encodes the equations (including the RHS values); we can perform so-called *row operations* and as long as we are consistent with what we do for the LHS and RHS components then what this system is describing will not change - a solution to the updated system will be the same as the solution to the original system. 

Our task is to update the system so we can easily read off the solution - of course this is exactly what we do via the substitution approach:

First we multiplied the second equation by 2, this yield the updated augmented matrix:

$$
\left[
  \begin{array}{rr|r}
    2 & 3 & 7 \\
    2 & -8 & 6 \\
  \end{array}
\right]
$$

We can use the following notation to describe this operation:

$$ Eq. (2) \leftarrow 2\times Eq. (2) $$

Note importantly that this does not change anything about what these pair of equations are telling us about the unknown solution vector $\pmb{x}$ which although it doesn't appear is implicilty defined by this augmented equation.

The next step was to subtract the first equation from the updated second ($ Eq. (2) \leftarrow Eq. (2) - Eq. (1) $):

$$
\left[
  \begin{array}{rr|r}
    2 & 3 & 7 \\
    0 & -11 & -1 \\
  \end{array}
\right]
$$

The square matrix that is now in the $A$ position of this augmented system is an example of an *upper-triangular* matrix - all entries below the diagonal are zero.  

For such a matrix we can perform back substitution - starting at the bottom to solve trivially for the final unknown ($y$ here which clearly takes the value $-1/-11$), and then using this knowledge working our way up to solve for each remaining unknown in turn, here just $x$ (solving $2x + 3\times (1/11) = 7$).

Note that we can perform the similar substitution if we had a lower triangular matrix, first finding the first unknown and then working our way forward through the remaining unknowns - hence in this case *forward substitution*.

Note that if we wished we could of course continue working on the augmented matrix to make the $A$ component diagonal: divide the second equation by 11 and multiply by 3 ($ Eq. (2) \leftarrow (3/11)\times Eq. (2) $) and add it to the first ($ Eq. (1) \leftarrow Eq. (1) +  Eq. (2) $):

$$
\left[
  \begin{array}{rr|r}
    2 & 0 & 7-3/11\\
    0 & -3 & -3/11 \\
  \end{array}
\right]
$$

and we can further make it the identity by dividing the rows by $2$ and $-3$ respectively ($ Eq. (1) \leftarrow (1/2)\times Eq. (1) $, $ Eq. (2) \leftarrow (-1/3)\times Eq. (2) $) :

$$
\left[
  \begin{array}{rr|r}
    1 & 0 & (7-3/11)/2 \\
    0 & 1 & 1/11 \\
  \end{array}
\right]
$$

Each of these augmented matrices encodes exactly the same information as the original matrix system in terms of the unknown vector $\pmb{x}$, and hence this is telling us that

$$ \pmb{x} = I \pmb{x} = \left[
  \begin{array}{c}
    (7-3/11)/2 \\
    1/11 \\
  \end{array}
\right]
$$

i.e. exactly the solution we found when we performed back substitution from the upper-triangular form of the augmented system.


### <span style="color:blue">Exercise 5.3: Gaussian elimination $3 \times 3$ example (by hand) </span>

Consider the system of linear equations

\begin{align*}
  2x + 3y - 4z &= 5 \\
  6x + 8y + 2z &= 3 \\
  4x + 8y - 6z &= 19
\end{align*}

write this in matrix form, form the corresponding augmented system and perform row operations until you get to upper-triangular form, find the solution using back substitution (**do this all with pen and paper**).

Write some code to check your answer using `sl.inv(A) @ b`.

You should find $x=-6$, $y=5$, $z=-1/2$.

I will also include the scheme to solve this below, but please only look at the scheme to do this if you are seriously stuck, and have given a good attempt at solving this and is still stuck. Thinking hard about something help reinforce the thing in your mind. 



```python

```


```python

```


```python

```

We begin by noting that Gaussian Elimination has two major parts for solving this equation. 

The 1st part is converting it to an upper triangle form, 

and the 2nd part is back substitution. 

We are trying to solve the system of equation given by:

\begin{align*}
  2x + 3y - 4z &= 5 \\
  6x + 8y + 2z &= 3 \\
  4x + 8y - 6z &= 19
\end{align*}



For Gaussian Elimination, the 1st part is obtaining the upper triangle form 

We begin by taking the first and second equations, so equations 1 and 2

Gaussian elimination works in sequence, so you'll see lots of "fist, second, third etc."

\begin{align*}
  2x + 3y - 4z &= 5 \\
  6x + 8y + 2z &= 3 \\
\end{align*}

We note that the first unknown we have is the unknown $x$ and so we begn with eliminating $x$

We note that the coefficient of unknown $x$ in the first equation, i.e. $2x+3y-4z=5$ is $2$

We note that the coefficient of unknown $x$ in the second equation, i.e. $6x+8y+2z=3$ is $6$

I understand that to get from $2$ of the first equation, to the $6$ of the second equation, I need to multiply $2$ by $3$ to obtain $6$. I get the $3$ from dividing the coefficient of the unknown $x$ in 2nd equation $6$ by the coefficient of the unknown $x$ in the 1st equation $2$, i.e. $6/2=3$
Thus I will multiply the 1st equation by $3$

After multiplying the first equation by 3, I will obtain the first and second equation as

\begin{align*}
  6x + 9y - 12z &= 15 \\
  6x + 8y + 2z &= 3 \\
\end{align*}

Then, I will subtract equation 1 from equation 2, so that the unknown $x$ is eliminated from equation 2. 
Subtracting equation 1 from equation 2, I obtain 

$$6x + 8y + 2z - (6x+9y-12z) = 3-15$$
$$-1y + 14z = -12$$

I will use the result of this subtraction to update my equation 2, thus, now I have an equation 2 that will not have the unknown $x$

\begin{align*}
  2x + 3y - 4z &= 5 \\
  0x - 1y + 14z &= -12 \\
  4x + 8y - 6z &= 19
\end{align*}

Thus, we have succesfully eliminated the unknown $x$ from one of the equations!

Then, we continue by taking the equations 1 and 3 

\begin{align*}
  2x + 3y - 4z &= 5 \\
  4x + 8y - 6z &= 19
\end{align*}

We note that the first unknown we have is the unknown $x$ and so we begn with eliminating $x$

We note that the coefficient of unknown $x$ in the first equation, i.e. $2x+3y-4z=5$ is $2$

We note that the coefficient of unknown $x$ in the third equation, i.e. $4x+8y-6z=19$ is $4$

I understand that to get from $2$ of the first equation, to the $4$ of the third equation, I need to multiply $2$ by $2$ to obtain $4$. I get the $2$ from dividing the coefficient of the unknown $x$ in 3rd equation $4$ by the coefficient of the unknown $x$ in the 1st equation $2$, i.e. $4/2=2$
Thus I will multiply the 1st equation by $2$

After multiplying the first equation by 2, I will obtain the first and second equation as

\begin{align*}
  4x + 6y - 8z &= 10 \\
  4x + 8y - 6z &= 19
\end{align*}

Then, I will subtract equation 1 from equation 3, so that the unknown $x$ is eliminated from equation 3. 
Subtracting equation 1 from equation 3, I obtain 

$$4x + 8y - 6z - (4x+6y-8z) = 19 - 10$$
$$2y + 2z = 9$$

I will use the result of this subtraction to update my equation 3, thus, now I have an equation 3 that will not have the unknown $x$

\begin{align*}
  2x + 3y - 4z &= 5 \\
  0x - 1y  + 14z &= -12 \\
  0x + 2y + 2z &= 9
\end{align*}

Thus, we have succesfully eliminated the unknown $x$ from another one of the equations!

Well, after we have done equation 1 equation 2, and then equation 1 equation 3, the only pair left we could use is equation 2 and equation 3. <span style="color:red">Note the sequential way we are doing this, so 1-2, then 1-3, then 2-3. If it would have 4 unknowns, and 4 equations that neither over nor underconstrains,  then you would do 1-2, 1-3, 1-4, 2-3, 2-4, 3-4. If 5 unknown, and 5 equations that neither over nor underconstrains, then 1-2, 1-3, 1-4, 1-5, 2-3, 2-4, 2-5, 3-4, 3-5, 4-5. You could extend this to even more unknowns, and more equations, assuming that the number of equations is enough to ensure than there is neither overconstraining nor underconstraning. You could do this in Python through a double for loop and editing the start and end of your loop, or use pandas library (you do not need to use pandas though)! Make sure you do the loop correctly before you begin writing the actual operations! </span>

After our previous updates, we have equation 2 and equation 3 as 

\begin{align*}
  0x - 1y  + 14z &= -12 \\
  0x + 2y + 2z &= 9
\end{align*}

We note that the first unknown we have is the unknown $y$(Note: $0x$ does not count)and so we begn with eliminating $y$

We note that the coefficient of unknown $y$ in the 2nd equation,  i.e. $-1y + 14z =-12$ is $-1$

We note that the coefficient of unknown $y$ in the 3rd equation, i.e. $2y+2z = 9$ is $2$

I understand that to get from $1$ of the 2nd equation, to the $2$ of the 3rd equation, I need to multiply $-1$ by $-2$ to obtain $2$. I get the $-2$ from dividing the coefficient of the unknown $y$ in 3rd equation $2$ by the coefficient of the unknown $y$ in the 2nd equation $-1$, i.e. $2/-1=-2$
Thus I will multiply the 1st equation by $-2$

After multiplying the first equation by $-2$, I will obtain the first and second equation as

\begin{align*}
  0x + 2y  - 28z &= 24 \\
  0x + 2y + 2z &= 9
\end{align*}

Then, I will subtract equation 2 from equation 3, so that the unknown $y$ is eliminated from equation 3. 
Subtracting equation 1 from equation 3, I obtain 

$$(2y + 2z) - (2y - 28z) = 24 - 9$$
$$ 0y - 30z = 15$$

I will use the result of this subtraction to update my equation 3, thus, now I have an equation 3 that will not have the unknown $y$

\begin{align*}
  2x + 3y - 4z &= 5 \\
  0x + y  - 14z &= 12 \\
  0x + 0y - 30z &= 15
\end{align*}

Thus, we have succesfully eliminated the unknown $y$ from one of the equations!
Now, our equations look like they have a triangle where the coefficients are $0$ and another triangle where the coefficients are $\neq 0$. Because the triangle where coefficients are $\neq 0$ is above the triangle where the coefficients are $0$, this form is called the Upper Triangle Form. 

Note that there is also the Lower Triangle Form, which is pretty much the mirror image of the Upper Triangle Form. The procedure to find the Lower Triangle Form is essentially the same as the procedure to find the Upper Triangle Form, but this time we work our way up, instead of the way down, since well it's like a mirror. So for a 3 unknown, 3 equation, instead of 1-2, 1-3, 2-3, you will go instead 3-2, 3-1, 2-1. 



The step following the obtainment of the Upper Triangle Form is Back Substitution. 

Let's start from the Upper Triangle Form

Since we are doing **back** substitution, we start from the last equation and go to the 1st equation, so for our case, equation 3, then 2 then 1. 
If you had Lower Triangle Form, well, your **back** substitution,(it's like reverse the **back** so you get, uhm, **no back**) you will start from the first equation and go to last equation, so, if you had 3 unknown, 3 equation, not overconstrain or underconstrain, then you go 1-2-3. 

Note that unlike when we were doing our Triangle Form which required 2 equations at once, back substitution only needs 1 equation at once. 

\begin{align*}
  2x + 3y - 4z &= 5 \\
  0x + y  - 14z &= 12 \\
  0x + 0y - 30z &= 15
\end{align*}

I don't know the value of any of the unknowns yet. 

From equation 3, I notice that 

$$ 0x + 0y - 30z = 15 $$

Since I don't know any unknowns yet, I cannot substitute. 

Which is basically 

$$-30z = 15$$

$$z = -15/30 = -1/2$$

So we have found one of the unknowns, and we know that $z = -1/2$

Then I will proceed to equation 2. 

I know the value of the unknown $z$ being $z=-1/2$

From equation 2, I notice that 

$$0x + y  - 14z = 12$$

I know that $z = -1/2$, so I can substitute, and I get

$$ +y - 14 \times -1/2 = 12 $$

$$ y +7 = 12 $$

$$ y = 5 $$

So we have found another one of the unknowns, and we know that $y = 5$ and $z = -1/2$

Then I will proceed to equation 1. 

I know the value of the unknown $z$ being $z=-1/2$ and the unknown $y$ being $y = 5$

From equation 1, I notice that 

$$ 2x + 3y - 4z = 5 $$

I know that $y = 5$ and $z = -1/2$, so I can substitute, and I get

$$ 2x + 3\times5 - 4 \times -1/2 = 5 $$

$$ 2x + 15 +2 =5 $$

$$ 2x = -12 $$

$$ x = -6$$

So we have found another one of the unknowns, and we know that $x = -6$ and $y = 5$ and $z = -1/2$

We can convert this into Matrix form $Ax = b$

$$
A=
\left(
  \begin{array}{rrr}
    A_{11} & A_{12} & A_{13} \\
    A_{21} & A_{22} & A_{23} \\
    A_{31} & A_{32} & A_{33}  \\
  \end{array}
\right)
$$

$$
A=
\left(
  \begin{array}{rrr}
    2 & 3 & -4 \\
    6 & 8 & 2 \\
    4 & 8 & -6  \\
  \end{array}
\right)
$$

and 

$$
x=
\left(
  \begin{array}{c}
    x_1 \\
    x_2 \\
    x_3 \\
  \end{array}
\right)
$$

$$
x=
\left(
  \begin{array}{c}
    x \\
    y \\
    z \\
  \end{array}
\right)
$$

and

$$
b=
\left(
  \begin{array}{c}
    b_1 \\
    b_2 \\
    b_3 \\
  \end{array}
\right)
$$

$$
b=
\left(
  \begin{array}{c}
    5 \\
    3 \\
    19 \\
  \end{array}
\right)
$$

Thus we write

$$
\left(
  \begin{array}{rrr}
    2 & 3 & -4 \\
    6 & 8 & 2 \\
    4 & 8 & -6  \\
  \end{array}
\right)
\left(
  \begin{array}{c}
    x \\
    y \\
    z \\
  \end{array}
\right)
=
\left(
  \begin{array}{c}
    5 \\
    3 \\
    19 \\
  \end{array}
\right)
$$

Do try to follow the same procedure outlined above to solve the equation, 

and see how would you write the above procedure in matrix form instead, 

and once you have finished, try to think about how to implement it in Python code. 

Separate your code into two parts, the part where you do the conversion to **Upper Triangle form** (or **Back Triangle Form** if you wish, noting the order of doing things), and the part where you do **Back Substitution**.  

## V.    Gaussian elimination - algorithm and code

Notice that we are free to perform the following operations on the augmented system without changing the corresponding solution:


* Exchange two rows (refer to the section on *partial pivoting* next week)


* Multiply a row by a non-zero constant ($Eq. (i)\leftarrow \lambda \times Eq.(i)$)


* Subtracting a (non-zero) multiple of one row with another ($Eq. (i)\leftarrow Eq. (i) - \lambda \times Eq.(j)$)



Note that the equation/row being subtracted is termed the *pivot*.


Let's consider the algorithm mid-way working on an arbitrary matrix system, i.e. assume that the first $k$ rows (i.e. above the horizontal dashed line in the matrix below) have already been transformed into upper-triangular form, while the equations/rows below are not yet in this form.

The augmented equation in this case can be assumed to look like

$$
\left[
  \begin{array}{rrrrrrrrr|r}
    A_{11} & A_{12} & A_{13} & \cdots & A_{1k} & \cdots & A_{1j} & \cdots & A_{1n} & b_1 \\
    0      & A_{22} & A_{23} & \cdots & A_{2k} & \cdots & A_{2j} & \cdots & A_{2n} & b_2 \\
    0      & 0      & A_{33} & \cdots & A_{3k} & \cdots & A_{3j} & \cdots & A_{3n} & b_3 \\
    \vdots & \vdots & \vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots & \vdots \\
    0      & 0      & 0      & \cdots & A_{kk} & \cdots & A_{kj} & \cdots & A_{kn} & b_k \\    
\hdashline
    \vdots & \vdots & \vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots & \vdots \\
    0      & 0      & 0      & \cdots & A_{ik} & \cdots & A_{ij} & \cdots & A_{in} & b_i \\
    \vdots & \vdots & \vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots & \vdots \\
    0      & 0      & 0      & \cdots & A_{nk} & \cdots & A_{nj} & \cdots & A_{nn} & b_n \\
\end{array}
\right]
$$

Remember that here as we are mid-way through the algorithm the $A$'s and $b$'s in the above are **not** the same as in the original system!

Our aim as a next step in the algorithm is to use row $k$ (the "pivot row") to *eliminate* $A_{ik}$, and we need to do this for all of the rows $i$ below the pivot, i.e. for all $i>k$.

Note that the zeros to the left of the leading term in the pivot row means that these operations will not mess up the fact that we have all the zeros we are looking for in the lower left part of the matrix.

To eliminate $A_{ik}$ for a single row $i$ we need to perform the operation 
$$ Eq. (i)\leftarrow Eq. (i) - \frac{A_{ik}}{A_{kk}} \times Eq.(k) $$

or equivalently

\begin{align}
A_{ij} &\leftarrow A_{ij} - \frac{A_{ik}}{A_{kk}} A_{kj}, \quad j=k,k+1,\ldots,n\\
b_i &\leftarrow b_i - \frac{A_{ik}}{A_{kk}} b_{k}
\end{align}

$j$ only needs to run from $k$ upwards as we can assume that the earlier entries in column $i$ have already been set to zero, and also that the corresponding terms from the pivot row are also zero (we don't need to perform operations that we know involve the addition of zeros!).

And to eliminate these entries for all rows below the pivot we need to repeat for all $i>k$.

### <span style="color:blue">Exercise 5.3: Gaussian elimination</span>

Write some code that takes a matrix $A$ and a vector $\pmb{b}$ and converts it into upper-triangular form using the above algorithm. For the $2 \times 2$ and $3\times 3$ examples from above compare the resulting $A$ and $\pmb{b}$ you obtain following elimination.



```python
def upper_triangle(A, b):
    """ A function to covert A into upper triangluar form through row operations.
    The same row operations are performed on the vector b.
    
    Note that this implementation does not use partial pivoting which is introduced below.
    
    Also note that A and b are overwritten, and hence we do not need to return anything
    from the function.
    """
    n = np.size(b)
    rows, cols = np.shape(A)
    # check A is square
    assert(rows == cols)
    # and check A has the same numner of rows as the size of the vector b
    assert(rows == n)

    # Loop over each pivot row - all but the last row which we will never need to use as a pivot
    for k in range(n-1):
        # Loop over each row below the pivot row, including the last row which we do need to update
        for i in range(.................
            ...........
            ...........

# Test our code on our 2x2 and 3x3 examples from above

A = .......
b = .......

upper_triangle(A, b)

# Here is a new trick for you - "pretty print"
from pprint import pprint

print('\nOur A matrix following row operations to transform it into upper-triangular form:')
pprint(A)
print('The correspondingly updated b vector:')
pprint(b)
```


```python

```

### Back substitution

Now that we have an augmented system in the upper-triangular form

$$
\left[
  \begin{array}{rrrrr|r}
    A_{11} & A_{12} & A_{13} & \cdots & A_{1n} &  b_1 \\
    0      & A_{22} & A_{23} & \cdots & A_{2n} &  b_2 \\
    0      & 0      & A_{33} & \cdots & A_{3n} &  b_3 \\
    \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\
    0      & 0      & 0      & \cdots & A_{nn} &  b_n \\    
\end{array}
\right]
$$

where the solution $\pmb{x}$ of the original system also satisfies $A\pmb{x}=\pmb{b}$ for the $A$ and $\pmb{b}$ in the above upper-triangular form (rather than the original $A$ and $\pmb{b}$!).

We can solve the final equation row to yield

$$x_n = \frac{b_n}{A_{nn}}$$

The second to last equation then yields (note we've introduced a comma in the subscripts here simply to make the more complex double indices easier to read)

\begin{align}
A_{n-1,n-1}x_{n-1} + A_{n-1,n}x_n &= b_{n-1}\\
\implies x_{n-1} = \frac{b_{n-1} - A_{n-1,n}x_n}{A_{n-1,n-1}}\\
\implies x_{n-1} = \frac{b_{n-1} - A_{n-1,n}\frac{b_n}{A_{nn}}}{A_{n-1,n-1}}
\end{align}

and so on to row $k$ which yields

\begin{align}
A_{k,k}x_{k} + A_{k,k+1}x_{k+1} +\cdots +  A_{k,n}x_n &= b_{k}\\
\iff A_{k,k}x_{k} + \sum_{j=k+1}^{n}A_{kj}x_j &= b_{k}\\
\implies x_{k} &= \left( b_k - \sum_{j=k+1}^{n}A_{kj}x_j\right)\frac{1}{A_{kk}}
\end{align}

### <span style="color:blue">Exercise 5.5: Back substitution</span>

Extend your code to perform back substitution and hence to obtain the final solution $\pmb{x}$.  Check against the solutions found earlier.  Come up with some random $n\times n$ matrices (you can use `np.random.rand` for that) and check your code against `sl.inv(A)@b` (remember to use the original $A$ and $\pmb{b}$ here of course!)


```python
# This function assumes that A is already an upper triangular matrix, 
# e.g. we have already run our upper_triangular function if needed.

def back_substitution(A, b):
    """ Function to perform back subsitution on the system Ax=b.
    
    Returns the solution x.
    
    Assumes that A is on upper triangular form.
    """
    n = np.size(b)
    # Check A is square and its number of rows and columns same as size of the vector b
    rows, cols = np.shape(A)
    assert(rows == cols)
    assert(rows == n)
    # We can/should check that A is upper triangular using np.triu which is the 
    # upper triangular part of a matrix - if A is already upper triangular, then
    # it should of course match the upper-triangular component of A!!
    assert(np.allclose(A, np.triu(A)))
    
    x = np.zeros(n)
    # start at the end (row n-1) and work backwards
    for k in range(............
        ..................

    return x


# This A is the upper triangular matrix carried forward
# from the Python box above, and b the correspondingly updated b vector.
A = ..........
b = ..........

# print the solution using our codes
x = back_substitution(A, b)
print('Our solution: ',x)  

# Reinitialise A and b !
# remember our functions overwrote them
A = np.array([[2., 3., -4.],
                 [6., 8., 2.],
                 [4., 8., -6.]])
b = np.array([5., 3., 19.])

# check our answer against what SciPy gives us by multiplying b by A inverse 
print('SciPy solution: ',sl.inv(A) @ b)

print('Success: ', np.allclose(x, sl.inv(A) @ b))
```


```python

```

## VI.    Gauss-Jordan elimination

Recall that for the augmented matrix example we did by hand above we continued past the upper-triangular form so that the augmented matrix had the identity matrix in the $A$ location. This algorithm has the name Gauss-Jordan elimination but note that it requires more operations than the conversion to upper-triangular form followed by back subsitution and so is only of academic interest.

### Matrix inversion

Note that if we were to form the augmented equation with the full identity matrix in the place of the vector $\pmb{b}$, i.e. $[A|I]$ and performed row operations exactly as above until $A$ is transformed into the identity matrix $I$, then we would be left with the inverse of $A$ in the original $I$ location, i.e.

$$ [A|I] \rightarrow [I|A^{-1}] $$

Ok, so the details above is sparse, and probably not quite sufficient. I suppose that some will get to coding right away, but let me give some hints to those who are, well, confused, and see if it helps alleviate the confusion. Those of you who don't need the hints can go straight to the next exercise. :)

I have typed a bunch of matrices, see if can find what the sequence of matrices mean. You should feel that the operations are oddly familiar. 

Note that due to my low dexterity stat, although I double checked, there may be typos in the steps below. 

$$
A=
\left(
  \begin{array}{rrr}
    2 & -1 & 0 \\
    -1 & 2 & -1 \\
    0  & -1 & 2  \\
  \end{array}
\right)
\;\;\;\;\;\;\;\;I=
\left(
  \begin{array}{rrr}
    1 & 0 & 0 \\
    0 & 1 & 0 \\
    0 & 0 & 1  \\
  \end{array}
\right)
$$

$$---------------------------------------------------------------------------------------------------------------------$$

$$
\left(
  \begin{array}{rrr}
    2 & -1 & 0 \\
    2 & -4 & +2 \\
    0  & -1 & 2  \\
  \end{array}
\right)
\;\;\;\;\;\;\;\;\;
\left(
  \begin{array}{rrr}
    1 & 0 & 0 \\
    0 & -2 & 0 \\
    0  & 0 & 1  \\
  \end{array}
\right)
$$

$$
\left(
  \begin{array}{rrr}
    2 & -1 & 0 \\
    0 & -3 & +2 \\
    0  & -1 & 2  \\
  \end{array}
\right)
\;\;\;\;\;\;\;\;\;
\left(
  \begin{array}{rrr}
    1 & 0 & 0 \\
    -1 & -2 & 0 \\
    0  & 0 & 1  \\
  \end{array}
\right)
$$

$$
\left(
  \begin{array}{rrr}
    2 & -1 & 0 \\
    0 & -3 & +2 \\
    0  & -3 & 6  \\
  \end{array}
\right)
\;\;\;\;\;\;\;\;\;
\left(
  \begin{array}{rrr}
    1 & 0 & 0 \\
    -1 & -2 & 0 \\
    0  & 0 & 3  \\
  \end{array}
\right)
$$

$$
\left(
  \begin{array}{rrr}
    2 & -1 & 0 \\
    0 & -3 & +2 \\
    0  & 0 & 4  \\
  \end{array}
\right)
\;\;\;\;\;\;\;\;\;
\left(
  \begin{array}{rrr}
    1 & 0 & 0 \\
    -1 & -2 & 0 \\
    1  & 2 & 3  \\
  \end{array}
\right)
$$

$$---------------------------------------------------------------------------------------------------------------------$$

$$
\left(
  \begin{array}{rrr}
    2 & -1 & 0 \\
    0 & -6 & +4 \\
    0  & 0 & 4  \\
  \end{array}
\right)
\;\;\;\;\;\;\;\;\;
\left(
  \begin{array}{rrr}
    1 & 0 & 0 \\
    -2 & -4 & 0 \\
    1  & 2 & 3  \\
  \end{array}
\right)
$$

$$
\left(
  \begin{array}{rrr}
    2 & -1 & 0 \\
    0 & -6 & 0 \\
    0  & 0 & 4  \\
  \end{array}
\right)
\;\;\;\;\;\;\;\;\;
\left(
  \begin{array}{rrr}
    1 & 0 & 0 \\
    -3 & -6 & -3 \\
    1  & 2 & 3  \\
  \end{array}
\right)
$$

$$
\left(
  \begin{array}{rrr}
    12 & -6 & 0 \\
    0 & -6 & 0 \\
    0  & 0 & 4  \\
  \end{array}
\right)
\;\;\;\;\;\;\;\;\;
\left(
  \begin{array}{rrr}
    6 & 0 & 0 \\
    -3 & -6 & -3 \\
    1  & 2 & 3  \\
  \end{array}
\right)
$$

$$
\left(
  \begin{array}{rrr}
    -12 & 0 & 0 \\
    0 & -6 & 0 \\
    0  & 0 & 4  \\
  \end{array}
\right)
\;\;\;\;\;\;\;\;\;
\left(
  \begin{array}{rrr}
    -9 & -6 & -3 \\
    -3 & -6 & -3 \\
    1  & 2 & 3  \\
  \end{array}
\right)
$$

$$---------------------------------------------------------------------------------------------------------------------$$

$$
\left(
  \begin{array}{rrr}
    1 & 0 & 0 \\
    0 & 1 & 0 \\
    0  & 0 & 1  \\
  \end{array}
\right)
\;\;\;\;\;\;\;\;\;
\left(
  \begin{array}{rrr}
    3/4 & 1/2 & 1/4 \\
    1/2 & 1 & 1/2 \\
    1/4 & 2/4 & 3/4  \\
  \end{array}
\right)
$$

$$---------------------------------------------------------------------------------------------------------------------$$

$$
A^{-1}=
\left(
  \begin{array}{rrr}
    3/4 & 1/2 & 1/4 \\
    1/2 & 1 & 1/2 \\
    1/4 & 2/4 & 3/4  \\
  \end{array}
\right)
$$


### <span style="color:blue">Exercise 5.6: Matrix inversion</span>

Try updating your code to construct the inverse matrix. Check your answer against the inverse matrix given by a built-in function.

Hint: Once you have performed your Gaussian elimination to transform $A$ into an upper triangular matrix, perform another elimination "from bottom to top" to transform $A$ into a diagonal matrix.


```python
# Updated version of the upper_triangular function that
# assumes that a matrix, B, is in the old vector location
# in the augmented system, and applies the same operations to
# B as to A
def upper_triangle2(A, B):
    # your code here

                
# Function which transforms the matrix into lower triangular 
# form - the point here is that if you give it a 
# matrix that is already in upper triangular form, then the 
# result will be a diagonal matrix
def lower_triangle2(A, B):
    # your code here

    
# Let's redefine A as our matrix above
A = # your code here
# and B is the identity of the corresponding size
B = # your code here

# transform A into upper triangular form (and perform the same operations on B)
# your code here

# now make this updated A lower triangular as well (the result should be diagonal)
# your code here

# final step is just to divide each row through by the value of the diagonal
# to end up with the identity in the place of A
# your code here

# print final A and B and check solution
# your code here
```


```python

```

### <span style="color:blue">Exercise 5.6: Zeros on the diagonal</span>

You may have noticed above that we have no way of guaranteeing that the $A_{kk}$ we divide through by in the Guassian elimination or back substitution algorithms is non-zero (or not very small which will also lead to computational problems).
Note also that we commented that we are free to exchange two rows in our augmented system - how could you use this fact to build robustness into our algorithms in order to deal with matrices for which our algorithms do lead to very small or zero $A_{kk}$ values?  

See if you can figure out how to do this - more on this next week!


```python

```

61e67d07d5732397649347b5f1a51162e70ea023
ipynb
Jupyter Notebook
notebooks/c_mathematics/numerical_methods/NM1_2020/NM1_2020_lecture5.ipynb
primer-computational-mathematics/book
305941b4f1fc4f15d472fd11f2c6e90741fb8b64
MIT
2020-08-02T07:32:14.000Z
2021-11-16T16:40:43.000Z
notebooks/c_mathematics/numerical_methods/NM1_2020/NM1_2020_lecture5.ipynb
primer-computational-mathematics/book
305941b4f1fc4f15d472fd11f2c6e90741fb8b64
MIT
2020-07-27T10:45:26.000Z
2020-08-12T15:09:14.000Z
notebooks/c_mathematics/numerical_methods/NM1_2020/NM1_2020_lecture5.ipynb
primer-computational-mathematics/book
305941b4f1fc4f15d472fd11f2c6e90741fb8b64
MIT
2020-08-05T13:57:32.000Z
2022-02-02T19:03:57.000Z
Qwen/Qwen-72B
1. YES
2. YES