Rigorous numerics for maxima and minima (one variable) Let $f:\mathbb{R}_0^+\to \mathbb{R}$ be defined by some combination of the four basic operations and square roots. (The argument of square-roots is assumed is to be non-negative, and the value of square roots is defined to be non-negative as well.) How can one find the global maximum and minimum of $f$ numerically, with full rigor?
(A rigorous numerical solution has to be given as an interval, e.g., [1.500001,1.500002],
within which the maximum (say) lies.)
Rather convincing but not quite rigorous method:
Plot $f'(x)$, see roughly where the zeroes are, and narrow them down using some standard
(e.g. SAGE, mathematica). Then verify this (non-rigorous) datum rigorously simply by
checking that $f'(x)$ has different signs at $x=x_0-\epsilon$ and at $x=x_0+\epsilon$,
where $x_0$ is an alleged zero of $f'(x)$. Compare the values of $f(x)$ at all $x_0$, and also at $x=0$ and as $x\to \infty$.
There is really only one way in which this procedure isn't rigorous: there is no guarantee that there aren't any zeroes of $f'(x)$ we have missed. (If $f$ is a polynomial, this can be easily dealt with, but $f$ is not necessarily a polynomial.)
Are there any (free) programs out there that for the minimum and maximum of $f$ rigorously?
 A: There is an interval algorithm for finding global extrema. See, e.g. section 5.5 of Jaulin et al. or 5.2 of Tucker for overviews, with some C++ source code.
A: You can convert your problem into a multivariate one, with polynomial equality constraints, by introducing new variables for each square root and division. Thus you end up with the problem of minimizing a polynomial function on a semialgebraic set (add various necessary conditions for the optimum, too, e.g. what you get from $f'(x)=0$; this won't need more variables to be introduced), which you indeed can do rigorously, i.e. you could do even more: find your minimum and maximum as real algebraic numbers, by quantifier elimination. 
If you're lucky then necessary conditions for the optimum give you a zero-dimensional ideal to deal with, with zeros being (potential) optima.
Well, this is in theory; in practice, you might run into problems, and instead opt for using a sum of squares-based approach of polynomial optimization, which turns into solving a sequence of semidefinite programming problems of increasing size, with optimal values approximating the real optimum with increasing precision.
