next up previous
Next: 2.2.4 Recursion Method Up: 2.2 Exact Methods Previous: 2.2.2 Householder Transformation


2.2.3 LU Decomposition

Task: Solve a linear system $\mbox{${\bf A}$} \cdot \vec{x}=\vec{b}$

Method (Banachiewicz, Cholesky and Crout): Factorize a matrix into two triangular matrices $\mbox{${\bf L}$}$(ower) and $\mbox{${\bf U}$}$(pper):

\begin{displaymath}
\mbox{${\bf A}$}=\mbox{${\bf L}$} \cdot \mbox{${\bf U}$}
\end{displaymath}

with

\begin{displaymath}
\mbox{${\bf L}$}=
\left(
\begin{array}{cccc}
l_{11} & 0 & ....
...\\
. & . & . & . \\
0 & . & 0 & u_{NN}
\end{array}\right)
\end{displaymath}

When this is done, write $\mbox{${\bf A}$} \cdot \vec{x}=\vec{b}$ as

\begin{displaymath}
\mbox{${\bf L}$} \cdot \left( \mbox{${\bf U}$} \cdot \vec{x}\right) = \vec{b}
\end{displaymath}

such that we have the two simpler systems of equations

\begin{displaymath}
\mbox{${\bf L}$} \cdot \vec{y} = \vec{b}
\end{displaymath}

and

\begin{displaymath}
\mbox{${\bf U}$} \cdot \vec{x} = \vec{y}
\end{displaymath}

Since $\mbox{${\bf L}$}$ and $\mbox{${\bf U}$}$ are triangular these equations are easy to solve.

Forward substitution:
$\displaystyle y_{1}$ $\textstyle =$ $\displaystyle \frac{1}{l_{11}} \, b_{1}$  
$\displaystyle y_{i}$ $\textstyle =$ $\displaystyle \frac{1}{l_{ii}} \, \left( b_{i} - \sum_{j=1}^{i-1} l_{ij} \,
y_{j} \right) \, ; \;\;\; i=2, \dots ,N$  



Back substitution:
$\displaystyle x_{N}$ $\textstyle =$ $\displaystyle \frac{1}{u_{NN}} \, y_{N}$  
$\displaystyle x_{i}$ $\textstyle =$ $\displaystyle \frac{1}{u_{ii}} \, \left( y_{i} - \sum_{j=i+1}^{N} u_{ij} \,
x_{j} \right) \, ; \;\;\; i=N\!-\!1, \dots ,1$  



But first: find the matrices $\mbox{${\bf L}$}$ and $\mbox{${\bf U}$}$!

Write $\mbox{${\bf L}$} \cdot \mbox{${\bf U}$}= \mbox{${\bf A}$}$ as

\begin{displaymath}
\sum_{k=1}^{N} l_{ik} \, u_{kj} = a_{ij} \, ; \; i=1, \dots N;
\; j=1, \dots N
\end{displaymath}

These are $N^{2}$ equations but $N^{2}+N$ unknowns $l_{ij}, u_{ij}$.

$\Longrightarrow$Choose $l_{ii}=1 \, (i=1, \dots N)$. Also, due to the triangular structure of $\mbox{${\bf L}$}$ and $\mbox{${\bf U}$}$ we may write
$\displaystyle {\rm for} \;\;\; i \leq j$ $\textstyle :$ $\displaystyle \; \sum_{k=1}^{i} l_{ik} \, u_{kj} = a_{ij} \;$  
$\displaystyle {\rm for} \;\;\; i > j$ $\textstyle :$ $\displaystyle \; \sum_{k=1}^{j} l_{ik} \, u_{kj} = a_{ij}$  



$\Longrightarrow$Crout's procedure:

\fbox{
\parbox{480pt}{
{\bf LU decomposition:}
For $j=1,2, \dots N$\ compute \\ ...
..._{ik} \,
u_{kj} \right) \, ; \;\;\; i=j+1, \dots ,N \nonumber
\end{eqnarray}}}



Side result: Determinant $\vert\mbox{${\bf A}$}\vert$

\begin{displaymath}
\left\vert \mbox{${\bf A}$} \right\vert = \left\vert \mbox{$...
... \mbox{${\bf U}$} \right\vert
= u_{11} \, u_{22} \dots u_{NN}
\end{displaymath}




EXAMPLE: Let

\begin{displaymath}
\mbox{${\bf A}$}=\mbox{$\left( \begin{array}{cc}1&2\\ \vspace{-9pt}\\ 3&4\end{array} \right)$}
\end{displaymath}

$\Longrightarrow$(Crout)
$\displaystyle j=1, \; i=1: \;\;\; u_{11}$ $\textstyle =$ $\displaystyle a_{11}=1$  
$\displaystyle j=1, \; i=2: \;\;\; l_{21}$ $\textstyle =$ $\displaystyle \frac{1}{u_{11}} \, a_{21}=3$  
       
$\displaystyle j=2, \; i=1: \;\;\; u_{12}$ $\textstyle =$ $\displaystyle a_{12}=2$  
$\displaystyle j=2, \; i=2: \;\;\; u_{22}$ $\textstyle =$ $\displaystyle a_{22}-l_{21} \, u_{12}=-2$  

and thus

$\displaystyle \mbox{$\left( \begin{array}{cc}1&2\\  \vspace{-9pt}\\  3&4\end{array} \right)$}$ $\textstyle =$ $\displaystyle \underbrace{\mbox{$\left( \begin{array}{cc}1&0\\  \vspace{-9pt}\\...
...begin{array}{cc}1&2\\  \vspace{-9pt}\\  0&-2\end{array} \right)$}} \nolinebreak$  
$\displaystyle \nolinebreak$   $\displaystyle \hspace{21pt} \mbox{${\bf L}$} \hspace{51pt} \mbox{${\bf U}$}$  


At each step $(j,i)$ the required elements $l_{ik}$, $u_{kj}$ are already available.

Advantage of LU decomposition:

The vector $\vec{b}$ has not been touched. Therefore we may use the factors $\mbox{${\bf L}$}$ and $\mbox{${\bf U}$}$ of a given matrix $\mbox{${\bf A}$}$ with different vectors $\vec{b}$.

To find the Inverse $\mbox{${\bf A}$}^{-1}$:

Solve $\mbox{${\bf A}$} \cdot \vec{x}_{j} = \vec{e}_{j}$, with the $N$ unit vectors $\vec{e}_{j}$; combine the column vectors $\vec{x}_{j}$ to find $\mbox{${\bf A}$}^{-1}$.
next up previous
Next: 2.2.4 Recursion Method Up: 2.2 Exact Methods Previous: 2.2.2 Householder Transformation
Franz J. Vesely Oct 2005
See also:
"Computational Physics - An Introduction," Kluwer-Plenum 2001