Processing Math: 5%
To print higher-resolution math symbols, click the
Hi-res Fonts for Printing button on the jsMath control panel.
No jsMath TeX fonts found -- using image fonts instead.
These may be slow and might not print well.
Use the jsMath control panel to get additional information.
jsMath Control PanelHide this Message

jsMath
  Franz J. Vesely > CompPhys Tutorial > 2 Linear Algebra
 
 





INDEX
1 FINITE DIFFERENCES
2 LINEAR ALGEBRA
3 STOCHASTICS
4 ORDINARY DE
5 PARTIAL DE
6 SIMULATION
7 QUANTUM M.
8 HYDRODYN.
APPENDICES
 
< >
Ch. 2 Sec. 1



2.2 Exact Methods

Solve Ax=b exactly:

2.2.1 Gauss Elimination and Back Substitution


        a11 a21    a12 a22     a1N aNN                  x1    xN           =         b1    bN            

Convert this to triangular form:

        a11 0   0 a12 a22         0  aNN                  x1    xN           =         b1    bN            



Then solve the system by Back Substitution.


2.2.2 Householder Transformation

Task: Strip column and/or row vectors of some of their rear elements.

Given A – A' triangular, tridiagonal, or otherwise simple.

Let a be a vector (such as a matrix column vector), and e_{1} a unit vector. Define an auxiliary vector b \equiv a \pm |a| e_{1} and normalize it.

Householder matrix:

P \equiv 1-2 b_{0} \cdot b_{0}^{T}

See for yourself that in the transformed vector a' \equiv P \cdot a all elements but the first are equal to zero.

Apply this transformation successively to simplify a matrix.

EXAMPLE: Let

A=\left( \begin{array}{ccc}1&2&3\\ \\4&5&6 \\ \\7&8&9\end{array} \right)

From a \equiv (1,4,7)^{T} construct the auxiliary vector b=(-7.124, 4, 7)^{T}, normalize to b_{0}=(-0.662, 0.372, 0.651)^{T}. The Householder matrix is

P_{1}=\left( \begin{array}{ccc}0.123&0.492&0.862\\ \\0.492&0.724&-0.484 \\ \\0.862&-0.484&0.153\end{array} \right)

Check: Multiplying A by P_{1} we have

P_{1} \cdot A= \left( \begin{array}{ccc}8.124&2.708&11.078\\ \\0&4.602&1.464 \\ \\0&-0.696&1.062\end{array} \right)

Next step: treat the lower right 2 \times 2 submatrix of P_{1} \cdot A. From a = (4.602, -0.696)^{T} construct a 2 \times 2 Householder matrix which is promoted to a 3 \times 3 matrix by adding a trivial first line and column:

P_{2} = \left( \begin{array}{ccc}1&0&0\\ \\0&0.989&-0.149 \\ \\0&-0.149&-0.989\end{array} \right)

Result:

P_{2} \cdot P_{1} \cdot A= \left( \begin{array}{ccc}8.124&2.708&11.078\\ \\0&4.655&1.289 \\ \\0&0&-1.269\end{array} \right)



2.2.3 LU Decomposition

Task: Solve a linear system A \cdot x=b

Method (Banachiewicz, Cholesky and Crout): Factorize a matrix into two triangular matrices L(ower) and U(pper):

A= L \cdot U

with

L= \left( \begin{array}{cccc} l_{11} & 0 & . & 0 \\ l_{21} & l_{22} & . & . \\ . & & . & 0 \\ l_{N1} & . & . & l_{NN} \end{array} \right) ; \; U= \left( \begin{array}{cccc} u_{11} & u_{12} & . & u_{1N} \\ 0 & u_{22} & . & . \\ . & . & . & . \\ 0 & . & 0 & u_{NN} \end{array} \right)

When this is done, write A \cdot x=b as

L \cdot \left( U \cdot x \right) = b

such that we have the two simpler systems of equations

L \cdot y = b

and

U \cdot x = y

Since L and U are triangular these equations are easy to solve.

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


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


But first: find the matrices L and U!

Write L \cdot U= A as

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

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 L and U we may write
\begin{eqnarray} {\rm for} \;\;\; i \leq j &:& \; \sum_{k=1}^{i} l_{ik} u_{kj} = a_{ij} \; \\ {\rm for} \;\;\; i > j &:& \; \sum_{k=1}^{j} l_{ik} u_{kj} = a_{ij} \end{eqnarray}


\Longrightarrow Crout's procedure:

LU decomposition: For j=1,2, \dots N compute
\begin{eqnarray} u_{1j} &=& a_{1j} \\ u_{ij} &=& a_{ij} - \sum_{k=1}^{i-1} l_{ik} u_{kj} ; \;\;\; i=2, \dots ,j \\ l_{ij} &=& \frac{1}{u_{jj}} \left( a_{ij} - \sum_{k=1}^{j-1} l_{ik} u_{kj} \right) ; \;\;\; i=j+1, \dots ,N \end{eqnarray}



Side result: Determinant \left| A \right| = \left| L \right| \cdot \left| U \right| = u_{11} u_{22} \dots u_{NN}




EXAMPLE: Let

A=\left( \begin{array}{cc}1&2\\ \\3&4\end{array} \right)

