Section 4.4 Solving a system
How to solve a system
\begin{equation*}
A\cdot v = b
\end{equation*}
after we do a LU decomposition of \(A\text{?}\) There are two steps. First, we re-write the system as
\begin{equation*}
U\cdot v = L^{-1}\cdot b.
\end{equation*}
There is no actual need to evaluate \(L^{-1}\text{,}\) it is enough to update the components of \(b\) every time we modify a row. Namely, we add to the line
A(j,k) = A(j,k) - L(j,i)*A(i,k); the line
b(j) = b(j) - L(j,i)*b(i);. This gives us automatically at the r.h.s. the vector \(c=L^{-1}\cdot b\text{.}\) System
\begin{equation*}
U\cdot v = c
\end{equation*}
is upper triangular and so, here the second point, it can be easily solved recursively. Let us go back to the \(3\times 3\) system in Subsection 4.2.2. Notice that
\begin{equation*}
\begin{pmatrix}
1&0&0\cr 4&1&0\cr 7&2&1\cr
\end{pmatrix}^{-1}
=
\begin{pmatrix}
1&0&0\cr -4&1&0\cr 1&-2&1\cr
\end{pmatrix}
\end{equation*}
and, as already pointed out,
\begin{equation*}
\begin{pmatrix}
1&0&0\cr -4&1&0\cr 1&-2&1\cr
\end{pmatrix}
\begin{pmatrix}
1\cr 1\cr 0\cr
\end{pmatrix}
=
\begin{pmatrix}
\phantom{-}1\cr -3\cr -1\cr
\end{pmatrix}.
\end{equation*}
So, after the Gaussian elimination process, the system is now
\begin{equation*}
\begin{aligned}
x+2y+3z&=1\cr
\phantom{2x+}-3y-6z&=-3\cr
\phantom{7x+y}-z&=-1\cr
\end{aligned}.
\end{equation*}
From the last we get immediately that \(z=1\text{,}\) by dividing the rhs component by the coefficient of \(z\text{.}\) Then we substitute \(z=1\) in the previous two equations, obtaining the \(2\times2\) system
\begin{equation*}
\begin{aligned}
x+2y&=-2\cr
\phantom{2x+}-3y&=3\cr
\end{aligned}.
\end{equation*}
As above, one finds immediately that \(y=-1\) and, by substituting this value in the first equation, the system reduces to the single equation
\begin{equation*}
x=0.
\end{equation*}
Hence the solution is the vector
\begin{equation*}
\begin{pmatrix}
\phantom{-}0\cr -1\cr \phantom{-}1\cr
\end{pmatrix}.
\end{equation*}
The code below implements exactly this algorithm for a general \(n\times n\) matrix. Below is the corresponding code using pivoting: