https://github.com/cran/pracma
Tip revision: 03698027c2d84118bd0c53c4a9a5b5d23676f388 authored by HwB on 01 October 2012, 00:00:00 UTC
version 1.2.0
version 1.2.0
Tip revision: 0369802
numderiv.R
##
## n u m d e r i v . R Richardson Numerical Derivative
##
numderiv <- function(f, x0, maxiter = 16, h = 1, ...,
accuracy = .Machine$double.eps)
{
if (length(x0) != 1 || !is.numeric(x0))
stop("Argument 'x0' must be a numeric scalar.")
fun <- match.fun(f)
f <- function(x) fun(x, ...)
if (length(f(x0)) != 1)
stop("Function 'f' must be a univariate function of one variable.")
tol <- accuracy^(2/3)
eps <- .Machine$double.eps
err <- 1; relerr <- 1
j <- 1
D <- matrix(0, nrow = maxiter, ncol = maxiter)
D[1, 1] <- (f(x0+h) - f(x0-h))/(2*h)
while (relerr > tol && err > tol && j < maxiter) {
h <- h/2
D[j+1, 1] <- (f(x0+h) - f(x0-h))/(2*h)
for (k in 1:j) {
D[j+1, k+1] <- D[j+1, k] + (D[j+1,k] - D[j,k]) / (4^k - 1)
}
err <- abs(D[j+1,j+1] - D[j,j])
relerr <- 2*err / (abs(D[j+1,j+1]) + abs(D[j,j]) + eps)
j <- j + 1
}
if (j >= maxiter)
warning("Maximum number of iterations reached.")
return(list(df = D[j, j], err = err, relerr = relerr, n = j))
}
numdiff <- function(f, x, maxiter = 16, h = 1, ...,
accuracy = .Machine$double.eps) {
if (!is.vector(x, mode = "numeric"))
stop("Argument 'x' must be a numeric vector.")
ndf <- function(xx) numderiv(f, xx, maxiter = maxiter, h = h, ...,
accuracy = .Machine$double.eps)$df
sapply(x, ndf)
}