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
(Added material, used \vec for vectors, still needs more factorization info)
(small typo fixed: first we solve intermediary vectory y and then x based on this y)
 
(9 intermediate revisions not shown)
Line 16: Line 16:
we will be able to find the solution to the original system by solving
we will be able to find the solution to the original system by solving
-
:<math>U\vec{x} = \vec{b}</math>
+
:<math>U\vec{x} = \vec{y}</math>
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
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
:<math>
:<math>
-
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</math>  
i = 1,2,\ldots,n</math>  
Line 27: Line 27:
:<math>
:<math>
-
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.</math>  
+
i = n,n-1,\ldots,1.</math>
 +
 
 +
 
 +
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 35: Line 38:
:    for k:=1 step until n do <br>
:    for k:=1 step until n do <br>
::    for i:=1 step until k-1 <br>
::    for i:=1 step until k-1 <br>
-
:::    <math>b_k=b_k-L_{ki}b_{i}</math> <br>
+
:::    <math>b_k=b_k-l_{ki}b_{i}</math> <br>
::    end loop (i) <br>
::    end loop (i) <br>
-
::    <math>b_{k}=b_{k}/L_{kk}</math><br>
+
::    <math>b_{k}=b_{k}/l_{kk}</math><br>
:      end loop (k)
:      end loop (k)
Backward substitution
Backward substitution
:    for k:=n stepdown until 1 do <br>
:    for k:=n stepdown until 1 do <br>
::    for i:=k+1 step until n <br>
::    for i:=k+1 step until n <br>
-
:::    <math>b_k=b_k-U_{ki}b_{i}</math> <br>
+
:::    <math>b_k=b_k-u_{ki}b_{i}</math> <br>
::    end loop (i) <br>
::    end loop (i) <br>
-
::    <math>x_{k}=b_{k}/U_{kk}</math><br>
+
::    <math>x_{k}=b_{k}/u_{kk}</math><br>
:      end loop (k)
:      end loop (k)
Line 51: 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.
 +
 +
Here a Fortran code fragment for LU decomposition written by  D. Partenov
 +
 +
do j = 1, n
 +
  do i = 1, j
 +
    if (i == 1) then
 +
      a(i,j) = a(i, j)
 +
    else
 +
      suma = 0.0
 +
      do k = 1, i-1
 +
        suma = suma + a(i,k)*a(k,j)
 +
      end do
 +
      a(i,j) = a(i,j) - suma
 +
    end if
 +
  end do
 +
  if ( j < n) then
 +
    do s = 1, n-j
 +
      i = j + s
 +
      if (j == 1) then
 +
        a(i,j) = a(i,j)/a(j,j)
 +
      else
 +
        suma = 0.0
 +
        do k = 1, j-1
 +
          suma = suma + a(i,k)*a(k,j)
 +
        end do
 +
        a(i,j) = (a(i,j) - suma)/a(j,j)
 +
      end if
 +
    end do
 +
  end if
 +
end do
 +
y(1) = b(1)
 +
do i = 2, n
 +
  suma = 0.0
 +
  do j = 1, i-1
 +
    suma = suma + a(i,j)*y(j)
 +
  end do
 +
  y(i) = b(i) - suma
 +
end do
 +
i = n
 +
j = n
 +
x(i) = y(i)/a(i,j)
 +
do s = 1, n -1
 +
  i = n - s
 +
  suma = 0.0
 +
  do j = i+1, n
 +
    suma = suma + a(i,j)*x(j)
 +
  end do
 +
  x(i) = (y(i) - suma)/a(i,i)
 +
end do
 +
 +
== References ==
 +
*{{reference-book|author=Golub and Van Loan|year=1996|title=Matrix Computations|rest= 3rd edition, The Johns Hopkins University Press, Baltimore}}
 +
*[http://en.wikipedia.org/wiki/LU_decomposition Wikipedia article ''LU decomposition]

Latest revision as of 09:22, 21 November 2011

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{y}

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.

Here a Fortran code fragment for LU decomposition written by D. Partenov 

do j = 1, n
 do i = 1, j
   if (i == 1) then
     a(i,j) = a(i, j)
   else
      suma = 0.0
     do k = 1, i-1
        suma = suma + a(i,k)*a(k,j)
     end do
     a(i,j) = a(i,j) - suma
   end if
 end do
 if ( j < n) then
   do s = 1, n-j
     i = j + s
     if (j == 1) then
       a(i,j) = a(i,j)/a(j,j)
     else
       suma = 0.0
       do k = 1, j-1
         suma = suma + a(i,k)*a(k,j)
       end do
       a(i,j) = (a(i,j) - suma)/a(j,j)
     end if
   end do
 end if
end do
y(1) = b(1)
do i = 2, n
 suma = 0.0
 do j = 1, i-1
   suma = suma + a(i,j)*y(j)
 end do
 y(i) = b(i) - suma
end do
i = n
j = n
x(i) = y(i)/a(i,j)
do s = 1, n -1
 i = n - s
 suma = 0.0
 do j = i+1, n
   suma = suma + a(i,j)*x(j)
 end do
 x(i) = (y(i) - suma)/a(i,i)
end do

References

Retrieved from "https://mail.cfd-online.com/Wiki/LU_decomposition"
My wiki