# chex36.r                         J F Monahan, March 2010
 #                                  
 # Example 3.6 Condition of Interpolation problem
 #
 # function to compute infinity norm
  inorm <- function(A) max(apply(abs(A),1,sum))
 # set n = size of matrix
  n <- 4                          # four is big enough
  print("size of matrix")
  print(n)
  for( method in(1:2)){            # loop on method
  if( method == 1) x <- (2*(1:n)-1-n)/(2*n)      # equally spaced abscissas
  if( method == 2) x <- cos((2*(1:n)-1)*pi/(2*n))# Chebyshev abscissas
 # rhs vector b -- simple function f(x)=1/(1+x^2)
  b <- 1/(1+x*x)
 # make matrix A
  A <- matrix(0,n,n)
  A[,1] <- 1
  for( j in(2:n) ){ A[,j] <- x*A[,(j-1)] }
  Ai <- solve(A)
 # condition number
  print("condition number, reciprocal")
  print(c(inorm(A)*inorm(Ai),1/(inorm(A)*inorm(Ai))))
  
 # solution
  Aib <- Ai %*% b
  print(Aib)
  normAib <- max(abs(Aib))
  print("k, norm soln, eps, norm difference")
 # perturb rhs and solve
  for ( k in(1:5) ){
  eps <- 10^(-8-k)
  bd <- b + eps                # b + epsilon
 # perturbed solution
  Aibd <- Ai %*% bd
  normdiff <- max(abs(Aib-Aibd))
  print(c(k,normAib,eps,normdiff))  
                         }     # end of k loop
                         }     # end of method loop
 # 
  rm(list=ls())