>  # 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)
 [1] 8.881784e-16 3.552714e-15 0.000000e+00 0.000000e+00 0.000000e+00
 [6] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
[11] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
[16] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
[21] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
[26] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
[31] 0.000000e+00
>  #
>  #  targets
>   digamma(x)                  # first derivative
[1] -1.96351
>   trigamma(x)                 # second derivative
[1] 4.934802
>  #  table
>   cbind(h,fp1,fp2,fpp1,fpp2)
                 h       fp1       fp2       fpp1       fpp2
 [1,] 2.500000e-01 -1.476336 -2.169483   5.545177   5.545177
 [2,] 1.250000e-01 -1.692284 -2.008978   5.067110   5.067110
 [3,] 6.250000e-02 -1.819352 -1.974565   4.966841   4.966841
 [4,] 3.125000e-02 -1.889025 -1.966255   4.942750   4.942750
 [5,] 1.562500e-02 -1.925627 -1.964195   4.936785   4.936785
 [6,] 7.812500e-03 -1.944403 -1.963681   4.935298   4.935298
 [7,] 3.906250e-03 -1.953914 -1.963553   4.934926   4.934926
 [8,] 1.953125e-03 -1.958702 -1.963521   4.934833   4.934833
 [9,] 9.765625e-04 -1.961103 -1.963513   4.934810   4.934810
[10,] 4.882812e-04 -1.962306 -1.963511   4.934804   4.934804
[11,] 2.441406e-04 -1.962908 -1.963510   4.934803   4.934803
[12,] 1.220703e-04 -1.963209 -1.963510   4.934802   4.934802
[13,] 6.103516e-05 -1.963359 -1.963510   4.934802   4.934802
[14,] 3.051758e-05 -1.963435 -1.963510   4.934802   4.934802
[15,] 1.525879e-05 -1.963472 -1.963510   4.934803   4.934803
[16,] 7.629395e-06 -1.963491 -1.963510   4.934801   4.934801
[17,] 3.814697e-06 -1.963501 -1.963510   4.934807   4.934807
[18,] 1.907349e-06 -1.963505 -1.963510   4.934814   4.934814
[19,] 9.536743e-07 -1.963508 -1.963510   4.934937   4.934937
[20,] 4.768372e-07 -1.963509 -1.963510   4.935059   4.935059
[21,] 2.384186e-07 -1.963509 -1.963510   4.937500   4.937500
[22,] 1.192093e-07 -1.963510 -1.963510   4.937500   4.937500
[23,] 5.960464e-08 -1.963510 -1.963510   4.968750   4.968750
[24,] 2.980232e-08 -1.963510 -1.963510   5.125000   5.125000
[25,] 1.490116e-08 -1.963510 -1.963510   6.000000   6.000000
[26,] 7.450581e-09 -1.963510 -1.963510   6.000000   6.000000
[27,] 3.725290e-09 -1.963510 -1.963510   8.000000   8.000000
[28,] 1.862645e-09 -1.963510 -1.963510   0.000000   0.000000
[29,] 9.313226e-10 -1.963510 -1.963510 128.000000 128.000000
[30,] 4.656613e-10 -1.963510 -1.963510 512.000000 512.000000
[31,] 2.328306e-10 -1.963510 -1.963510   0.000000   0.000000
>  #  done
>   rm(list=ls())
>   q()