Hi, There is a nice article by Fornberg and Sloan (Acta Numerica 1994) on various order accuracy (Taylor-series based) approximations for different order derivatives. I had coded a couple of these in R for first and second order derivatives, with truncation errors of orders 2 and 4.
Here is the code and a simple example demonstrating its accuracy for a sharply oscillating function: ############################################ # A function to compute highly accurate first- and second-order derivatives # From Fornberg and Sloan (Acta Numerica, 1994, p. 203-267; Table 1, page 213) deriv <- function(x, fun, h=NULL, order=1, accur=4) { macheps <- .Machine$double.eps if (order==1) { if(is.null(h)) h <- macheps^(1/3)* abs(x) ifelse (accur==2, w <- c(-1/2,1/2), w <- c(1/12,-2/3, 2/3,-1/12)) ifelse (accur==2, xd <- x + h*c(-1,1), xd <- x + h*c(-2,-1,1,2)) return(sum(w*fun(xd))/h) } else if (order==2) { if(is.null(h)) h <- macheps^(1/4)* abs(x) ifelse (accur==2, w <- c(1,-2,1), w <- c(-1/12,4/3,-5/2,4/3,-1/12)) ifelse (accur==2, xd <- x + h*c(-1,0,1), xd <- x + h*c(-2,-1,0,1,2)) return(sum(w*fun(xd))/h^2) } } ############################################ func1 <- function(x){ sin(10*x) - exp(-x) } ############################################# > curve(func1,from=0,to=5) > x <- 2.04 # first order derivative > numd1 <- deriv(x,f=func1) > exact <- 10*cos(10*x) + exp(-x) > c(numd1,exact,numd1/exact-1) [1] 3.335371e-01 3.335371e-01 1.981793e-11 # second order derivative > numd2 <- deriv(x,f=func1,order=2) > exact <- -100*sin(10*x) - exp(-x) > c(numd2,exact,numd2/exact-1) [1] -1.001093e+02 -1.001093e+02 -2.300948e-11 Hope this is helpful. Best, Ravi. -------------------------------------------------------------------------- Ravi Varadhan, Ph.D. Assistant Professor, The Center on Aging and Health Division of Geriatric Medicine and Gerontology Johns Hopkins University Ph: (410) 502-2619 Fax: (410) 614-9625 Email: [EMAIL PROTECTED] -------------------------------------------------------------------------- > -----Original Message----- > From: [EMAIL PROTECTED] [mailto:r-help- > [EMAIL PROTECTED] On Behalf Of Peter Dalgaard > Sent: Thursday, May 05, 2005 8:51 PM > To: Uzuner, Tolga > Cc: 'r-help@stat.math.ethz.ch'; 'Berton Gunter' > Subject: Re: [R] Numerical Derivative / Numerical Differentiation of unkno > wnfunct ion > > "Uzuner, Tolga" <[EMAIL PROTECTED]> writes: > > > Ah... I searched for half an hour for this function... you know, the > > help function in R could really be a lot better... > > > > But wait a minute... looking at this, it appears you have to pass in > > an expression. What if it is an unknown function, where you only > > have a handle to the function, but you cannot see it's > > implementation ? Will this work then ? > > > > -----Original Message----- > > From: Berton Gunter [mailto:[EMAIL PROTECTED] > > Sent: 05 May 2005 23:34 > > To: 'Uzuner, Tolga'; r-help@stat.math.ethz.ch > > Subject: RE: [R] Numerical Derivative / Numerical Differentiation of > > unknown funct ion > > > > > > But... > > > > See ?numericDeriv which already does it via a C call and hence is much > > faster (and probably more accurate,too). > > > > The expression passed to numericDeriv can easily be a call to .C or > similar. > > Actually, numericDeriv can get you in trouble if the function is not > smooth enough. It basically just calculates (f(a+d)-f(a))/d where d is > on the order of 1e-7 * a for each parameter. Sometimes a larger d and > a higher order approximation is need to avoid getting stuck in the > rough. > > (Yes, Bill, I do remember that you wanted an R News Programmer's Niche > item from me on this...) > > -- > O__ ---- Peter Dalgaard Blegdamsvej 3 > c/ /'_ --- Dept. of Biostatistics 2200 Cph. N > (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 > ~~~~~~~~~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907 > > ______________________________________________ > 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 ______________________________________________ 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