CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Wiki > LU decomposition

LU decomposition

From CFD-Wiki

Revision as of 12:09, 14 February 2006 by Alex Bakharev (Talk | contribs)
Jump to: navigation, search

Contents

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.


LU decomposition essentially stores the operations of Gaussian elimination in "higher-level" form (see Golub and Van Loan), so repeated solutions using the same left-hand side are computed without repetition of operations that are independent of the right-hand side.

Algorithm

Add factorization here
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.

References

My wiki