In [1]:

using SymPy

In [2]:

# Create an n x n (default 3), i,j elimination matrix
function E(i,j,n=3)
    E = sympy"eye"(n)           # start with the identity
    E[i,j] = - symbols("m$i$j")  # insert the negative multiplier
    E
end

Out[2]:

E (generic function with 2 methods)

In [3]:

# Create a symbolic 3x3 A matrix
  A = [symbols("A$i$j") for i=1:3, j=1:3]

Out[3]:

\begin{bmatrix}A_{11}&A_{12}&A_{13}\\A_{21}&A_{22}&A_{23}\\A_{31}&A_{32}&A_{33}\end{bmatrix}

In [4]:

E(2,1)

Out[4]:

\begin{bmatrix}1&0&0\\- m_{21}&1&0\\0&0&1\end{bmatrix}

In [5]:

using Interact

In [7]:

@manipulate for i=1:3, j=1:3
     E(i,j)
end

Out[7]:

\begin{bmatrix}1&0&0\\0&- m_{22}&0\\0&0&1\end{bmatrix}

In [8]:

E(2,1)

Out[8]:

\begin{bmatrix}1&0&0\\- m_{21}&1&0\\0&0&1\end{bmatrix}

In [9]:

inv(E(2,1))

Out[9]:

\begin{bmatrix}1&0&0\\m_{21}&1&0\\0&0&1\end{bmatrix}

subtract m21 times the first row from the second row

In [10]:

E(2,1) * A

Out[10]:

\begin{bmatrix}A_{11}&A_{12}&A_{13}\\- A_{11} m_{21} + A_{21}&- A_{12} m_{21} + A_{22}&- A_{13} m_{21} + A_{23}\\A_{31}&A_{32}&A_{33}\end{bmatrix}

do the row operation twice

In [11]:

E(2,1)^2 * A

Out[11]:

\begin{bmatrix}A_{11}&A_{12}&A_{13}\\- 2 A_{11} m_{21} + A_{21}&- 2 A_{12} m_{21} + A_{22}&- 2 A_{13} m_{21} + A_{23}\\A_{31}&A_{32}&A_{33}\end{bmatrix}

subtract m32 times the second row from the third row

In [12]:

E(3,2) * A

Out[12]:

\begin{bmatrix}A_{11}&A_{12}&A_{13}\\A_{21}&A_{22}&A_{23}\\- A_{21} m_{32} + A_{31}&- A_{22} m_{32} + A_{32}&- A_{23} m_{32} + A_{33}\end{bmatrix}

Note that after computing E(3,2) * A, the first row is untouched so we can apply E(2,1) without any interference from row 2

In [13]:

E(2,1) * E(3,2)  * A

Out[13]:

\begin{bmatrix}A_{11}&A_{12}&A_{13}\\- A_{11} m_{21} + A_{21}&- A_{12} m_{21} + A_{22}&- A_{13} m_{21} + A_{23}\\- A_{21} m_{32} + A_{31}&- A_{22} m_{32} + A_{32}&- A_{23} m_{32} + A_{33}\end{bmatrix}

However, in the below, row 2 has changed

In [14]:

E(2,1) * A

Out[14]:

\begin{bmatrix}A_{11}&A_{12}&A_{13}\\- A_{11} m_{21} + A_{21}&- A_{12} m_{21} + A_{22}&- A_{13} m_{21} + A_{23}\\A_{31}&A_{32}&A_{33}\end{bmatrix}

so if apply E(3,2)
(meaning: subtract m32 times the second row from the third row)
it will happen with the UPDATED row 2

In [15]:

E(3,2) * E(2,1) * A

Out[15]:

\begin{bmatrix}A_{11}&A_{12}&A_{13}\\- A_{11} m_{21} + A_{21}&- A_{12} m_{21} + A_{22}&- A_{13} m_{21} + A_{23}\\A_{11} m_{21} m_{32} - A_{21} m_{32} + A_{31}&A_{12} m_{21} m_{32} - A_{22} m_{32} + A_{32}&A_{13} m_{21} m_{32} - A_{23} m_{32} + A_{33}\end{bmatrix}

