https://github.com/cran/pracma
Tip revision: 708a2ad382a163d1eef5af0665e3ae2aad200ced authored by HwB on 21 March 2013, 00:00:00 UTC
version 1.4.5
version 1.4.5
Tip revision: 708a2ad
fzero.R
##
## f z e r o . R
##
fzero <- function(f, x0, ..., maxiter = 100, tol = .Machine$double.eps^(1/2)) {
if (!is.numeric(x0) || length(x0) > 2)
stop("Argument 'x0' must be a scalar or a vector of length 2.")
if (length(x0) == 2) {
zero <- zeroin(f, x0, ..., tol = tol)
} else {
if (x0 != 0) dx <- x0/50
else dx <- 1/50
sqrt2 <- sqrt(2)
a <- b <- x0
fa <- fb <- f(x0, ...)
if (fa == 0) return(list(x = a, fval = fa))
iter <- 0
while (fa * fb > 0 && iter < maxiter) {
iter <- iter + 1
dx <- sqrt2 * dx
a <- a - dx
fa <- f(a, ...)
if (fa * fb <= 0) break
b <- b + dx
fb <- f(b, ...)
}
if (iter == maxiter) {
warning("Maximum number of iterations exceeded; no zero found.")
return(list(x = NA, fval = NA))
}
zero <- zeroin(f, c(a, b), ..., tol = tol)
}
x.zero <- zero$root
f.zero <- zero$f.root
return(list(x = x.zero, fval = f.zero))
}
zeroin <- function(f, interval, ..., tol = .Machine$double.eps^(1/2)) {
if (!is.numeric(interval) || length(interval) != 2)
stop("Argument 'interval' must be a numeric vector of length 1.")
a <- interval[1]; fa <- f(a, ...); nf <- 1
if (fa == 0)
return(list(root = a, f.root = fa, f.calls = nf, estim.prec = 0.0))
b <- interval[2]; fb <- f(b, ...); nf <- nf + 1
if (fb == 0)
return(list(root = b, f.root = fb, f.calls = nf, estim.prec = 0.0))
if (sign(fa) == sign(fb))
stop("Function 'f' must change sign on the interval.")
# Initialize
cc <- a
fc <- fa
d <- b - cc
e <- d
# Main loop
while (fb != 0) {
if (sign(fa) == sign(fb)) {
a <- cc; fa <- fc
d <- b - cc; e <- d
}
if (abs(fa) < abs(fb)) {
cc <- b; b <- a; a <- cc
fc <- fb; fb <- fa; fa <- fc
}
# Convergence test and possible exit
m <- 0.5 * (a - b)
tol <- 2.0 * tol * max(abs(b), 1.0)
if (abs(m) <= tol || (fb == 0.0)) {
break
}
# Choose bisection or interpolation
if (abs(e) < tol || abs(fc) <= abs(fb)) {
# Bisection
d <- m
e <- m
} else {
# Interpolation
s <- fb/fc
if (a == cc) {
# Linear interpolation (secant)
p <- 2.0 * m * s
q <- 1.0 - s
} else {
# Inverse quadratic interpolation
q <- fc/fa
r <- fb/fa
p <- s * (2.0 * m * q * (q - r) - (b - cc) * (r - 1.0))
q <- (q - 1.0) * (r - 1.0) * (s - 1.0)
}
if (p > 0) q <- -q else p <- -p
# Is interpolated point acceptable
if (2.0*p < 3.0*m*q - abs(tol*q) && p < abs(0.5*e*q)) {
e <- d
d <- p/q
} else {
d <- m
e <- m
}
}
# Next point
cc <- b
fc <- fb
if (abs(d) > tol) {
b <- b + d
} else {
b <- b - sign(b-a) * tol
}
fb <- f(b, ...); nf <- nf + 1
}
return(list(root = b, f.root = fb, f.calls = nf, estim.prec = m))
}