Hilbert.java

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

/******************************************************************************
 *  Compilation:  javac -classpath .:jama.jar Hilbert.java
 *  Execution:    java -classpath .:jama.jar Hilbert N
 *  Dependencies: jama.jar
 *
 *  Compute the N-by-N Hilbert matrix and invert it. Hilbert
 *  matrix is ill-conditioned. Even though 100-by-100 Hilbert
 *  matrix is invertible mathematically, it is not numerically.
 *
 *  http://math.nist.gov/javanumerics/jama/
 *  http://math.nist.gov/javanumerics/jama/Jama-1.0.1.jar
 *
 *  % java -classpath .:jama.jar Hilbert 5
 *  1.000000  0.500000  0.333333  0.250000  0.200000
 *  0.500000  0.333333  0.250000  0.200000  0.166667
 *  0.333333  0.250000  0.200000  0.166667  0.142857
 *  0.250000  0.200000  0.166667  0.142857  0.125000
 *  0.200000  0.166667  0.142857  0.125000  0.111111
 *  condition number = 476607.2502414388
 *  error = 2.0520474208751693E-11
 *
 *  % java -classpath .:jama.jar Hilbert 10
 *  condition number = 1.6024258443197799E13
 *  error = 2.557804618845694E-4
 *
 *  % java -classpath .:jama.jar Hilbert 50
 * condition number = 5.2942010286715781E18
 * error = 358.80602399259806
 *
 *  % java -classpath .:jama.jar Hilbert 100
 *  Exception in thread "main" java.lang.RuntimeException: Matrix is singular.
 *         at Jama.LUDecomposition.solve(LUDecomposition.java:282)
 *         at Jama.Matrix.solve(Matrix.java:815)
 *         at Jama.Matrix.inverse(Matrix.java:833)
 *         at Hilbert.main(Hilbert.java:52)
 *
 ******************************************************************************/

import Jama.Matrix;

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

      // create a symmetric positive definite matrix
      double[][] a = new double[N][N];
      for (int i = 0; i < N; i++)
          for (int j = 0; j < N; j++)
              a[i][j] = 1.0 / (i + j + 1);
      Matrix A = new Matrix(a);
      Matrix B = A.inverse();
      Matrix I = Matrix.identity(N, N);

      if (N < 7) A.print(8, 6);
      StdOut.println("condition number = " + A.cond());
      StdOut.println("error = " + A.times(B).minus(I).normInf());
   }

}