# splint.s                         J F Monahan, June 2001
 # interpolatory cubic splines -- two methods
 # 1) natural spline
 # 2) with derivatives at endpoints
 npts <- 100                  # evaluation points for comparison
 n <- 7                       # interpolation points
 xk <- 4*( 2*seq(n)-n-1 )/(n-1) 
 "interpolation points, aka 'knots' "
 xk
 "function values at knots"
 z <- 1/(1 + xk*xk)
 z 
# natural spline 
 splstn <- function(xk,z,n) {
   h <- c( 0, (xk[2:n] - xk[1:(n-1)]) )
   s <- c( 0, (z[2:n] - z[1:(n-1)] ) / h[2:n] )
   d <- 6*( s[3:n] - s[2:(n-1)] ) / ( h[3:n] + h[2:(n-1)] )
   M <- c(0, d, 0 )
   Dd <- rep(2,n)
   Dsb <- c( NA, h[2:(n-1)] / ( h[3:n] + h[2:(n-1)]) , 0 )
   Dsp <- c( 0, h[3:n] / ( h[3:n] + h[2:(n-1)]), NA )
 # LU factorization of (Dsb, Dd, Dsp)
   for( k in 2:n ) {
     Dsb[k] <- Dsb[k] / Dd[k-1]
     Dd[k] <- Dd[k] - Dsb[k]*Dsp[k-1]
                   }
 # solve system
   for( k in 2:n ) M[k] <- M[k] - Dsb[k]*M[k-1]
   M[n] <- M[n]/Dd[n]
   for( k in (n-1):1 ) M[k] <- (M[k] - Dsp[k]*M[k+1])/Dd[k]
   return( list(xk=xk,z=z,n=n,M=M,h=h) )
                            } # end of function splstn
# spline with derivative condition
 splstd <- function(xk,z,n,d1,dn) {
   h <- c( 0, (xk[2:n] - xk[1:(n-1)]) )
   s <- c( 0, (z[2:n] - z[1:(n-1)] ) / h[2:n] )
   d <- 6*( s[3:n] - s[2:(n-1)] ) / ( h[3:n] + h[2:(n-1)] )
   M <- c(0, d, 0 )
   M[1] <- 6*(s[2]-d1)/h[2]         # s(2) and h(2)
   M[n] <- 6*(dn-s[n])/h[n]     # s(n) and h(n)
   Dd <- rep(2,n)
   Dsb <- c( NA, h[2:(n-1)] / ( h[3:n] + h[2:(n-1)]) , 1 )
   Dsp <- c( 1, h[3:n] / ( h[3:n] + h[2:(n-1)]), NA )
 # LU factorization of (Dsb, Dd, Dsp)
   for( k in 2:n ) {
     Dsb[k] <- Dsb[k] / Dd[k-1]
     Dd[k] <- Dd[k] - Dsb[k]*Dsp[k-1]
                   }
 # solve system
   for( k in 2:n ) M[k] <- M[k] - Dsb[k]*M[k-1]
   M[n] <- M[n]/Dd[n]
   for( k in (n-1):1 ) M[k] <- (M[k] - Dsp[k]*M[k+1])/Dd[k]
   return( list(xk=xk,z=z,n=n,M=M,h=h) )
                            } # end of function splstd
# evaluate spline at x vector
  splev <- function(x,splinfo) {
     j <-  1 + as.numeric( cut(x,splinfo$xk) )
     du <- splinfo$xk[j] - x
     dl <- x - splinfo$xk[(j-1)]
     Mu <- splinfo$M[j]
     Ml <- splinfo$M[(j-1)]
     hj <- splinfo$h[j]   
     ( (Ml*du*du*du + Mu*dl*dl*dl)/6
                   + (splinfo$z[(j-1)]-Ml*hj*hj/6 )*du
                   + (splinfo$z[j]-Mu*hj*hj/6 )*dl )/hj
                                }  # end of function splev 
 #
 # first get info on natural spline
   splnat <- splstn(xk,z,n)
   "natural spline -- second derivs at knots"
   splnat$M
   splnat$n
 # now get info on spline with deriv condition
   splder <- splstd(xk,z,n,-8/289,8/289) # deriv at endpoints
   "spline(with deriv cond) -- second derivs at knots"
   splder$M
   splder$n
 #
 # compare two methods
 xi <- 4*( 2*seq(npts)-npts-1 )/(npts-1)
 fxi <- 1/(1 + xi*xi )
 fni <- splev(xi,splnat)   # natural spline
 fdi <- splev(xi,splder)   # with deriv condition

 xy <- paste( 'xi', format(xi),'fxi',format(fxi),
               'nat',format(fni),'der',format(fdi) )
 cat(xy,sep="\n")
 rm(list=ls())