The row interpretation of matrix muliply correctly acounts for m32 times the second row and m21 x m32 times the first row -- the combined effect

In [16]:

E(3,2) * E(2,1)

Out[16]:

\begin{bmatrix}1&0&0\\- m_{21}&1&0\\m_{21} m_{32}&- m_{32}&1\end{bmatrix}

In [17]:

E(2,1) * E(3,2)

Out[17]:

\begin{bmatrix}1&0&0\\- m_{21}&1&0\\0&- m_{32}&1\end{bmatrix}

Let's move to 5x5 matrices

In [18]:

E(2,1,5)

Out[18]:

\begin{bmatrix}1&0&0&0&0\\- m_{21}&1&0&0&0\\0&0&1&0&0\\0&0&0&1&0\\0&0&0&0&1\end{bmatrix}

In [19]:

E1 = E(2,1,5) * E(3,1,5) * E(4,1,5)* E(5,1,5)

Out[19]:

\begin{bmatrix}1&0&0&0&0\\- m_{21}&1&0&0&0\\- m_{31}&0&1&0&0\\- m_{41}&0&0&1&0\\- m_{51}&0&0&0&1\end{bmatrix}

In [20]:

E(3,1,5) * E(2,1,5) * E(5,1,5)* E(4,1,5)  # Why does the order not matter

Out[20]:

\begin{bmatrix}1&0&0&0&0\\- m_{21}&1&0&0&0\\- m_{31}&0&1&0&0\\- m_{41}&0&0&1&0\\- m_{51}&0&0&0&1\end{bmatrix}

In [23]:

E1 # subtracts multiples of the first row

Out[23]:

\begin{bmatrix}1&0&0&0&0\\- m_{21}&1&0&0&0\\- m_{31}&0&1&0&0\\- m_{41}&0&0&1&0\\- m_{51}&0&0&0&1\end{bmatrix}

In [24]:

inv(E1) # INVERSE means add back the same multiples of the first row

Out[24]:

\begin{bmatrix}1&0&0&0&0\\m_{21}&1&0&0&0\\m_{31}&0&1&0&0\\m_{41}&0&0&1&0\\m_{51}&0&0&0&1\end{bmatrix}

In [25]:

E(3,2,5)

Out[25]:

\begin{bmatrix}1&0&0&0&0\\0&1&0&0&0\\0&- m_{32}&1&0&0\\0&0&0&1&0\\0&0&0&0&1\end{bmatrix}

In [26]:

inv(E(3,2,5))

Out[26]:

\begin{bmatrix}1&0&0&0&0\\0&1&0&0&0\\0&m_{32}&1&0&0\\0&0&0&1&0\\0&0&0&0&1\end{bmatrix}

In [27]:

E1 * E(3,2,5)

Out[27]:

\begin{bmatrix}1&0&0&0&0\\- m_{21}&1&0&0&0\\- m_{31}&- m_{32}&1&0&0\\- m_{41}&0&0&1&0\\- m_{51}&0&0&0&1\end{bmatrix}

In [28]:

inv(E1) * inv(E(3,2,5))

Out[28]:

\begin{bmatrix}1&0&0&0&0\\m_{21}&1&0&0&0\\m_{31}&m_{32}&1&0&0\\m_{41}&0&0&1&0\\m_{51}&0&0&0&1\end{bmatrix}

In [21]:

E1 * E(3,2,5)  # Why does this have a simple looking answer?

Out[21]:

\begin{bmatrix}1&0&0&0&0\\- m_{21}&1&0&0&0\\- m_{31}&- m_{32}&1&0&0\\- m_{41}&0&0&1&0\\- m_{51}&0&0&0&1\end{bmatrix}

In [22]:

inv(E1) * inv(E(3,2,5)) # Why does this have an even slightly simpler looking answer?

Out[22]:

\begin{bmatrix}1&0&0&0&0\\m_{21}&1&0&0&0\\m_{31}&m_{32}&1&0&0\\m_{41}&0&0&1&0\\m_{51}&0&0&0&1\end{bmatrix}

