# numdif.r                      October 2007, 2009
 # 
 # numerical differentiation demonstration
 #
 # compute first and second derivatives of log(gamma(x)) at x=1/2
 #
  h <- 2**(-(2:32))
  x <- 1/2
  fx <- lgamma(x)
  fxph <- lgamma(x+h)
  fxmh <- lgamma(x-h)
 # now numerical derivatives
  fp1 <- (fxph - fx)/h        # one sided -- first difference
  fp2 <- (fxph - fxmh)/(2*h)  # central difference
 # second derivatives -- how are these different?
  fpp1 <- (fxph - 2*fx + fxmh)/(h*h)
  fpp2 <- ((fxph - fx) - (fx - fxmh))/(h*h)
  abs(fpp1-fpp2)
 #
 #  targets
  digamma(x)                  # first derivative
  trigamma(x)                 # second derivative
 #  table
  cbind(h,fp1,fp2,fpp1,fpp2)
 #  done
  rm(list=ls())
  q()