# R/complexstep.R In pracma: Practical Numerical Math Functions

```##
##  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 appear to accept complex arguments.")

if (!is.complex(fx0hi) || !is.double(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)
}

grad_csd <- function(f, x0, h = 1e-20, ...) {
fun <- match.fun(f)
f <- function(x) fun(x, ...)

z <- f(x0)
n <- length(x0); m <- length(z)
if (m > 1)
stop("Funktion 'f' does not return a scalar; call 'Jacobian_csd'.")

G <- rep(NA, n)
for (k in 1:n) {
x1 <- x0
x1[k] <- x1[k] + h * 1i
G[k] <- Im(f(x1)) / h
}
return(G)
}

jacobian_csd <- 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)
}

hessian_csd <- function(f, x0, h = 1e-20, ...) {
fun <- match.fun(f)
f <- function(x) fun(x, ...)

z <- f(x0)
n <- length(x0); m <- length(z)
if (m > 1)
stop("Funktion 'f' does not return a scalar, as needed for Hessian.")

H <- matrix(NA, nrow = n, ncol = n)
for (i in 1:n) {
fi <- function(x) {
x[i] <- x[i] + h*1i
Im(f(x)) / h
}
for (j in 1:n) {
ff <- function(x) {
xx <- x0
xx[j] <- x
fi(xx)
}
H[i, j] <- numderiv(ff, x0[j])\$df
}
}
return(H)
}

laplacian_csd <- function(f, x0, h = 1e-20, ...) {
fun <- match.fun(f)
f <- function(x) fun(x, ...)

z <- f(x0)
n <- length(x0); m <- length(z)
if (m > 1)
stop("Funktion 'f' does not return a scalar, as needed for Lapacian.")

L <- rep(NA, n)
for (i in 1:n) {
fi <- function(x) {
x[i] <- x[i] + h*1i
Im(f(x)) / h
}
ff <- function(x) {
xx <- x0
xx[i] <- x
fi(xx)
}
L[i] <- numderiv(ff, x0[i])\$df
}
return(sum(L))
}
```

## Try the pracma package in your browser

Any scripts or data that you put into this service are public.

pracma documentation built on March 18, 2022, 5:12 p.m.