>  # chex35.r                         J F Monahan, August 2007
>  #                                  revised September 2009
>  # Example 3.5 Condition of Hilbert matrix
>  #
>  # set n = size of matrix
>   for (n in(3:8)) {               # not too big, but 2 is a problem 
+   print("size of matrix")
+   print(n)
+  # Hilbert matrix is easy, but not exactly represented
+   A <- matrix(0,n,n)              # initialize to zero
+   "Hilbert matrix"
+   A <- 1/(row(A)+col(A)-1)
+  # inverse of matrix has integer values -- exact representation
+  # save some constants to make calculation easier
+   Is <- 1:n                       # needs n > 2 to work
+   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]
+                       } }   
+   print("inverse of Hilbert matrix")
+   print(B)
+  # is it really the inverse?
+   "A %*% B"
+   isitI <- A %*% B
+   print(c(max(abs(diag(isitI)-1)),max(abs(isitI[diag(n)==0]))))
+  # compute p = inf norms (row sum norm)
+   normA <- max( apply(abs(A),1,sum) )
+   normB <- max( apply(abs(B),1,sum) )
+   print("infinity norm of Hilbert matrix and inverse")
+   print(c(normA,normB))
+   print("true kappa, reciprocal, and estimate")
+   tkappa <- normA*normB
+   print(c(tkappa,1/tkappa,kappa(B,norm="I",exact=FALSE)))
+                                      } # end of loop on n
[1] "size of matrix"
[1] 3
[1] "inverse of Hilbert matrix"
     [,1] [,2] [,3]
[1,]    9  -36   30
[2,]  -36  192 -180
[3,]   30 -180  180
[1] 0 0
[1] "infinity norm of Hilbert matrix and inverse"
[1]   1.833333 408.000000
[1] "true kappa, reciprocal, and estimate"
[1] 7.480000e+02 1.336898e-03 1.944313e+01
[1] "size of matrix"
[1] 4
[1] "inverse of Hilbert matrix"
     [,1]  [,2] [,3]  [,4]
[1,]   16  -120  240  -140
[2,]  120 -1200 2700 -1680
[3,]  240 -2700 6480 -4200
[4,]  140 -1680 4200 -2800
[1] 2750 4800
[1] "infinity norm of Hilbert matrix and inverse"
[1]     2.083333 13620.000000
[1] "true kappa, reciprocal, and estimate"
[1] 2.837500e+04 3.524229e-05 1.088215e+02
[1] "size of matrix"
[1] 5
[1] "inverse of Hilbert matrix"
      [,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
[1] 0.000000e+00 3.637979e-12
[1] "infinity norm of Hilbert matrix and inverse"
[1] 2.283333e+00 4.132800e+05
[1] "true kappa, reciprocal, and estimate"
[1] 9.436560e+05 1.059708e-06 6.099623e+02
[1] "size of matrix"
[1] 6
[1] "inverse of Hilbert matrix"
     [,1]    [,2]    [,3]     [,4]    [,5]     [,6]
[1,]   36    -630    3360    -7560    7560    -2772
[2,]  630  -14700   88200  -211680  220500   -83160
[3,] 3360  -88200  564480 -1411200 1512000  -582120
[4,] 7560 -211680 1411200 -3628800 3969000 -1552320
[5,] 7560 -220500 1512000 -3969000 4410000 -1746360
[6,] 2772  -83160  582120 -1552320 1746360  -698544
[1] 1466432 2787120
[1] "infinity norm of Hilbert matrix and inverse"
[1]        2.45 11865420.00
[1] "true kappa, reciprocal, and estimate"
[1] 2.907028e+07 3.439939e-08 3.421450e+03
[1] "size of matrix"
[1] 7
[1] "inverse of Hilbert matrix"
       [,1]     [,2]      [,3]      [,4]       [,5]       [,6]      [,7]
[1,]     49    -1176      8820    -29400      48510     -38808     12012
[2,]  -1176    37632   -317520   1128960   -1940400    1596672   -504504
[3,]   8820  -317520   2857680 -10584000   18711000  -15717240   5045040
[4,] -29400  1128960 -10584000  40320000  -72765000   62092800 -20180160
[5,]  48510 -1940400  18711000 -72765000  133402500 -115259760  37837800
[6,] -38808  1596672 -15717240  62092800 -115259760  100590336 -33297264
[7,]  12012  -504504   5045040 -20180160   37837800  -33297264  11099088
[1] 0.000000e+00 4.656613e-10
[1] "infinity norm of Hilbert matrix and inverse"
[1] 2.592857e+00 3.799650e+08
[1] "true kappa, reciprocal, and estimate"
[1] 9.851949e+08 1.015028e-09 2.132766e+04
[1] "size of matrix"
[1] 8
[1] "inverse of Hilbert matrix"
       [,1]      [,2]      [,3]        [,4]       [,5]        [,6]       [,7]
[1,]     64     -2016     20160      -92400     221760     -288288     192192
[2,]   2016    -84672    952560    -4656960   11642400   -15567552   10594584
[3,]  20160   -952560  11430720   -58212000  149688000  -204324120  141261120
[4,]  92400  -4656960  58212000  -304920000  800415000 -1109908800  776936160
[5,] 221760 -11642400 149688000  -800415000 2134440000 -2996753760 2118916800
[6,] 288288 -15567552 204324120 -1109908800 2996753760 -4249941696 3030051024
[7,] 192192 -10594584 141261120  -776936160 2118916800 -3030051024 2175421248
[8,]  51480  -2882880  38918880  -216216000  594594000  -856215360  618377760
           [,8]
[1,]     -51480
[2,]   -2882880
[3,]  -38918880
[4,] -216216000
[5,] -594594000
[6,] -856215360
[7,] -618377760
[8,] -176679360
[1] 1155536384 2201223024
[1] "infinity norm of Hilbert matrix and inverse"
[1] 2.717857e+00 1.246305e+10
[1] "true kappa, reciprocal, and estimate"
[1] 3.387279e+10 2.952222e-11 1.501650e+05
>  # 
>   rm(list=ls())
>