The following function is one of the most difficult test functions for
univariate root finding I know of,
julia> function fn(x::Real)
return log(x) + x^2/(2*exp(1)) - 2 * x/sqrt(exp(1)) + 1
end
whose exact root between 1.0 and 3.4 is sqrt(e) = 1.6487212707001282... .
Applying your nice version of bisect_root to this function returns
julia> bisect_root(fn, 1.0, 3.4)
(1.648707085670482,1.6487070856704822)
with an absolute error of about 1.5e-05 ! This innocent looking test
function appears to be so flat that only an arbitrary precision solver will
locate this root more exactly.
Hans Werner
BTW It would be nice to have Ridders' algorithm available in Julia, too,
about which the "Numerical Recipes" say:
"*In both reliability and speed, Ridders's method is generally competitive
with the more highly developed and better established (but more
complicated) method of Van Wijngaarden, Dekker, and Brent [...].*"
My Julia version of Ridders' method would be as follows, perhaps it can be
improved with your proposed x1 < (x1 + x2)/2 < x2 :
#-- --------------------------------------------------------------------
function ridders(f::Function, a::Real, b::Real;
maxiter::Integer = 1000, tol::Real = eps())
x1 = a; x2 = b
f1 = f(x1); f2 = f(x2)
if f1 * f2 > 0
error("f(a) and f(b) must have different signs.")
elseif f1 == 0
return x1
elseif f2 == 0
return f2
end
niter = 2
while niter < maxiter && abs(x1 - x2) > tol
xm = (x1 + x2)/2.0
fm = f(xm)
if fm == 0; return xm; end
x3 = xm + (xm - x1) * sign(f1 - f2) * fm / sqrt(fm^2 - f1 * f2)
f3 = f(x3)
niter += 2
if f3 == 0; return x3; end
if fm * f3 < 0
x1 = xm; f1 = fm
x2 = x3; f2 = f3
elseif f1 * f3 < 0
x2 = x3; f2 = f3
elseif f2 * f3 < 0
x1 = x3; f1 = f3
else
error("Inform the maintainer: you should never get here.")
end
end
return (x1 + x2) / 2.0
end
#-- --------------------------------------------------------------------
