R help - Jan 2007 - what is the purpose of an error message in uniroot?

Hi all,

This is probably a blindingly obvious question: Why does it matter in
the uniroot function whether the f() values at the end points that you
supply are of the same sign?

For example:

f <- function(x,y) {y-x^2+1}

#this gives a warning
uniroot(f,interval=c(-5,5),y=0)
Error in uniroot(f, interval=c(-5, 5), y = 0) : f() values at end
points not of opposite sign

#this doesn't give a warning
uniroot(f,interval=c(.1,5),y=0)
$root
[1] 1

$f.root
[1] 1.054e-05

$iter
[1] 9

$estim.prec
[1] 6.104e-05

If I comment out the two lines of script in the uniroot function that
produce this warning and create a new function, call it uniroot2,
everything works as I'd like. But for didactic purposes, why did the
creators of uniroot want the f() values at endpoints to be of opposite
sign?

Thanks in advance,

Matt

mckellercran at gmail.com wrote:
> This is probably a blindingly obvious question:
        Yes, it is.
> Why does it matter in the uniroot function whether the f() values at
> the end points that you supply are of the same sign?
        Plot some graphs.

        Think about the *name* of the function --- *uni*root.

        Does that ring any bells?

        And how do you know there *is* a root in the interval
        in question?  Try your ``uniroot2'' on f(x) = 1+x^2
        and the interval [-5,5].

        To belabour the point --- if the f() values are of the
        same sign, then there are 0, or 2, or 4, or ....
        roots in the interval in question.

        The ***only chance*** you have of there being a unique
        root is if the f() values are of opposite sign.

        The algorithm used and the precision estimates returned
        presumably depend on the change of sign.  You can get
        answers --- sometimes --- if the change of sign is not
        present, but the results could be seriously misleading.

        Without the opposite sign requirement the user will often
        wind up trying to do something impossible or getting
        results about which he/she is deluded.

                                cheers,

                                        Rolf Turner
                                        rolf at math.unb.ca

P. S.  If the f() values are of the same sign, uniroot() DOES
NOT give a warning!  It gives an error.

                                        R. T.

Matt,

Some time back I didn't like the uniroot restriction either so I wrote a
short function "manyroots" that breaks an interval into many shorter
intervals and looks for a single root in each of them.  This function is NOT
guaranteed to find all roots in an interval even if you specify many
subintervals.  Graphing is always a good idea.  There is always a chance the
function dips to or below the axis and back up in an arbitrarily short interval.
For example, f(x) = x^2.  If one of your subintervals doesn't happen to end
at 0, you're out of luck.  (This is in contrast to uniroot working on a
continuous function that is positive at one end of an interval and negative at
the other end:  at least one root is guaranteed and uniroot will find it given
enough iterations.)  Furthermore, if you are working with a polynomial, just use
polyroot.

Chris

(Sorry for the dearth of code comments in the following)

manyroots <- function(f, interval, ints=1, maxlen=NULL,
  lower = min(interval), upper = max(interval),
  tol = .Machine$double.eps^0.25, maxiter = 1000, ...)
{
    if (!is.numeric(lower) || !is.numeric(upper) || lower >        upper)
        stop("lower < upper  is not fulfilled")
    if (is.infinite(lower) || is.infinite(upper))
        stop("Interval must have finite length")
    if (!is.null(maxlen))
        ints <- ceiling((upper-lower)/maxlen)
    if (!is.numeric(ints) || length(ints)>1 || floor(ints)!=ints ||
ints<1)
        stop("ints must be positive integer")

    ends <- seq(lower, upper, length=ints+1)
    fends <- numeric(length(ends))
    for (i in seq(along=ends)) fends[i] <- f(ends[i], ...)

    zeros <- iters <- prec <- rep(NA, ints)

    for (i in seq(ints)) {
        cat(i, ends[i], ends[i+1], fends[i], fends[i+1], "\n")
        if (fends[i] * fends[i+1] > 0) {
#            cat("f() values at end points not of opposite sign\n")
            next;
        }
        if (fends[i] == 0 & i>1) {
#            cat("this was found in previous iteration\n")
            next;
        }
        
        val <- .Internal(zeroin(function(arg) f(arg, ...), ends[i],
            ends[i+1], tol, as.integer(maxiter)))
        if (as.integer(val[2]) == maxiter) {
            warning("Iteration limit (", maxiter, ") reached in
interval (",
              ends[i], ",", ends[i+1], ").")
        }
        zeros[i] <- val[1]
        iters[i] <- val[2]
        prec[i] <- val[3]
    }
    zeros <- as.vector(na.omit(zeros))
    fzeros <- numeric(length(zeros))
    for (i in seq(along=zeros)) fzeros[i] <- f(zeros[i], ...)

    list(root = zeros, f.root = fzeros,
        iter = as.vector(na.omit(iters)),
        estim.prec = as.vector(na.omit(prec)))
}

gg <- function(x) x*(x-1)*(x+1)
manyroots(gg, c(-4,4), 13, maxiter=200, tol=10^-10)

hh <- function(x,x2) x^2-x2
manyroots(hh, c(-10, 10), maxlen=.178, x2=9)

manyroots(sin, c(-4,20), maxlen=.01)
#but
ss <- function(x) sin(x)^2
manyroots(ss, c(-4,20), maxlen=.01)
plot(ss, -4,20)
abline(h=0)



From: rolf at math.unb.ca

mckellercran at gmail.com wrote:
> > This is probably a blindingly obvious question:
        Yes, it is.
> > Why does it matter in the uniroot function whether the f() values at
> > the end points that you supply are of the same sign?
        Plot some graphs.

        Think about the *name* of the function --- *uni*root.

        Does that ring any bells?

        And how do you know there *is* a root in the interval
        in question?  Try your ``uniroot2'' on f(x) = 1+x2
        and the interval [-5,5].

        To belabour the point --- if the f() values are of the
        same sign, then there are 0, or 2, or 4, or ....
        roots in the interval in question.

Rolf,
Only if f is continuous.... (of course finding roots of discontinuous functions
is a greater challenge)

        The ***only chance*** you have of there being a unique
        root is if the f() values are of opposite sign.

        The algorithm used and the precision estimates returned
        presumably depend on the change of sign.  You can get
        answers --- sometimes --- if the change of sign is not
        present, but the results could be seriously misleading.

        Without the opposite sign requirement the user will often
        wind up trying to do something impossible or getting
        results about which he/she is deluded.

                                cheers,

                                        Rolf Turner
                                        rolf at math.unb.ca

P. S.  If the f() values are of the same sign, uniroot() DOES
NOT give a warning!  It gives an error.

                                        R. T.




-- 
Christopher Andrews, PhD
SUNY Buffalo, Department of Biostatistics
242 Farber Hall, candrews at buffalo.edu, 716 829 2756