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

LU decomposition

From CFD-Wiki

(Difference between revisions)
Jump to: navigation, search
(Description)
m (Added reference)
Line 31: Line 31:
-
LU decomposition essentially stores the operations of Gaussian elimination in "higher-level" form, so repeated solutions using the same left-hand side are computed without repetition of operations that are independent of the right-hand side.
+
LU decomposition essentially stores the operations of Gaussian elimination in "higher-level" form (see [[#References|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 ==
== Algorithm ==
Line 54: Line 54:
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.
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 ==
 +
{{reference-book|author=Golub and Van Loan|year=1996|title=Matrix Computations|rest= The Johns Hopkins University Press, 3rd edition}}

Revision as of 22:51, 18 December 2005

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

Golub and Van Loan (1996), Matrix Computations, The Johns Hopkins University Press, 3rd edition.

My wiki