swh:1:snp:81eadaa089e8253d8469bcef66aa332632c6c669
Tip revision: c1688b374d201c13fb40b4dda2d2a89e34b94ec6 authored by Hans W. Borchers on 23 January 2021, 09:10:02 UTC
version 2.3.3
version 2.3.3
Tip revision: c1688b3
shubert.R
##
## s h u b e r t . R Shubert-Piyawskii Method
##
shubert <- function(f, a, b, L, crit = 1e-04, nmax = 1000) {
stopifnot(is.numeric(a), is.numeric(b),
length(a) == 1, length(b) == 1)
fun <- match.fun(f); f <- function(x) fun(x)
nf <- 0
x0 <- (a + b)/2.0
y0 <- f(x0); nf <- nf + 1
xmax <- x0; ymax <- y0
fmax <- y0 + L*(b - a)/2.0
X <- numeric(nmax); Z <- numeric(nmax)
X[1] <- b; Z[1] <- y0 + fmax
X[2] <- a; Z[2] <- y0 + fmax
n <- 2
while ((fmax - ymax) > crit && n < nmax) {
xn <- X[n]; zn <- Z[n]
yn <- f(xn); nf <- nf + 1
if (yn > ymax) {
xmax <- xn; ymax <- yn
}
zL <- (zn + yn)/2.0; zR <- zL
xL <- xn - (zn - yn)/2.0/L
xR <- xn + (zn - yn)/2.0/L
i1 <- 0; i2 <- 0
if (xL >= a && xL <= b) i1 <- 1
if (xR >= a && xR <= b) i2 <- 1
if (i1 == 1 && i2 == 0) {
X[n] <- xL; Z[n] <- zL
} else if (i1 == 0 && i2 == 1) {
X[n] <- xR; Z[n] <- zR
} else if (i1 == 1 && i2 == 1) {
X[n] <- xL; Z[n] <- zL
X[n+1] <- xR; Z[n+1] <- zR
n <- n + 1
}
Zo <- order(Z[1:n])
Z[1:n] <- Z[Zo]; X[1:n] <- X[Zo]
fmax <- Z[n]
}
return(list(xopt = xmax, fopt = fmax, nopt = n))
}