wimlib.net Git - wimlib/blob - tools/log2_interpolation.r

   1 #
   2 # This R language script finds a degree 2 interpolating polynomial for f =
   3 # log2(x) over the interval [1, 2).
   4 #
   5
   6 f = function(x) { log2(x) }
   7
   8 get.chebyshev.points = function(a, b, n)
   9 {
  10     (a+b)/2 + (b-a)/2 * cos(((2*(1:n)-1)*pi)/(2*n))
  11 }
  12
  13 build.vandermonde.matrix = function(x)
  14 {
  15     n = length(x)
  16     V = matrix(nrow=n, ncol=n)
  17     for (j in 1:n)
  18         V[,j] = x^(j-1)
  19     return(V)
  20 }
  21
  22 evaluate.polynomial = function(coeffs, x)
  23 {
  24     y = coeffs[length(coeffs)]
  25     for (i in (length(coeffs)-1):1)
  26         y = y*x + coeffs[i]
  27     return(y)
  28 }
  29
  30 x.plot = seq(1, 2, length=1000)
  31 x.chebychev = get.chebyshev.points(1, 2, 3)
  32 V = build.vandermonde.matrix(x.chebychev)
  33 coeffs.a = solve(V, f(x.chebychev))
  34 coeffs.a = c(coeffs.a)
  35 polynomial.interp = function(x) { evaluate.polynomial(coeffs.a, x) }
  36 cat("Coefficients of degree 2 interpolating polynomial:\n")
  37 options(digits=10)
  38 cat(coeffs.a, "\n")
  39 pdf("polynomial-interp.pdf")
  40 plot(x.plot, f(x.plot), col="black", type="l", xlab="x", ylab="y",
  41          main="f(x) and interpolating polynomial approximation")
  42 points(x.chebychev, f(x.chebychev), pch=19, col="red")
  43 lines(x.plot, polynomial.interp(x.plot), col="blue")
  44 legend("topleft", pch=c(NA, NA, 19),
  45            col=c("black", "blue", "red"), lty=c(1,1,0),
  46        legend=c("f(x)", "Interpolating polynomial", "Chebychev points"),
  47        inset=0.07)