R version 2.9.0 (2009-04-17)
Copyright (C) 2009 The R Foundation for Statistical Computing
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 an HTML browser interface to help.
Type 'q()' to quit R.

>  # 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")
[1] "size of matrix"
>   print(n)
[1] 4
>   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
[1] "condition number, reciprocal"
[1] 133.83333333   0.00747198
           [,1]
[1,]  0.9981033
[2,]  0.0000000
[3,] -0.8632244
[4,]  0.0000000
[1] "k, norm soln, eps, norm difference"
[1] 1.000000e+00 9.981033e-01 1.000000e-09 9.999999e-10
[1] 2.000000e+00 9.981033e-01 1.000000e-10 9.999990e-11
[1] 3.000000e+00 9.981033e-01 1.000000e-11 1.000000e-11
[1] 4.000000e+00 9.981033e-01 1.000000e-12 9.999779e-13
[1] 5.000000e+00 9.981033e-01 1.000000e-13 1.001421e-13
[1] "condition number, reciprocal"
[1] 18.63688432  0.05365704
              [,1]
[1,]  9.411765e-01
[2,] -2.775558e-17
[3,] -4.705882e-01
[4,]  5.551115e-17
[1] "k, norm soln, eps, norm difference"
[1] 1.000000000 0.941176471 0.000000001 0.000000001
[1] 2.0000000000 0.9411764706 0.0000000001 0.0000000001
[1] 3.000000e+00 9.411765e-01 1.000000e-11 1.000000e-11
[1] 4.000000e+00 9.411765e-01 1.000000e-12 9.999779e-13
[1] 5.000000e+00 9.411765e-01 1.000000e-13 1.001421e-13
>  # 
>   rm(list=ls())
>