\Longrightarrow (Crout)
\begin{eqnarray} j=1, \; i=1: \;\;\; u_{11}&=&a_{11}=1 \\ j=1, \; i=2: \;\;\; l_{21}&=&\frac{1}{u_{11}} a_{21}=3 \\ && \\ j=2, \; i=1: \;\;\; u_{12}&=&a_{12}=2 \\ j=2, \; i=2: \;\;\; u_{22}&=&a_{22}-l_{21} u_{12}=-2 \end{eqnarray}

and thus \begin{eqnarray} \left( \begin{array}{cc}1&2\\ \\3&4\end{array} \right)&=&\underbrace{\left( \begin{array}{cc}1&0\\ \\3&1\end{array} \right)} \cdot \underbrace{\left( \begin{array}{cc}1&2\\ \\0&-2\end{array} \right)} \\ && \;\;\;L \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; U \end{eqnarray}
At each step (j,i) the required elements l_{ik}, u_{kj} are already available.

Advantage of LU decomposition:

The vector b has not been touched. Therefore we may use the factors L and U of a given matrix A with different vectors b.

To find the Inverse A^{-1}: Solve A \cdot x_{j} = e_{j} with the N unit vectors e_{j}\;; combine the column vectors x_{j} to find A^{-1}.


2.2.4 Recursion Method

Find solution x if A is tri-diagonal (maybe after application of Householder).

With
\begin{eqnarray} A &\equiv& \left( \begin{array}{cccccc} \beta_{1}& \gamma_{1} & 0 & . & . & 0 \\ \alpha_{2} & \beta_{2} & \gamma_{2} & 0 & . & 0 \\ 0 & \alpha_{3} & \beta_{3} & \gamma_{3} & 0 & . \\ . & . & . & . & . & . \\ . & . & . & \alpha_{N-1} & \beta_{N-1} & \gamma_{N-1} \\ . & . & . & 0 & \alpha_{N} & \beta_{N} \end{array} \right) \end{eqnarray}

the system of equations reads
\begin{eqnarray} \beta_{1}x_{1} + \gamma_{1}x_{2}&=&b_{1} \\ \alpha_{i}x_{i-1}+\beta_{i}x_{i} + \gamma_{i}x_{i+1}&=&b_{i} ; \;\; i=2,\dots,N\!-\!1 \\ \alpha_{N}x_{N\!-\!1} + \beta_{N}x_{N}&=&b_{N} \end{eqnarray}

Introducing auxiliary variables g_{i} and h_{i} by the recursive ansatz

\begin{eqnarray} x_{i+1}&=&g_{i} x_{i} + h_{i} ; \; i=1, \dots , N\!-\!1 \end{eqnarray}

we find the "downward recursion formulae"

\textstyle \begin{eqnarray} g_{N-1}=\frac{-\alpha_{N}}{\beta_{N}} \; & ; & \; h_{N-1}=\frac{b_{N}}{\beta_{N}} \\ g_{i-1}=\frac{-\alpha_{i}}{\beta_{i}+\gamma_{i} g_{i}} \; & ; & \; h_{i-1}=\frac{b_{i}-\gamma_{i} h_{i}}{\beta_{i}+\gamma_{i} g_{i}} ; \;\; i=N-1, \dots,2 \end{eqnarray}


Having arrived at g_{1} and h_{1} we insert the known values of g_{i}, h_{i} in the "upward recursion formulae"

\begin{eqnarray} x_{1}&=&\frac{b_{1}-\gamma_{1} h_{1}}{\beta_{1}+\gamma_{1} g_{1}} \\ x_{i+1}&=&g_{i}x_{i}+h_{i} ; \;\; i=1, \dots,N\!-\!1 \end{eqnarray}


(The equation for the starting value x_{1} follows from \beta_{1} x_{1} + \gamma_{1}x_{2}=b_{1} and x_{2}=g_{1}x_{1}+h_{1}.)


EXAMPLE: In A \cdot x = b, let

A \equiv \left( \begin{array}{cccc} \beta_{1}& \gamma_{1} & 0 & 0 \\ \alpha_{2} & \beta_{2} & \gamma_{2} & 0 \\ 0 & \alpha_{3} & \beta_{3} & \gamma_{3} \\ 0 & 0 & \alpha_{4} & \beta_{4} \end{array} \right) = \left( \begin{array}{cccc} 2 & 1 & 0 & 0 \\ 2 & 3 & 1 & 0 \\ 0 & 1 & 4 & 2 \\ 0 & 0 & 1 & 3 \end{array} \right) \;\; \;\; {\rm and} \;\;\; b = \left( \begin{array}{c} 1 \\2 \\3 \\4 \end{array} \right)


Downward recursion: g_{3}=-\alpha_{4}/\beta_{4}=-1/3, h_{3}= b_{4}/\beta_{4}=4/3, and

\begin{eqnarray} i=3: \;\; g_{2}= -3/10 \; &,& \;\; h_{2}= 1/10 \\ i=2: \;\; g_{1}= -20/27 \; &,& \;\; h_{1}= 19/27 \end{eqnarray}

Upward recursion: x_{1} = 8/34, and
\begin{eqnarray} i=1: \;\; x_{2} &=& 9/17 \\ i=2: \;\; x_{3} &=& -1/17 \\ i=3: \;\; x_{4} &=& 23/17 \end{eqnarray}



vesely 2005-10-10

< >