Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- # 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))
Add Comment
Please, Sign In to add comment