diff --git a/README.md b/README.md index 3700b0cc..cae4f615 100644 --- a/README.md +++ b/README.md @@ -121,4 +121,4 @@ Some additional documentation can be read [here](http://nbviewer.ipython.org/url ## Polynomials Special methods for finding roots of polynomials have been moved to -the `PolynomialZeros` package and its `polyroots(f, domain)` function. +the `PolynomialZeros` package and its `poly_roots(f, domain)` function. diff --git a/src/Roots.jl b/src/Roots.jl index 71d645f1..948d2827 100644 --- a/src/Roots.jl +++ b/src/Roots.jl @@ -166,7 +166,7 @@ fzeros(p) = Base.depwarn(""" Calling fzeros with just a polynomial is deprecated. Either: * Specify an interval to search over: fzeros(p, a, b). - * Use the `realroots` function from `PolynomialZeros` + * Use the `poly_roots(p, Over.R)` call from `PolynomialZeros` * Use `Polynomials` or `PolynomialRoots` and filter. For example, ``` diff --git a/src/bracketing.jl b/src/bracketing.jl index b1165af5..bb3f5e77 100644 --- a/src/bracketing.jl +++ b/src/bracketing.jl @@ -52,6 +52,64 @@ function _middle(x::Real, y::Real) end +""" + + `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) + + if a > b + b,a = a, b + end + + + 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 + + while a < m < b + fm = sign(f(m)) + + if iszero(fm) || isnan(fm) || isinf(fm) + return m + elseif fa * fm < 0 + b,fb=m,fm + else + a,fa=m,fm + end + m = _middle(a,b) + end + return m +end + + #### ## find_zero interface. We need to specialize for T<:Float64, and BigSomething const BigSomething = Union{BigFloat, BigInt} @@ -64,9 +122,21 @@ type A42 <: AbstractBisection end ## As we don't have the `A42` algorithm implemented through `find_zero`, we adjust here. function find_zero{M<:AbstractBisection, T<:Real}(f, x0::Vector{T}, method::M; maxevals::Int=50, verbose::Bool=false, kwargs...) x = sort(float(x0)) + + if isinf(x[1]) + x[1] = nextfloat(x[1]) + end + if isinf(x[2]) + x[2] = prevfloat(x[2]) + end + if eltype(x) <: Float64 - prob, options = derivative_free_setup(method, DerivativeFree(f), x; verbose=verbose, maxevals=maxevals, kwargs...) - find_zero(prob, method, options) + if verbose + prob, options = derivative_free_setup(method, DerivativeFree(f, f(x0[1])), x; verbose=verbose, maxevals=maxevals, kwargs...) + find_zero(prob, method, options) + else # avoid overhead of generic calling method + bisection64(f, x0[1], x0[2])::eltype(x) + end else a42(f, x[1], x[2]; maxeval=maxevals, verbose=verbose) end @@ -91,11 +161,12 @@ end ## The tolerances can be set to 0, in which case, the termination occurs when `nextfloat(x0) = x2`. ## The bracket `[a,b]` must be bounded. -function init_state{T <: Float64}(method::Bisection, fs, x::Vector{T}, bracket) +function init_state{T <: Float64,R}(method::Bisection, fs::DerivativeFree{R}, x::Vector{T}, bracket) x0, x2 = sort(x[1:2]) isinf(x0) && (x0 = nextfloat(x0)) isinf(x2) && (x2 = prevfloat(x2)) - @compat y0, y2 = fs.f.([x0, x2]) + y0::R = fs.f(x0) + y2::R = fs.f(x2) sign(y0) * sign(y2) > 0 && throw(ArgumentError(bracketing_error)) @@ -154,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: @@ -415,6 +488,7 @@ end """ + Searches for zeros of `f` in an interval [a, b]. Basic algorithm used: @@ -427,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) @@ -441,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 @@ -456,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 diff --git a/src/find_zero.jl b/src/find_zero.jl index 0c2ef36b..d3f59c89 100644 --- a/src/find_zero.jl +++ b/src/find_zero.jl @@ -13,29 +13,32 @@ @compat abstract type UnivariateZeroMethod end # container for callable objects; not really necessary, but has some value. -@compat abstract type CallableFunction end -immutable DerivativeFree <: CallableFunction +@compat abstract type CallableFunction{R} end +immutable DerivativeFree{R} <: CallableFunction{R} f + fx0::R end -immutable FirstDerivative <: CallableFunction +immutable FirstDerivative{R} <: CallableFunction{R} f fp + fx0::R end -immutable SecondDerivative <: CallableFunction +immutable SecondDerivative{R} <: CallableFunction{R} f fp fpp + fx0::R end ## allows override for automatic derivatives, see Newton -function callable_function(m::UnivariateZeroMethod, f) - !isa(f, Tuple) && return DerivativeFree(f) - length(f) == 1 && return DerivativeFree(f[1]) - length(f) == 2 && return FirstDerivative(f[1], f[2]) - SecondDerivative(f[1], f[2], f[3]) +function callable_function(m::UnivariateZeroMethod, f, x0) + !isa(f, Tuple) && return DerivativeFree(f, f(x0)) + length(f) == 1 && return DerivativeFree(f[1], f(x0)) + length(f) == 2 && return FirstDerivative(f[1], f[2], f(x0)) + SecondDerivative(f[1], f[2], f[3], f(x0)) end @@ -146,14 +149,15 @@ function derivative_free_setup{T<:AbstractFloat}(method::Any, fs::CallableFuncti ) x = float(x0) prob = UnivariateZeroProblem(fs, x, isa(bracket, Nullable) ? bracket : Nullable(convert(Vector{T}, bracket))) - options = univariate_zero_options(;xabstol=xabstol, - xreltol=xreltol, - abstol=abstol, - reltol=reltol, - maxevals=maxevals, - maxfnevals=maxfnevals, - verbose=verbose, - kwargs...) + # options = univariate_zero_options(;xabstol=xabstol, + # xreltol=xreltol, + # abstol=abstol, + # reltol=reltol, + # maxevals=maxevals, + # maxfnevals=maxfnevals, + # verbose=verbose, + # kwargs...) + options = UnivariateZeroOptions(xabstol, xreltol, abstol, reltol, maxevals, maxfnevals, verbose) prob, options end @@ -399,7 +403,7 @@ find_zero(f, 1.0, Steffensen()) # also Order2() ``` """ find_zero{T<:Number}(f, x0::Union{T,Vector{T}}, method::UnivariateZeroMethod; kwargs...) = - find_zero(method, callable_function(method, f), x0; kwargs...) + find_zero(method, callable_function(method, f, float(x0)), x0; kwargs...) ## some defaults for methods find_zero{T <: Number}(f, x0::T; kwargs...) = find_zero(f, x0, Order0(); kwargs...) @@ -418,7 +422,7 @@ end ## old interface for fzero ## old keyword arguments (see ?fzero) handled in univariate_zero_options -function derivative_free{T <: AbstractFloat}(f, x0::T; order::Int=0, +@noinline function derivative_free{T <: AbstractFloat}(f, x0::T; order::Int=0, kwargs...) if order == 0 diff --git a/src/newton.jl b/src/newton.jl index 092c6736..a339987b 100644 --- a/src/newton.jl +++ b/src/newton.jl @@ -7,11 +7,11 @@ type Newton <: UnivariateZeroMethod end -function callable_function(method::Newton, f::Tuple) +function callable_function(method::Newton, f::Tuple, x0) length(f) == 1 && return FirstDerivative(f[1], D(f[1])) - FirstDerivative(f[1], f[2]) + FirstDerivative(f[1], f[2], f[1](x0)) end -callable_function(method::Newton, f::Any) = FirstDerivative(f, D(f)) +callable_function(method::Newton, f::Any, x0) = FirstDerivative(f, D(f), x0) function init_state{T}(method::Newton, fs, x0::T) state = UnivariateZeroState(x0, @@ -143,12 +143,12 @@ newton(f, fp, x0; kwargs...) = find_zero((f, fp), x0, Newton(); kwargs...) type Halley <: UnivariateZeroMethod end -function callable_function(method::Halley, f::Tuple) - length(f) == 1 && return SecondDerivative(f[1], D(f[1]), D(f[1],2)) - length(f) == 2 && return SecondDerivative(f[1], f[2], D(f[2],1)) - SecondDerivative(f[1], f[2], f[3]) +function callable_function(method::Halley, f::Tuple, x0) + length(f) == 1 && return SecondDerivative(f[1], D(f[1]), D(f[1],2), f[1](x0)) + length(f) == 2 && return SecondDerivative(f[1], f[2], D(f[2],1), f[1](x0)) + SecondDerivative(f[1], f[2], f[3], f[1](x0)) end -callable_function(method::Halley, f) = SecondDerivative(f, D(f), D(f, 2)) +callable_function(method::Halley, f, x0) = SecondDerivative(f, D(f), D(f, 2), f(x0)) function update_state{T}(method::Halley, fs, o::UnivariateZeroState{T}, options::UnivariateZeroOptions) diff --git a/test/test_fzero.jl b/test/test_fzero.jl index b546dc9d..989bd125 100644 --- a/test/test_fzero.jl +++ b/test/test_fzero.jl @@ -19,7 +19,7 @@ end fn, xstar, x0, br = x -> x^5 - x - 1, 1.1673039782614187, 1.0, [1.0, 2.0] @test fzero(fn, x0, order=1) ≈ xstar -@test_throws Roots.ConvergenceFailed fzero(fn, x0, order=1, maxeval=2) +@test_throws Roots.ConvergenceFailed fzero(fn, x0, order=1, maxevals=2) #@test norm(fzero(fn, x0, order=1, ftol=1e-2) - xstar) > 1e-5 #@test norm(fzero(fn, x0, order=1, xtol=1e-2) - xstar) > 1e-10