LU decomposition

(Difference between revisions)
 Revision as of 01:10, 17 December 2005 (view source)Jasond (Talk | contribs) (Added material, used \vec for vectors, still needs more factorization info)← Older edit Revision as of 22:38, 18 December 2005 (view source)Jasond (Talk | contribs) m (→Description)Newer edit → Line 21: Line 21: :$:[itex] - y_i = {1 \over {L_{ii} }}\left( {b_i - \sum\limits_{j = 1}^{i-1} {L_{ij} y_j } } \right), + y_i = {1 \over {l_{ii} }}\left( {b_i - \sum\limits_{j = 1}^{i-1} {l_{ij} y_j } } \right), i = 1,2,\ldots,n$ i = 1,2,\ldots,n[/itex] Line 27: Line 27: :$:[itex] - x_i = {1 \over {U_{ii} }}\left( {y_i - \sum\limits_{j = i + 1}^n {U_{ij} x_j } } \right), + x_i = {1 \over {u_{ii} }}\left( {y_i - \sum\limits_{j = i + 1}^n {u_{ij} x_j } } \right), - i = n,n-1,\ldots,1.$ + i = n,n-1,\ldots,1.[/itex] == Algorithm == == Algorithm ==

Description

Consider the system of equations $A\vec{x}=\vec{b}$, where $A$ is an $n\times n$ nonsingular matrix. $A$ may be decomposed into an lower triangular part $L$ and an upper triangular part $U$ that will lead us to a direct procedure for the solution of the original system. This decomposition procedure is quite useful when more than one right-hand side (more than one $\vec{b}$) is to be used.

The algorithm is relatively straightforward - first, we determine the upper and lower triangular parts:

$A = LU.$

Then,

$A \vec{x} = (LU) \vec{x} = L(U\vec{x})=L\vec{y},$

where $y=Ux$. Once we solve the system

$L\vec{y} = \vec{b},$

we will be able to find the solution to the original system by solving

$U\vec{x} = \vec{b}$

The first solution is a foward substitution, while the second solution is a backward substitution. Both can be done efficiently once the factorization is available. The forward substitution may be expressed as

$y_i = {1 \over {l_{ii} }}\left( {b_i - \sum\limits_{j = 1}^{i-1} {l_{ij} y_j } } \right), i = 1,2,\ldots,n$

and the backward substitution process may be expressed as

$x_i = {1 \over {u_{ii} }}\left( {y_i - \sum\limits_{j = i + 1}^n {u_{ij} x_j } } \right), i = n,n-1,\ldots,1.$

Algorithm

Forward substitution

for k:=1 step until n do
for i:=1 step until k-1
$b_k=b_k-L_{ki}b_{i}$
end loop (i)
$b_{k}=b_{k}/L_{kk}$
end loop (k)

Backward substitution

for k:=n stepdown until 1 do
for i:=k+1 step until n
$b_k=b_k-U_{ki}b_{i}$
end loop (i)
$x_{k}=b_{k}/U_{kk}$
end loop (k)

Important Considerations

As with Gaussian elimination, LU decomposition is probably best used for relatively small, relatively non-sparse systems of equations (with small and non-sparse open to some interpretation). For larger and/or sparse problems, it would probably be best to either use an iterative method or use a direct solver package (e.g. DSCPACK) as opposed to writing one of your own.

If one has a single lefthand-side matrix and many right-hand side vectors, then LU decomposition would be a good solution procedure to consider. At the very least, it should be faster than solving each system separately with Gaussian elimination.