Skip to content

Commit

Permalink
adjust docs, adjust tolerance names on unexported function
Browse files Browse the repository at this point in the history
  • Loading branch information
jverzani committed Jun 22, 2017
1 parent 479ed2a commit b2a6899
Showing 1 changed file with 39 additions and 8 deletions.
47 changes: 39 additions & 8 deletions src/bracketing.jl
Original file line number Diff line number Diff line change
Expand Up @@ -52,9 +52,33 @@ function _middle(x::Real, y::Real)
end


## Unexported but much speedier version of bisection method for float64
## guaranteed to terminate
function bisection64(f, a::Real, b::Real)
"""
`Roots.bisection64(f, a, b)` (unexported)
* `f`: a callable object, like a function
* `a`, `b`: Real values specifying a *bracketing* interval (one with
`f(a) * f(b) < 0`). These will be converted to `Float64` values.
Runs the bisection method using midpoints determined by a trick
leveraging 64-bit floating point numbers. After ensuring the
intermediate bracketing interval does not straddle 0, the "midpoint"
is half way between the two values onces converted to unsigned 64-bit
integers. This means no more than 64 steps will be taken, fewer if `a`
and `b` already share some bits.
The process is guaranteed to return a value `c` with `f(c)` one of
`0`, `Inf`, or `NaN`; *or* one of `f(prevfloat(c))*f(c) < 0` or
`f(c)*f(nextfloat(c)) > 0` holding.
This function is a bit faster than the slightly more general
`find_zero(f, [a,b], Bisection())` call.
Due to Jason Merrill.
"""
function bisection64(f, a::Number, b::Number)

a,b = Float64(a), Float64(b)

Expand All @@ -65,6 +89,8 @@ function bisection64(f, a::Real, b::Real)

m = _middle(a,b)
fa, fb = sign(f(a)), sign(f(b))

fa * fb > 0 && throw(ArgumentError(bracketing_error))
iszero(fa) || isnan(fa) || isinf(fa) && return a
iszero(fb) || isnan(fb) || isinf(fb) && return b

Expand Down Expand Up @@ -199,6 +225,8 @@ end

"""
`Roots.a42(f, a, b; kwargs...)` (not exported)
Finds the root of a continuous function within a provided
interval [a, b], without requiring derivatives. It is based on algorithm 4.2
described in: 1. G. E. Alefeld, F. A. Potra, and Y. Shi, "Algorithm 748:
Expand Down Expand Up @@ -460,6 +488,7 @@ end


"""
Searches for zeros of `f` in an interval [a, b].
Basic algorithm used:
Expand All @@ -472,10 +501,12 @@ If there are many zeros relative to the number of points, the process
is repeated with more points, in hopes of finding more zeros for
oscillating functions.
Called by `fzeros` or `Roots.find_zeros`.
"""
function find_zeros(f, a::Real, b::Real, args...;
no_pts::Int=100,
ftol::Real=10*eps(), reltol::Real=10*eps(),
abstol::Real=10*eps(), reltol::Real=10*eps(), ## should be abstol, reltol as used.
kwargs...)

a, b = a < b ? (a,b) : (b,a)
Expand All @@ -486,13 +517,13 @@ function find_zeros(f, a::Real, b::Real, args...;
## Look in [ai, bi)
for i in 1:(no_pts+1)
ai,bi=xs[i:i+1]
if isapprox(f(ai), 0.0, rtol=reltol, atol=ftol)
if isapprox(f(ai), 0.0, rtol=reltol, atol=abstol)
push!(rts, ai)
elseif sign(f(ai)) * sign(f(bi)) < 0
push!(rts, find_zero(f, [ai, bi], Bisection()))
else
try
x = find_zero(f, ai + (0.5)* (bi-ai), Order8(); maxevals=10, reltol=ftol, xreltol=reltol)
x = find_zero(f, ai + (0.5)* (bi-ai), Order8(); maxevals=10, abstol=abstol, reltol=reltol)
if ai < x < bi
push!(rts, x)
end
Expand All @@ -501,11 +532,11 @@ function find_zeros(f, a::Real, b::Real, args...;
end
end
## finally, b?
isapprox(f(b), 0.0, rtol=reltol, atol=ftol) && push!(rts, b)
isapprox(f(b), 0.0, rtol=reltol, atol=abstol) && push!(rts, b)

## redo if it appears function oscillates alot in this interval...
if length(rts) > (1/4) * no_pts
return(find_zeros(f, a, b, args...; no_pts = 10*no_pts, ftol=ftol, reltol=reltol, kwargs...))
return(find_zeros(f, a, b, args...; no_pts = 10*no_pts, abstol=abstol, reltol=reltol, kwargs...))
else
return(sort(rts))
end
Expand Down

0 comments on commit b2a6899

Please sign in to comment.