# R/zeroin.R In rootoned: Roots of one-dimensional functions in R-only code

#### Documented in zeroin

```zeroin <- function(f, interval,
tol = .Machine\$double.eps^0.5, maxiter = 1000, trace=FALSE,...){
# trace <- FALSE # to allow for debugging
# epsilon <- .Machine\$double.eps^0.5
epsilon <- .Machine\$double.eps
if (trace) cat("epsilon= ",epsilon,"\n")

# uniroot <- function(f, interval, ...,
#             lower = min(interval), upper = max(interval),
#             f.lower = f(lower, ...), f.upper = f(upper, ...),
#             tol = .Machine\$double.eps^0.25, maxiter = 1000)

#      a zero of the function  f(x)  is computed in the interval ax,bx .
#
#  input..
#
#  ax     left endpoint of initial interval
#  bx     right endpoint of initial interval
#  f      function subprogram which evaluates f(x) for any x in
#         the interval  ax,bx
#  tol    desired length of the interval of uncertainty of the
#         final result (.ge.0.)
#
#  output..
#
#  zeroin abscissa approximating a zero of  f  in the interval ax,bx

ax <- interval[1]
bx <- interval[2]
fa <- f(ax, ...)
fb <- f(bx, ...)
if (trace) cat("Start zeroin: f(",ax,")=",fa,"   f(",bx,")=",fb,"\n")

a <- ax
b <- bx
c <- a
fc <- fa;
wtol <- tol # to avoid changing tol
# First test if we have found a root at an endpoint
if(fa == 0.0) {
wtol <- 0.0
maxit <- 0
if (trace) cat("fa is 0\n")
return(list(root=a, froot=0, wtol=0.0, iter=2))
}
if(fb ==  0.0) {
wtol <- 0.0
maxit <- 0
if (trace) cat("fb is 0\n")
return(list(root=b, froot=0, wtol=0.0, iter=2))
}
maxit <- maxiter + 2 # count evaluations as maxiter-maxit

while(maxit > 0)		## Main iteration loop	*/
{
if (trace) cat("Top of iteration, maxit=",maxit,"\n")

prev_step <- b-a		## Distance from the last but one
##   to the last approximation	*/
#double tol_act;			## Actual tolerance		*/
#double p;			## Interpolation step is calcu- */
#double q;			## lated in the form p/q; divi-
#				 * sion operations is delayed
#				 * until the last moment	*/
#double new_step;		## Step at this iteration	*/

if( abs(fc) < abs(fb) )
{				## Swap data for b to be the	*/
if (trace) cat("fc smaller than fb\n")
a <- b
b <- c
c <- a	## best approximation		*/
fa<-fb
fb<-fc
fc<-fa
}
tol_act <- 2*epsilon*abs(b) + tol/2
if (trace) cat("tol_act= ",tol_act,"\n")

new_step <- (c-b)/2 # bisection
if (trace) cat("new_step= ",new_step,"\n")

if( (abs(new_step) <= tol_act) || (fb == 0.0 ))
{
if (trace) cat("DONE! -- small new_step or fb=0\n")
wtol = abs(c-b)
return(list(root=b, froot=fb, rtol=wtol, iter=maxiter-maxit)) ## Acceptable approx. is found	*/
}

## Decide if the interpolation can be tried	*/
if( (abs(prev_step) >= tol_act)	## If prev_step was large enough*/
&& (abs(fa) > abs(fb)) ) {	## and was in true direction,
## Interpolation may be tried	*/
#    register double t1,cb,t2;
if (trace) cat("prev_step larger than tol_act and fa bigger than fb\n")

cb <- c-b
if( a==c ) {		## If we have only two distinct	*/
## points linear interpolation	*/
t1 <- fb/fa		## can only be applied		*/
p <- cb*t1
q <- 1.0 - t1
if (trace) cat("a == c\n")
}
else {			## Quadric inverse interpolation*/
if (trace) cat("a != c\n")
q <- fa/fc
t1 <- fb/fc
t2 <- fb/fa
p <- t2 * ( cb*q*(q-t1) - (b-a)*(t1-1.0) )
q <- (q-1.0) * (t1-1.0) * (t2-1.0)
}
if( p>0.0 )	{	## p was calculated with the */
if (trace) cat("p>0\n")
q <- -q			## opposite sign; make p positive */
} else {			## and assign possible minus to	*/
if (trace) cat("! p>0\n")
p <- -p			## q				*/
}
if( (p < (0.75*cb*q-abs(tol_act*q)/2))  ## If b+p/q falls in [b,c]*/
&& (p < abs(prev_step*q/2)) ) {	## and isn't too large	*/
if (trace) cat("p satisfies conditions for changing new_step\n")
new_step <- p/q			## it is accepted
}					# * If p/q is too large then the
# * bisection procedure can
# * reduce [b,c] range to more
# * extent */
}

if( abs(new_step) < tol_act) {	## Adjust the step to be not less*/
if (trace) cat("new_step smaller than tol_act\n")
if( new_step > 0.0 )	## than tolerance		*/
new_step <- tol_act
else
new_step <- -tol_act
}
a <- b
fa <- fb			## Save the previous approx. */
b <- b + new_step
fb <- f(b, ...)
maxit <- maxit-1
## Do step to a new approxim. */
if( ((fb > 0) && (fc > 0)) || ((fb < 0) && (fc < 0)) ) {
if (trace) cat("make c to have sign opposite to b\n")
## Adjust c for it to have a sign opposite to that of b */
c <- a
fc <- fa
}

}
## failed! */
if (trace) cat("Failed!\n")
wtol <- abs(c-b)
maxit <- -1
return(list(root=b, froot=NA, rtol=wtol, iter=maxit)) ## Acceptable approx. is found	*/
}
```

## Try the rootoned package in your browser

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

rootoned documentation built on Aug. 29, 2018, 3 a.m.