In [29]:

E(3,2,5) * E1 # Why doesn't this have a simple looking answer?

Out[29]:

\begin{bmatrix}1&0&0&0&0\\- m_{21}&1&0&0&0\\m_{21} m_{32} - m_{31}&- m_{32}&1&0&0\\- m_{41}&0&0&1&0\\- m_{51}&0&0&0&1\end{bmatrix}

In [30]:

E2 = prod( E(i,2,5) for i=3:5)

Out[30]:

\begin{bmatrix}1&0&0&0&0\\0&1&0&0&0\\0&- m_{32}&1&0&0\\0&- m_{42}&0&1&0\\0&- m_{52}&0&0&1\end{bmatrix}

In [31]:

E3 = prod( E(i,3,5) for i=4:5)

Out[31]:

\begin{bmatrix}1&0&0&0&0\\0&1&0&0&0\\0&0&1&0&0\\0&0&- m_{43}&1&0\\0&0&- m_{53}&0&1\end{bmatrix}

In [32]:

E4 = E(5,4,5)

Out[32]:

\begin{bmatrix}1&0&0&0&0\\0&1&0&0&0\\0&0&1&0&0\\0&0&0&1&0\\0&0&0&- m_{54}&1\end{bmatrix}

In [33]:

E1 * E2 * E3 * E4 # Why is this simple?

Out[33]:

\begin{bmatrix}1&0&0&0&0\\- m_{21}&1&0&0&0\\- m_{31}&- m_{32}&1&0&0\\- m_{41}&- m_{42}&- m_{43}&1&0\\- m_{51}&- m_{52}&- m_{53}&- m_{54}&1\end{bmatrix}

In [34]:

L = inv(E1) * inv(E2)  * inv(E3) * inv(E4) # Why is this simple?  This is the L Matrix!!

Out[34]:

\begin{bmatrix}1&0&0&0&0\\m_{21}&1&0&0&0\\m_{31}&m_{32}&1&0&0\\m_{41}&m_{42}&m_{43}&1&0\\m_{51}&m_{52}&m_{53}&m_{54}&1\end{bmatrix}

In [35]:

E4 * E3 * E2 * E1 # Why is this a mess?  Good thing we don't need it in Gaussian Elimination

Out[35]:

\begin{bmatrix}1&0&0&0&0\\- m_{21}&1&0&0&0\\m_{21} m_{32} - m_{31}&- m_{32}&1&0&0\\- m_{21} \left(m_{32} m_{43} - m_{42}\right) + m_{31} m_{43} - m_{41}&m_{32} m_{43} - m_{42}&- m_{43}&1&0\\- m_{21} \left(- m_{32} \left(m_{43} m_{54} - m_{53}\right) + m_{42} m_{54} - m_{52}\right) - m_{31} \left(m_{43} m_{54} - m_{53}\right) + m_{41} m_{54} - m_{51}&- m_{32} \left(m_{43} m_{54} - m_{53}\right) + m_{42} m_{54} - m_{52}&m_{43} m_{54} - m_{53}&- m_{54}&1\end{bmatrix}

In [36]:

inv(L)  # Right this is the inv(L) , not how the inverse of a triangular matrix isn't so pretty in general

Out[36]:

\begin{bmatrix}1&0&0&0&0\\- m_{21}&1&0&0&0\\m_{21} m_{32} - m_{31}&- m_{32}&1&0&0\\m_{21} m_{42} - m_{41} - m_{43} \left(m_{21} m_{32} - m_{31}\right)&m_{32} m_{43} - m_{42}&- m_{43}&1&0\\m_{21} m_{52} - m_{51} - m_{53} \left(m_{21} m_{32} - m_{31}\right) - m_{54} \left(m_{21} m_{42} - m_{41} - m_{43} \left(m_{21} m_{32} - m_{31}\right)\right)&m_{32} m_{53} - m_{52} - m_{54} \left(m_{32} m_{43} - m_{42}\right)&m_{43} m_{54} - m_{53}&- m_{54}&1\end{bmatrix}

How do we solve Ly=b for the unknown y given b?¶

In [37]:

Out[37]:

\begin{bmatrix}1&0&0&0&0\\m_{21}&1&0&0&0\\m_{31}&m_{32}&1&0&0\\m_{41}&m_{42}&m_{43}&1&0\\m_{51}&m_{52}&m_{53}&m_{54}&1\end{bmatrix}

Notice that the exact answer gets messy:

In [38]:

b(i) = symbols("b$i")

Out[38]:

b (generic function with 1 method)

In [39]:

[b(i) for i=1:5]

Out[39]:

\begin{bmatrix}b_{1}\\b_{2}\\b_{3}\\b_{4}\\b_{5}\end{bmatrix}

In [40]:

y = L\[b(i) for i=1:5]

Out[40]:

\begin{bmatrix}b_{1}\\- b_{1} m_{21} + b_{2}\\- b_{1} m_{31} + b_{3} - m_{32} \left(- b_{1} m_{21} + b_{2}\right)\\- b_{1} m_{41} + b_{4} - m_{42} \left(- b_{1} m_{21} + b_{2}\right) - m_{43} \left(- b_{1} m_{31} + b_{3} - m_{32} \left(- b_{1} m_{21} + b_{2}\right)\right)\\- b_{1} m_{51} + b_{5} - m_{52} \left(- b_{1} m_{21} + b_{2}\right) - m_{53} \left(- b_{1} m_{31} + b_{3} - m_{32} \left(- b_{1} m_{21} + b_{2}\right)\right) - m_{54} \left(- b_{1} m_{41} + b_{4} - m_{42} \left(- b_{1} m_{21} + b_{2}\right) - m_{43} \left(- b_{1} m_{31} + b_{3} - m_{32} \left(- b_{1} m_{21} + b_{2}\right)\right)\right)\end{bmatrix}

Forward Substitution¶

but if you compute the answer the obvious way, it is really straightforward.

In [41]:

y[1] = b(1)

Out[41]:

$$b_{1}$$

L[2,1]*y[1] + y[2] = b(2) implies

In [42]:

y[2] = b(2) - L[2,1]*y[1]

Out[42]:

$$- b_{1} m_{21} + b_{2}$$

L[3,1]y[1] + L[3,2]y[2] + y[3] = b(3) implies

In [43]:

y[3] = b(3) - L[3,1]*y[1] - L[3,2]*y[2]

Out[43]:

$$- b_{1} m_{31} + b_{3} - m_{32} \left(- b_{1} m_{21} + b_{2}\right)$$

Putting this all together

In [44]:

y = [b(i) for i=1:5]
for i = 1:5, j=1:i-1  
    y[i] -= L[i,j]*y[j]
end

In [45]:

Out[45]:

\begin{bmatrix}b_{1}\\- b_{1} m_{21} + b_{2}\\- b_{1} m_{31} + b_{3} - m_{32} \left(- b_{1} m_{21} + b_{2}\right)\\- b_{1} m_{41} + b_{4} - m_{42} \left(- b_{1} m_{21} + b_{2}\right) - m_{43} \left(- b_{1} m_{31} + b_{3} - m_{32} \left(- b_{1} m_{21} + b_{2}\right)\right)\\- b_{1} m_{51} + b_{5} - m_{52} \left(- b_{1} m_{21} + b_{2}\right) - m_{53} \left(- b_{1} m_{31} + b_{3} - m_{32} \left(- b_{1} m_{21} + b_{2}\right)\right) - m_{54} \left(- b_{1} m_{41} + b_{4} - m_{42} \left(- b_{1} m_{21} + b_{2}\right) - m_{43} \left(- b_{1} m_{31} + b_{3} - m_{32} \left(- b_{1} m_{21} + b_{2}\right)\right)\right)\end{bmatrix}

Back Substitution¶

For upper triangular matrices with arbitrary diagonal there is an analagous algorithm known as "back substitution." Since the diagonal is not 1, there is one more divide by the diagonal in the obvious alagorithm. In the solution to Ax=b, the pivots are on the diagonal and we divide by those.

In [ ]: