PartialPivotStability.java

Below is the syntax highlighted version of PartialPivotStability.java from §9.5 Numerical Solutions to Differential Equations.

/******************************************************************************
 *  Compilation:  javac -classpath .:jama.jar PartialPivotStability.java
 *  Execution:    java  -classpath .:jama.jar PartialPivotStability N
 *
 *  Illustrates an N-by-N example where Gaussian elimination with
 *  partial pivoting is unstable.
 *
 *  % java -classpath .:jama.jar PartialPivotStability 50
 *  condition number = 20.18837490619234
 *  inf norm of residual = 4.2537157178712537E-4
 *
 *  java -classpath .:jama.jar PartialPivotStability 100
 *  condition number = 40.53375002113499
 *  inf norm of residual = 1.62641692125309E8
 *
 *  % java -classpath .:jama.jar PartialPivotStability 200
 *  condition number = 81.24512055332623
 *  inf norm of residual = 3.632547887822304E35
 *
 ******************************************************************************/

import Jama.Matrix;

public class PartialPivotStability {
   public static void main(String[] args) {
      int N = Integer.parseInt(args[0]);

      // entries for N-by-N matrix A
      double[][] array = new double[N][N];
      for (int i = 0; i < N; i++)
          for (int j = 0; j < i; j++)
              array[i][j] = -0.9;
      for (int i = 0; i < N; i++)
          array[i][i] = 1.1;
      for (int i = 0; i < N; i++)
          array[i][N-1] = 1;

      // entries for N-by-1 vector b
      double[][] barray = new double[N][1];
      barray[0][0] = 1;

      Matrix A = new Matrix(array);
      Matrix b = new Matrix(barray);
      Matrix x = A.solve(b);
      Matrix residual = A.times(x).minus(b);
      double rnorm = residual.normInf();
      StdOut.println("condition number = " + A.cond());
      StdOut.println("inf norm of residual = " + rnorm);
   }

}