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}} ______________________________________________ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html