To be able to do LU decomposition in basic Sage (not numpy or scipy, say), you have to make matrices over exact rings. AA is the set of exact algebraic reals (i.e., reals that are roots to polynomial equations with integer or rational coefficients).


Now we solve for $A\mathbf{x}=\mathbf{b}$ using this decomposition. We start by solving $L\mathbf{z} = \mathbf{b}$. I.e., we find a vector $\mathbf{z}=\begin{pmatrix} z_1\\z_2\\z_3\end{pmatrix}$ so that $L\mathbf{z} = \mathbf{b}$.
[ 1 0 0] [1/2 1 0] [ 1/2 3/5 1] [ 1 0 0] [1/2 1 0] [ 1/2 3/5 1] 

Since $L$ is lower triangular we can easily find the $\mathbf{z}$ we are looking for
1 5/2 4 1 5/2 4 
Now we find $\mathbf{x}=\begin{pmatrix} x_1\\x_2\\x_3\end{pmatrix}$ so that $U\mathbf{x}=\mathbf{z}$, using the $\mathbf{z}$ we just found using $L$.
[ 2 7 3] [ 0 5/2 7/2] [ 0 0 2/5] [ 2 7 3] [ 0 5/2 7/2] [ 0 0 2/5] 
Since $U$ is upper triangular you can find $\mathbf{x}$ by back substitution.
67 15 10 67 15 10 
Let's check our solution.
(1, 2, 3) (1, 2, 3) 

A way to think about finding the matrix inverse of $A$ is that you want the matrix $B$ so that $A\mathbf{c}_1= \begin{pmatrix} 1\\0\\\vdots\\ 0\end{pmatrix}$ where $c_1$ is the first column of the of $B$. In general, $A\mathbf{c}_i=\mathbf{e}_i$ where $c_i$ is the $i$th column of $B$ and $e_i$ is the vector with zeros everywhere except in the $i$th row. If in turn you can solve for $c_1,c_2,\dots$ you can construct the matrix $B$ which will be the inverse of $A$. Let's see an example to make this believable.

Let's find the LU decomposition.

Let's construct the vector $\mathbf{e}_1$.

Let's solve for $L\mathbf{z}=\mathbf{e}_1$ as we did above.
1 1/2 1/5 1 1/2 1/5 
Let's solve for $\mathbf{c}_1$ as above.
3 1/2 1/2 3 1/2 1/2 
Notice that the three numbers I just got make up the first column of the inverse of $A$. Proceeding in this way, you'd get the second column and the third column of the inverse.
[ 3 11 16] [1/2 5/2 7/2] [ 1/2 3/2 5/2] [ 3 11 16] [1/2 5/2 7/2] [ 1/2 3/2 5/2] 
