Lower-Upper Decomposition


One disadvantage of the Gauss-Jordan and Gaussian elimination methods is that when you reduce the A matrix, you also modify the b vector. If you want to apply the original A matrix to a different b vector, either you need to keep track of every change made to the first b vector or you have to reduce the A matrix again. It is true that with Gauss-Jordan elimination you can multiply the inverse A matrix with the new b vector, but this approach is more prone to round-off error than if the new b vector is modified along with the A matrix.

The LU decomposition method is different from the Gauss-Jordan and Gaussian methods in that the reduction of the A matrix is performed independently of any b vector. This approach avoids the coupling problems we've just discussed. The reduced A matrix can be applied to any number of b vectors. It is for this reason and because it is as efficient as Gaussian elimination that LU decomposition is the most commonly used method for solving systems of equations.

The LU decomposition method is based on the fact that the A matrix can be decomposed into the product of upper and lower triangular matrices.

Equation 19.11

graphics/19equ11.gif


Eq. (19.11) converts the original system of equations into two triangular systems of equations that can be more easily and rapidly solved . The lower and upper triangular matrices take the form of Eq. (19.12).

Equation 19.12

graphics/19equ12.gif


The l ij and u ij elements are computed column by column. For a given column number j , the u ij and l ij elements are computed according to the following equations ”

Equation 19.13

graphics/19equ13.gif


Equation 19.14

graphics/19equ14.gif


The a ij terms in Eq. (19.13) and Eq. (19.14) are elements of the original A matrix. Because the l ij and u ij elements reside in different regions of the matrix, it is not necessary to allocate memory for separate L and U matrices. As we did with the implementation of Gaussian and Gauss-Jordan elimination, we will implement LU decomposition such that the A matrix is overwritten column by column with the l ij and u ij elements. This process is known as Crout's algorithm and results in the modified A matrix in Eq. (19.5).

Equation 19.15

graphics/19equ15.gif


Once the A matrix is placed into upper and lower diagonal form, the solution vector can be obtained by solving the systems of equations shown in Eq. (19.16) and Eq. (19.17).

Equation 19.16

graphics/19equ16.gif


Equation 19.17

graphics/19equ17.gif


Equation (19.16) is solved using forward substitution, Eq. (19.18).

Equation 19.18

graphics/19equ18.gif


Once the y vector is known, the solution vector can be obtained by solving Eq. (19.17) using back substitution, Eq. (19.19).

Equation 19.19

graphics/19equ19.gif


One final consideration about implementing LU decomposition is pivoting. Pivoting is essential to maintain stability of the algorithm, to avoid zeros on the diagonal of the upper diagonal matrix, for instance. Pivoting is done differently for LU decomposition than for Gaussian or Gauss-Jordan elimination in that the original A matrix is not pivoted. What is pivoted is the combined LU matrix shown in Eq. (19.15). The pivoting is performed column by column as the matrix is formed . The l ij elements are first computed without dividing them by the u jj term . The magnitude of the initial u jj diagonal element is then compared to the magnitude of the undivided l ij elements. The row containing the largest magnitude element is placed on the diagonal and is then used to divide the l ij elements.

The luDecomp() method implements the pivoted LU decomposition algorithm we have just described. Because the pivoting is done differently in this case, the pivoting code is built in to the luDecomp() method. The original A matrix is decomposed into lower and upper components with pivoting performed as each column of the LU matrix is formed. The solution vector is then computed using forward and backward substitution. The luDecomp() method source code is ”

 public static void luDecomp(double a[][], double b[]) {   double temp;   int i,j,k,m;   double tempRow[];   int numRows = a.length;   int numCols = a[0].length;   //  The a[][] matrix is overwritten into one that   //  holds the lower and upper diagonal matrices   for(j=0; j<numCols; ++j) {     //  Upper diagonal elements     for(i=0; i<=j; ++i) {       for(k=0; k<i; ++k) {         a[i][j] -= a[i][k]*a[k][j];       }     }     //  Lower diagonal elements     for(i=j+1; i<numRows; ++i) {       for(k=0; k<j; ++k) {         a[i][j] -= a[i][k]*a[k][j];       }     }     //  Determine the pivot element for the current     //  column and rearrange the rows accordingly     m=j;     for(i=j+1; i<numRows; ++i) {       if (Math.abs(a[i][j]) > Math.abs(a[m][j])){         m=i;       }     }     if (m != j) {       tempRow = a[j];       a[j] = a[m];       a[m] = tempRow;       temp = b[j];       b[j] = b[m];       b[m] = temp;     }     //  Divide lower diagonal elements by diagonal value     for(i=j+1; i<numRows; ++i) {       a[i][j] /= a[j][j];     }   }  // end of LU matrix load   //  Use forward and backward substitution to solve for   //  the unknowns. First the forward substitution.   for(i=1; i<numRows; ++i) {     for(j=0; j<i; ++j) {       b[i] -= a[i][j]*b[j];     }   }   //  And then the backward substitution.   b[numRows-1] = b[numRows-1]/a[numRows-1][numRows-1];   for(i=numRows-2; i>=0; i) {     for(j=i+1; j<numCols; ++j) {       b[i] -= a[i][j]*b[j];     }     b[i] /= a[i][i];   }   return; } 


Technical Java. Applications for Science and Engineering
Technical Java: Applications for Science and Engineering
ISBN: 0131018159
EAN: 2147483647
Year: 2003
Pages: 281
Authors: Grant Palmer

flylib.com © 2008-2017.
If you may any questions please contact us: flylib@qtcs.net