Gaussian Elimination

A lot of times when solving a system of equations you don't need or don't care about the inverse A matrix. You only want the vector of unknowns. In this situation, there is no need to reduce the A matrix all the way down to the identity matrix. The Gaussian elimination algorithm reduces the A matrix to an upper (or lower if you wish) diagonal matrix, as in Eq. (19.9).

Equation 19.9

graphics/19equ09.gif

The reduction of the A matrix to upper diagonal form is performed column-by-column from left to right. Whatever multiplications and additions are done to the A matrix to place it into upper diagonal form are also performed on the b vector (the primes indicate that the b vector elements have changed). Once the A matrix is reduced to upper diagonal form, the solution vector can be computed using back substitution, as shown in Eq. (19.10).

Equation 19.10

graphics/19equ10.gif

Gaussian elimination uses approximately two- thirds fewer operations than Gauss-Jordan elimination and as such is a faster algorithm. Like Gauss-Jordan, the Gaussian elimination algorithm is unstable unless the original A matrix is pivoted. Gaussian elimination can be used to determine the inverse A matrix by solving for a series of b vectors equal to each column of the identity matrix. This process requires roughly the same number of operations as determining the inverse using Gauss-Jordan elimination.

The Gaussian elimination procedure is implemented by the gaussian() method of the EqnSolver class. The method is somewhat shorter and simpler than the gaussJordan() method previously described, but the basic process is similar. The A matrix is pivoted and then reduced to an upper diagonal form. The solution vector is then solved for using back substitution.

The source code for the gaussian() method is shown here.

 public static void gaussian(double a[][], double b[]) {   int numRows = a.length;   int numCols = a[0].length;   int index[][] = new int[numRows][2];   //  Perform an implicit partial pivoting of the   //  a[][] array and b[]  array.   partialPivot(a, b, index);   //  Turn the a[][] into an upper-diagonal by   //  subtracting a multiple of the current row   //  from those below it. Do the same to the   //  b[]  array.   for(int i=0; i<numRows; ++i) {     b[i] /= a[i][i];     for(int j=numCols-1; j>=i; j) {       a[i][j] /= a[i][i];     }     for(int k=i+1; k<numRows; ++k) {       b[k] -= a[k][i]*b[i];       for(int m=i+1; m<numCols; ++m) {         a[k][m] -= a[k][i]*a[i][m];       }     }   }   //  Solve for b[]  array with back substitution   //  The diagonal elements of the a[][] matrix   //  were previously normalized to 1.   for(int i=numRows-2; i>=0; i) {     for(int j=i+1; j<numRows; ++j) {       b[i] -= a[i][j]*b[j];     }   }   return; }

A program to test the gaussian() method is presented in the "Testing the EqnSolver Class Methods" section later in this chapter.