swh:1:snp:81eadaa089e8253d8469bcef66aa332632c6c669
Tip revision: 113add4c7a9bfd37466a9b41cf2dc060b7b87091 authored by HwB on 06 December 2012, 00:00:00 UTC
version 1.3.1
version 1.3.1
Tip revision: 113add4
complexstep.R
##
## c o m p l e x s t e p . R Complex Step Derivation
##
complexstep <- function(f, x0, h = 1e-20, ...) {
stopifnot(is.numeric(x0), is.numeric(h), h < 1e-15)
fun <- match.fun(f)
f <- function(x) fun(x, ...)
try(fx0hi <- f(x0 + h*1i))
if (inherits(fx0hi, "try-error"))
stop("Function 'f' does not accept complex arguments.")
if (!is.complex(fx0hi) || !is.real(f(x0))) {
# apply Richardson's method
f_csd <- numderiv(f, x0, h = 0.1)$df
warning("Some maginary part is zero: applied Richardson instead.")
} else {
# apply complex-step method
f_csd <- Im(fx0hi) / h # Im(f(x0 + h * 1i)) / h
}
return(f_csd)
}
complexstepJ <- function(f, x0, h = 1e-20, ...) {
fun <- match.fun(f)
f <- function(x) fun(x, ...)
z <- f(x0)
n <- length(x0); m <- length(z)
J <- matrix(NA, nrow = m, ncol = n)
for (k in 1:n) {
x1 <- x0
x1[k] <- x1[k] + h * 1i
J[ , k] <- Im(f(x1)) / h
}
# drop matrix dimensions if n = m = 1
if (m == 1 && n == 1) J <- J[,]
return(J)
}