R code for polynomial approximation example

poly-approx.r

# test polynomial interpolation code by interpolating points on a semi-circle 
# using a different number of points each time.

# contents:
# 1) make.circle(degree) is a function that generates 
#    n=degree+1 data points (x,y) on a semicircle. this is the "data" we 
#    want to interpolate.
# 2) pinterp() makes polynomial interpolation on "data" of degree "degree" 
#    using degree+1 parameters 
# 3) run() is a wrapper that runs the interpolation

# note: polynomial interpolation uses n=degree+1 points and T=degree+1 parameters.
#
# output:
# plot of the interpolation overlaid by the used data points

rm(list=ls(all=TRUE))
setwd("~/dropbox/Rstuff/splines")

require(plotrix)
# this package contains a function to draw a circle. 
# you could generate any other data as well.

# task: run this approximation for 3,7,11 and 15 sample points.
degrees <- c(2,6,10,14)


##################
# define functions
##################
make.circle <- function(degree){
# makes a vector of degree+1 points on a semi-circle
    plot(-2:2,-2:2,type="n",xlab="",ylab="",main="make.circle output. disregard.")
    vals <- draw.circle(0,0,radius=1,nv=2*(degree)+1)
    data <- array(unlist(vals),dim=c(2*degree+1,2))
    data <- data[data[ ,2]>=0, ]    # take only upper half of circle§
    return(data)
}

pinterp <- function(degree,params,data,t){
# polynomial interpolation on data at point t

# has a list "x" whose entries are each a 
# matrix for points at t at approximation of degree
# deg. x gets shorter with increasing degree. at final entry
# degree+1, x is (1,2)
    x <- vector("list",degree+1)
    x[[1]] <- data  # matrix 1 holds data for degree zero
                    # each row of x[[j]] 
                    # represents (x,y) coords of a point
    n <- degree+1   # number of points we're interpolating on
    for (ideg in 2:n){  # recursively build interpolations of increasing degree
        nn <- n-ideg+1
        x[[ideg]] <- array(data=NA,dim=c(nn,2))
        for (i in 1:nn){
            x[[ideg]][i, ] <- 
            (params[i+ideg-1] - t)/(params[i+ideg-1] - params[i])*x[[ideg-1]][i, ] + 
            (t - params[i])/(params[i+ideg-1] - params[i])*x[[ideg-1]][i+1, ]
        }
    }
    return(x[[degree+1]])   # return only final degree values,i.e. result
}

run <- function(degrees,t){
    nd <- length(degrees)
    rval <- vector("list",nd)
    origdata <- vector("list",nd)
    for (j in 1:nd){
        degree <- degrees[j]
        params <- seq(0,5,length=degree+1)
        nt <- length(t)
        rval[[j]] <- array(data=NA,dim=c(nt,2))
        data <- make.circle(degree)
        origdata[[j]] <- data
        for (it in 1:nt) rval[[j]][it, ] <- pinterp(degree,params,data,t[it])
    }
    return(list(interpol=rval,orig=origdata))
}

###################
# produce output
###################

outrun <- run(degrees,t=seq(1,5,length=20))
rval <- outrun$interpol
orig <- outrun$orig
png(file="graphs/ex1.png")
par(mfrow=c(2,2))
for (j in 1:length(degrees)){
    plot(x=orig[[j]][ ,1],y=orig[[j]][ ,2],type="o",
    main=paste("degree: ",degrees[j],", points: ",
    degrees[j]+1),xlab="x",ylab="y",ylim=c(0,1.5))
    lines(x=rval[[j]][ ,1],y=rval[[j]][ ,2],col="red")    
}
dev.off()
par(mfrow=c(1,1))