Skip to content

Numerical Derivative / Numerical Differentiation of unknown funct ion

3 messages · Uzuner, Tolga, Bert Gunter, Paul Gilbert

Uzuner, Tolga
# 
Hi,

I have been trying to do numerical differentiation using R. 

I found some old S code using Richardson Extrapolation which I managed to get
to work.

I am posting it here in case anyone needs it.


########################################################################
richardson.grad <- function(func, x, d=0.01, eps=1e-4, r=6, show=F){
# This function calculates a numerical approximation of the first
#   derivative of func at the point x. The calculation
#   is done by Richardson's extrapolation (see eg. G.R.Linfield and
J.E.T.Penny
#   "Microcomputers in Numerical Analysis"). The method should be used if
#   accuracy, as opposed to speed, is important.
#
#  *  modified by Paul Gilbert from orginal code by XINGQIAO LIU.
# CALCULATES THE FIRST ORDER DERIVATIVE 
#     VECTOR using a Richardson  extrapolation.
#
#  GENERAL APPROACH
#     --  GIVEN THE FOLLOWING INITIAL VALUES:
#             INTERVAL VALUE D, NUMBER OF ITERATIONS R, AND
#             REDUCED FACTOR V.
#      - THE FIRST ORDER aproximation to the DERIVATIVE WITH RESPECT TO Xi IS
#
#           F'(Xi)={F(X1,...,Xi+D,...,Xn) - F(X1,...,Xi-D,...,Xn)}/(2*D)
#       
#     --  REPEAT r TIMES  with successively smaller D  and 
#          then apply Richardson extraplolation
#
#  INPUT
#       func    Name of the function.
#       x       The parameters of func.
#       d       Initial interval value (real) by default set to 0.01*x or
#               eps if x is 0.0.
#       r       The number of Richardson improvement iterations.
#       show    If T show intermediate results.
#  OUTPUT
#
#       The gradient vector.

  v <- 2               # reduction factor.
  n <- length(x)       # Integer, number of variables.
  a.mtr <- matrix(1, r, n) 
  b.mtr <- matrix(1, (r - 1), n)
#------------------------------------------------------------------------
# 1 Calculate the derivative formula given in 'GENERAL APPROACH' section.
#   --  The first order derivatives are stored in the matrix a.mtr[k,i], 
#        where the indexing variables k for rows(1 to r),  i for columns
#       (1 to n),  r is the number of iterations, and n is the number of
#       variables.
#-------------------------------------------------------------------------  

  h <- abs(d*x)+eps*(x==0.0)
  for(k in 1:r)  { # successively reduce h                
     for(i in 1:n)  {
         x1.vct <- x2.vct <- x
         x1.vct[i]  <- x[i] + h[i]
         x2.vct[i]  <- x[i] - h[i]
         if(k == 1) a.mtr[k,i] <- (func(x1.vct) - func(x2.vct))/(2*h[i])
         else{
           if(abs(a.mtr[(k-1),i])>1e-20)
                 # some functions are unstable near 0.0              
                 a.mtr[k,i] <- (func(x1.vct)-func(x2.vct))/(2*h[i])
            else  a.mtr[k, i] <- 0
          }
      }
     h <- h/v     # Reduced h by 1/v.
    }	
   if(show)  {
        cat("\n","first order approximations", "\n")		
        print(a.mtr, 12)
    }

#------------------------------------------------------------------------
# 1 Applying Richardson Extrapolation to improve the accuracy of 
#   the first and second order derivatives. The algorithm as follows:
#
#   --  For each column of the 1st and 2nd order derivatives matrix a.mtr,
#       say, A1, A2, ..., Ar, by Richardson Extrapolation, to calculate a
#       new sequence of approximations B1, B2, ..., Br used the formula
#
#          B(i) =( A(i+1)*4^m - A(i) ) / (4^m - 1) ,  i=1,2,...,r-m
#
#             N.B. This formula assumes v=2.
#
#   -- Initially m is taken as 1  and then the process is repeated 
#      restarting with the latest improved values and increasing the 
#      value of m by one each until m equals r-1
#
# 2 Display the improved derivatives for each
#   m from 1 to r-1 if the argument show=T.
#
# 3 Return the final improved  derivative vector.
#-------------------------------------------------------------------------

  for(m in 1:(r - 1)) {		
     for(i in 1:(r - m)) b.mtr[i,]<- (a.mtr[(i+1),]*(4^m)-a.mtr[i,])/(4^m-1)
#     a.mtr<- b.mtr
#     a.mtr<- (a.mtr[2:(r+1-m),]*(4^m)-a.mtr[1:(r-m),])/(4^m-1)
     if(show & m!=(r-1) )  {
        cat("\n","Richarson improvement group No. ", m, "\n")		
        print(a.mtr[1:(r-m),], 12)
      }
   }
a.mtr[length(a.mtr)]
}

## try it out
richardson.grad(function(x){x^3},2)
########################################################################################


Regards,
Tolga Uzuner


==============================================================================
This message is for the sole use of the intended recipient. ...{{dropped}}
Bert Gunter
# 
But...

See ?numericDeriv which already does it via a C call and hence is much
faster (and probably more accurate,too).



-- Bert Gunter
Genentech Non-Clinical Statistics
South San Francisco, CA
 
"The business of the statistician is to catalyze the scientific learning
process."  - George E. P. Box
Paul Gilbert
# 
A current version of this code is in the dseplus bundle in the devel
section of CRAN. It does Richardson's extrapolation, which gives an
accurate numerical estimate at the expense of speed. (It does a very
large number of function evaluations.) That may not be what you want.

Paul Gilbert
Uzuner, Tolga wrote: