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-11x <- 2.04

# first order derivative

numd1 <- deriv(x,f=func1)

exact <- 10*cos(10*x) + exp(-x)

c(numd1,exact,numd1/exact-1)

# 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-11exact <- -100*sin(10*x) - exp(-x)

c(numd2,exact,numd2/exact-1)

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: rvaradhan at jhmi.edu

--------------------------------------------------------------------------

-----Original Message-----

From: r-help-bounces at stat.math.ethz.ch [mailto:r-help-

bounces at stat.math.ethz.ch] On Behalf Of Peter Dalgaard

Sent: Thursday, May 05, 2005 8:51 PM

To: Uzuner, Tolga

Cc: 'r-help at stat.math.ethz.ch'; 'Berton Gunter'

Subject: Re: [R] Numerical Derivative / Numerical Differentiation of unkno

wnfunct ion

"Uzuner, Tolga" <tolga.uzuner@csfb.com> writes:

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

~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk) FAX: (+45) 35327907

______________________________________________

R-help at 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

From: r-help-bounces at stat.math.ethz.ch [mailto:r-help-

bounces at stat.math.ethz.ch] On Behalf Of Peter Dalgaard

Sent: Thursday, May 05, 2005 8:51 PM

To: Uzuner, Tolga

Cc: 'r-help at stat.math.ethz.ch'; 'Berton Gunter'

Subject: Re: [R] Numerical Derivative / Numerical Differentiation of unkno

wnfunct ion

"Uzuner, Tolga" <tolga.uzuner@csfb.com> 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:gunter.berton at gene.com]

Sent: 05 May 2005 23:34

To: 'Uzuner, Tolga'; r-help at 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 orhelp 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:gunter.berton at gene.com]

Sent: 05 May 2005 23:34

To: 'Uzuner, Tolga'; r-help at 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).

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

~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk) FAX: (+45) 35327907

______________________________________________

R-help at 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