Raw File
fderiv.R
##
##  f d e r i v . R  Numerical Differentiation
##

fderiv <- function(f, x, n = 1, h = 0,
                   method = "central", ...) {
    if (length(x) == 0) return(c())
    if (!is.numeric(x))
        stop("Argument 'x' must be a number or a numeric vector.")
    if (! n %in% c(0, 1, 2, 3, 4))
        stop("The order of the derivative, 'n', can only be 0, 1, 2, 3, or 4.")

    fun <- match.fun(f)
    f <- function(x) fun(x, ...)

    fx <- f(x)
    if (length(fx) != length(x))
        stop("Function 'f' must first be vectorized: Vectorize(f).")

    if (n == 0) return(f(x))

    if (h == 0) {
        H <- .Machine$double.eps^c(1/3, 1/4, 1/5, 1/6)
        h <- H[n]
    }

    if (method == "central") {
        if (n == 1) {
            df <- (f(x+h) - f(x-h)) / (2*h)
        } else if (n == 2) {
            df <- (f(x+h) - 2*f(x) + f(x-h)) / h^2
        } else if (n == 3) {
            df <- (f(x+2*h) - 2*f(x+h) + 2*f(x-h) - f(x-2*h)) / (2*h^3)
        } else if (n == 4) {
            df <- (f(x+2*h) - 4*f(x+h) + 6*f(x) - 4*f(x-h) + f(x-2*h)) / h^4
        }

    } else if (method == "forward") {
        if (n == 1) {
            df <- (-f(x+2*h) + 4*f(x+h) - 3*f(x)) / (2*h)
        } else if (n == 2) {
            df <- (-f(x+3*h) + 4*f(x+2*h) - 5*f(x+h) + 2*f(x)) / h^2
        } else if (n == 3) {
            df <- (-3*f(x+4*h) + 14*f(x+3*h) - 24*f(x+2*h) + 18*f(x+h) - 5*f(x)) / (2*h^3)
        } else if (n == 4) {
            df <- (-2*f(x+5*h) + 11*f(x+4*h) - 24*f(x+3*h) + 26*f(x+2*h) - 14*f(x+h) + 3*f(x)) / h^4
        }

    } else if (method == "backward") {
        if (n == 1) {
            df <- (3*f(x) - 4*f(x-h) + f(x-2*h)) / (2*h)
        } else if (n == 2) {
            df <- (2*f(x) - 5*f(x-h) + 4*f(x-2*h) - f(x-3*h)) / h^2
        } else if (n == 3) {
            df <- (5*f(x) - 18*f(x-h) + 24*f(x-2*h) - 14*f(x-3*h) + 3*f(x-4*h)) / (2*h^3)
        } else if (n == 4) {
            df <- (3*f(x) - 14*f(x-h) + 26*f(x-2*h) - 24*f(x-3*h) + 11*f(x-4*h) - 2*f(x-5*h)) / h^4
        }
        
    } else
        stop("Unknown 'method'; use 'central', 'forward' or 'backward' instead.")

    return(df)
}
back to top