R : Copyright 2005, The R Foundation for Statistical Computing
Version 2.1.1  (2005-06-20), ISBN 3-900051-07-0

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for a HTML browser interface to help.
Type 'q()' to quit R.

>  # chex35.r                         J F Monahan, August 2007
>  # Example 3.5 Condition of Hilbert matrix
>  #
>  # set n = size of matrix
>   n <- 5                          # a nice size
>   "size of matrix"
[1] "size of matrix"
>   n
[1] 5
>  # Hilbert matrix is easy, but not exactly represented
>   A <- matrix(0,n,n)              # initialize to zero
>     for (i in 1:n) {              # fill element by element
+       for (j in 1:n) { A[i,j] <- 1/(i+j-1) } }
>   "Hilbert matrix"
[1] "Hilbert matrix"
>   A
          [,1]      [,2]      [,3]      [,4]      [,5]
[1,] 1.0000000 0.5000000 0.3333333 0.2500000 0.2000000
[2,] 0.5000000 0.3333333 0.2500000 0.2000000 0.1666667
[3,] 0.3333333 0.2500000 0.2000000 0.1666667 0.1428571
[4,] 0.2500000 0.2000000 0.1666667 0.1428571 0.1250000
[5,] 0.2000000 0.1666667 0.1428571 0.1250000 0.1111111
>  # inverse of matrix has integer values -- exact representation
>  # save some constants to make calculation easier
>   Is <- 1:n
>   bi <- choose(Is+n-1,Is-1)       # save these (use with i)
>   bj <- choose(n,Is)              # and these (use with j)
>   B <- matrix(0,n,n)              # initialize to zero
>     sign <- -1
>     for (i in 1:n) {
+       for (j in 1:n) {
+         sign <- -sign              # (-1)**(i+j)
+         B[i,j] <- sign*j*choose(i+j-2,i-1)*bi[i]*choose(j+n-1,n-i)*bj[j]
+                       } }   
>   "inverse of Hilbert matrix"
[1] "inverse of Hilbert matrix"
>   B
      [,1]   [,2]    [,3]    [,4]   [,5]
[1,]    25   -300    1050   -1400    630
[2,]  -300   4800  -18900   26880 -12600
[3,]  1050 -18900   79380 -117600  56700
[4,] -1400  26880 -117600  179200 -88200
[5,]   630 -12600   56700  -88200  44100
>  # is it really the inverse
>   "A %*% B"
[1] "A %*% B"
>   A %*% B
              [,1]          [,2]          [,3]          [,4]          [,5]
[1,]  1.000000e-00 -1.398881e-13  6.297185e-13  2.659206e-12 -1.329603e-12
[2,] -2.003953e-14  1.000000e+00 -2.343903e-12  2.635225e-12 -1.317613e-12
[3,]  2.342571e-14 -1.274536e-13  1.000000e+00 -2.938094e-12  5.595524e-13
[4,]  0.000000e+00 -2.273737e-13  9.094947e-13  1.000000e+00  0.000000e+00
[5,] -1.809664e-14  3.050893e-13 -3.499423e-13  2.362555e-12  1.000000e-00
>  # compute p = inf norms (row sum norm)
>   normA <- max( apply(abs(A),1,sum) )
>   "infinity norm of Hilbert matrix"
[1] "infinity norm of Hilbert matrix"
>   normA
[1] 2.283333
>   normB <- max( apply(abs(B),1,sum) )
>   "infinity norm of inverse of Hilbert matrix"
[1] "infinity norm of inverse of Hilbert matrix"
>   normB
[1] 413280
>   "kappa and reciprocal"
[1] "kappa and reciprocal"
>   kappa <- normA*normB
>   kappa
[1] 943656
>   1/kappa
[1] 1.059708e-06
>  # 
>   rm(list=ls())
>