diff --git a/.github/workflows/downstream.yml b/.github/workflows/downstream.yml new file mode 100644 index 00000000..ce2c25d1 --- /dev/null +++ b/.github/workflows/downstream.yml @@ -0,0 +1,49 @@ +name: IntegrationTest +on: + push: + branches: [master] + tags: [v*] + pull_request: + +jobs: + test: + name: ${{ matrix.package.repo }}/${{ matrix.package.group }}/${{ matrix.julia-version }} + runs-on: ${{ matrix.os }} + env: + GROUP: ${{ matrix.package.group }} + strategy: + fail-fast: false + matrix: + julia-version: [1] + os: [ubuntu-latest] + package: + - {user: lrnv, repo: Copulas.jl, group: All} + steps: + - uses: actions/checkout@v6 + - uses: julia-actions/setup-julia@v2 + with: + version: ${{ matrix.julia-version }} + arch: x64 + - uses: julia-actions/julia-buildpkg@latest + - name: Clone Downstream + uses: actions/checkout@v6 + with: + repository: ${{ matrix.package.user }}/${{ matrix.package.repo }} + path: downstream + - name: Load this and run the downstream tests + shell: julia --color=yes --project=downstream {0} + run: | + using Pkg + try + # force it to use this PR's version of the package + Pkg.develop(PackageSpec(path=".")) # resolver may fail with main deps + Pkg.update() + Pkg.test() # resolver may fail with test time deps + catch err + err isa Pkg.Resolve.ResolverError || rethrow() + # If we can't resolve that means this is incompatible by SemVer and this is fine + # It means we marked this as a breaking change, so we don't need to worry about + # Mistakenly introducing a breaking change, as we have intentionally made one + @info "Not compatible with this release. No problem." exception=err + exit(0) # Exit immediately, as a success + end \ No newline at end of file diff --git a/Project.toml b/Project.toml index 733fb085..cc3541a6 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "Roots" uuid = "f2b01f46-fcfa-551c-844a-d8ac1e96c665" -version = "3.0.0" +version = "3.0.6" [deps] Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697" @@ -32,7 +32,7 @@ ChainRulesCore = "1" ChainRulesTestUtils = "1" CommonSolve = "0.1, 0.2" ForwardDiff = "0.10, 1" -IntervalRootFinding = "0.5, 0.6" +IntervalRootFinding = "0.5, 0.6, 0.7" JET = "0.9, 0.10, 0.11" JSON = "0.21, 1" Measurements = "2.11" diff --git a/docs/src/index.md b/docs/src/index.md index 5374614c..3edc2e28 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -6,7 +6,7 @@ Documentation for [Roots.jl](https://github.com/JuliaMath/Roots.jl) ## About `Roots` is a `Julia` package for finding zeros of continuous -scalar functions of a single real variable using floating point numbers. That is solving ``f(x)=0`` for ``x`` adjusting for floating-point idiosyncrasies. +scalar functions of a single real variable using floating point numbers. That is, solving ``f(x)=0`` for ``x`` adjusting for floating-point idiosyncrasies. The `find_zero` function provides the primary interface. It supports various algorithms through the @@ -21,9 +21,9 @@ specification of a method. These include: computations allows. Other methods include `A42`, `AlefeldPotraShi`, `Roots.Brent`, `Roots.Chandrapatlu`, `Roots.ITP`, `Roots.Ridders`, `Roots.ModAB`, and ``12``-flavors of - `FalsePosition`. The default bracketing method for + `FalsePosition` and a few in `RegulaFalsi`. The default bracketing method for the basic floating-point types is `Bisection` , as it is more robust to some inputs, - but `A42` and `AlefeldPotraShi` typically converge in a few + but `Roots.ModAB`, `A42` and `AlefeldPotraShi` typically converge in a few iterations and are more performant. @@ -32,7 +32,7 @@ specification of a method. These include: (the Steffensen method), `Order5`, `Order8`, and `Order16`. The number indicates, roughly, the order of convergence. The `Order0` method is the default, and the most robust, as it finishes off with - a bracketing method when a bracket is encountered, The higher order + a bracketing method when a bracket is encountered, The higher-order methods promise higher order (faster) convergence, though don't always yield results with fewer function calls than `Order1` or `Order2`. `Roots.Sidi` is a family of methods. diff --git a/docs/src/reference.md b/docs/src/reference.md index 7e113103..d28ef1ee 100644 --- a/docs/src/reference.md +++ b/docs/src/reference.md @@ -142,6 +142,7 @@ Roots.Chandrapatla Roots.Ridders Roots.ITP Roots.ModAB +Roots.RegulaFalsi FalsePosition Roots.LithBoonkkampIJzermanBracket Roots.BracketedHalley @@ -227,9 +228,12 @@ good example. | Bracketing | A42 | ``(2 + 7^{1/2})`` | ``3,4`` |``(2 + 7^{1/2})^{1/3} = 1.6686\dots`` | | Bracketing | AlefeldPotraShi | | ``3,4`` | ``1.618\dots`` | | Bracketing | Brent | ``\leq 1.89\dots`` | ``1`` | ``\leq 1.89\dots`` | -| Bracketing | ITP | ``\leq \varphi`` | ``1`` | ``\leq \varphi`` | +| Bracketing | ITP | ``\leq \varphi`` | ``1`` | ``\leq \varphi`` | | Bracketing | Ridders | ``1.83\dots`` | ``2`` | ``1.225\dots`` | -| Bracketing | FalsePosition | ``1.442\dots`` | ``1`` | ``1.442\dots`` | +| Bracketing | RegulaFalsi{:classic} | ``1`` | ``1`` | ``1`` | +| Bracketing | RegulaFalsi{:Illinois} | ``1.442\dots`` | ``1`` | ``1.442\dots`` | +| Bracketing | RegulaFalsi{:AndersonBjork} | ``1.681\dots`` | ``1`` | ``1.681\dots`` | +| Bracketing | RegulaFalsi{:Ford4} | ``1.681\dots`` | ``1`` | ``1.681\dots`` | | Bracketing | ModAB | ``≈1.7\dots`` | ``1`` | ``1.7\dots`` | | Bracketing | LithBoonkkampIJzermanBracket | ``2.91`` | ``3`` | ``1.427\dots`` | | Robust | King | ``\varphi=1.618\dots`` | ``2`` | ``1.272\dots`` | @@ -247,7 +251,7 @@ Identifying when an algorithm converges or diverges requires specifications of t In the case of exact bisection, convergence is mathematically guaranteed. For floating point numbers, either an *exact* zero is -found, or the bracketing interval can be subdivided into ``[a_n,b_n]`` +found, or the bracketing interval can *not* be subdivided into ``[a_n,b_n]`` with ``a_n`` and ``b_n`` being adjacent floating point values. That is ``b_n-a_n`` is as small as possible in floating point numbers. This can be considered a stopping criteria in ``\Delta x``. For early termination @@ -263,7 +267,7 @@ determine closeness. Relying on the closeness of two ``x`` values will not be adequate for all problems, as there are examples where the difference ``\Delta_n=|x_n-x_{n-1}|`` can be quite small, ``0`` even, yet ``f(x_n)`` is -not near a ``0``. As such, for non-bracketing methods, a check on the +not near a ``0``. As well, a final step which should make $\Delta_n$ small enough, might get derailed due to floating point issues. As such, for non-bracketing methods and for some non-strict bracketing methods, a check on the size of ``f(x_n)`` is also used. As we find floating point approximations to ``\alpha``, the zero, we must consider values small when ``f(\alpha(1+\epsilon))`` is small. By Taylor's approximation, we @@ -311,6 +315,7 @@ Roots.bisection Roots.muller Roots.newton Roots.dfree +Roots.modab ``` diff --git a/docs/src/roots.md b/docs/src/roots.md index af2de905..90e8603a 100644 --- a/docs/src/roots.md +++ b/docs/src/roots.md @@ -18,7 +18,7 @@ julia> using Roots, ForwardDiff Consider the polynomial function ``f(x) = x^5 - x + 1/2``. As a polynomial, its roots, or zeros, could be identified with the `roots` function of the `Polynomials` package. However, even that function uses a numeric method to identify the values, as no solution with radicals is available. That is, even for polynomials, non-linear root finders are needed to solve ``f(x)=0``. (Though polynomial root-finders can exploit certain properties not available for general non-linear functions.) -The `Roots` package provides a variety of algorithms for this task. In this quick overview, only the default ones are illustrated. +The `Roots` package provides a variety of algorithms for this task. In this quick overview, only the default ones are emphasized. For the function ``f(x) = x^5 - x + 1/2`` a simple plot over ``[-2,2]``will show a zero somewhere **between** ``-1.5`` and ``-0.5`` and two zeros near ``0.6``. ("Between", as the continuous function has different signs at ``-1.5`` and ``-0.5``.) @@ -48,9 +48,9 @@ true ``` -This finds the answer to the left of the starting point. To get the other nearby zero, a starting point closer to the answer can be used. +This finds the answer to the left of the starting point. To get the other nearby zero, a starting point closer to that answer can be used. -However, an initial graph might convince one that any of the up-to-``5`` real roots will occur between ``-2`` and ``2``. The `find_zeros` function uses heuristics and a few of the algorithms to identify all zeros between the specified range. Here the method successfully identifies all ``3``: +However, an initial graph might convince one that any of the up-to-``5`` real roots will occur between ``-2`` and ``2``. The `find_zeros` function uses heuristics and a few of the algorithms to attempt to identify all zeros between the specified range. Here the method successfully identifies all ``3``: ``` julia> find_zeros(f, -2, 2) @@ -70,10 +70,10 @@ different signs at $a$ and $b$. If $f$ is a continuous function this ensures ([Bolzano](https://en.wikipedia.org/wiki/Intermediate_value_theorem)) there will be a zero in the interval $[a,b]$. If $f$ is not -continuous, then there must be a point $c$ in $[a,b]$ where the function +continuous, then there must be at least one value $c$ in $[a,b]$ where the function "jumps" over $0$. -Such values can be found, up to floating point +Such a value can be found, up to floating point round off. That is, given `f(a) * f(b) < 0`, a value `c` with `a < c < b` can be found where either `f(c) == 0.0` or `f(prevfloat(c)) * f(c) < 0` or `f(c) * f(nextfloat(c)) < 0`. @@ -89,17 +89,17 @@ condition and `Bisection()` to specify the algorithm: ```jldoctest roots julia> f(x) = cos(x) - x; -julia> x = find_zero(f, (0, pi/2), Bisection()) +julia> c = find_zero(f, (0, pi/2), Bisection()) 0.7390851332151607 -julia> x, f(x) +julia> c, f(c) (0.7390851332151607, 0.0) ``` -For this function we see that `f(x)` is `0.0`. +For this function we see that `f(c)` is `0.0`. -Functions may be parameterized. The following is a similar function as above, still having ``(0, \pi/2)`` as a bracket when ``p>0``. By passing in values of `p` to `find_zero`, different, related problems may be solved. +Functions may be parameterized. The following uses a similar function as above, still having ``(0, \pi/2)`` as a bracket when ``p>0``. By passing in values of `p` to `find_zero`, different, related problems may be solved. ```jldoctest roots julia> g(x, p=1) = cos(x) - x/p; @@ -126,33 +126,33 @@ pattern for `find_zero` is `find_zero(f, x0, M; kwargs...)`, where `kwargs` can specify details about parameters for the problem or tolerances for the solver. In this call `Bisection()` is not specified, as it will be the default (as the initial value is -not specified as a number is over `Float64` values: +not specified as a number and is over `Float64` values): ```jldoctest roots julia> f(x) = sin(x); -julia> x = find_zero(f, (pi/2, 3pi/2)) +julia> c = find_zero(f, (pi/2, 3pi/2)) 3.141592653589793 -julia> x, f(x) +julia> c, f(c) (3.141592653589793, 1.2246467991473532e-16) ``` -This value of `x` does not exactly produce a zero, however, it is as close as can be: +This value of `c` does not exactly produce a zero, however, it is as close as can be: ```jldoctest roots -julia> f(prevfloat(x)) * f(x) < 0.0 || f(x) * f(nextfloat(x)) < 0.0 +julia> f(prevfloat(c)) * f(c) < 0.0 || f(c) * f(nextfloat(c)) < 0.0 true ``` -That is, at `x` the function is changing sign. +That is, at `c` the function is changing sign. From a mathematical perspective, a zero is guaranteed for a *continuous* function. However, the computer algorithm doesn't assume continuity, it just looks for changes of sign. As such, the algorithm -will identify discontinuities, not just zeros. For example: +will identify discontinuities, not just zeros. For example: ```jldoctest roots julia> find_zero(x -> 1/x, (-1, 1)) @@ -181,8 +181,8 @@ iterations. For the problem of finding a zero of `sin` in the interval `(big(3), big(4))`, the default bisection takes ``252`` iterations, whereas the `A42` method takes ``4``. -The algorithms of Alefeld, Potra, and Shi and the well known algorithm -of Brent, also start with a bracketing algorithm. For many problems +The algorithms of Alefeld, Potra, and Shi; the `Roots.ModAB` algorithm; the well-known algorithm +of Brent; and others, also start with a bracketing algorithm. For many problems these will take far fewer steps than the bisection algorithm to reach convergence. These may be called directly. For example, @@ -195,11 +195,11 @@ julia> find_zero(sin, (3,4), A42()) By default, bisection will converge to machine tolerance. This may provide more accuracy than desired. A tolerance may be specified to terminate early, thereby utilizing fewer resources. For example, the following -``4`` steps to reach accuracy to $1/16$ (without specifying `xatol` it uses +uses ``4`` steps to reach accuracy to $1/16$ (without specifying `xatol` it uses ``53`` steps): ```jldoctest roots -julia> rt = find_zero(sin, (3.0, 4.0), xatol=1/16) +julia> rt = find_zero(sin, (3.0, 4.0); xatol=1/16) 3.125 julia> rt - pi @@ -231,10 +231,10 @@ the zero is found through: ```jldoctest roots julia> f(x) = cos(x) - x; -julia> x = find_zero(f , 1) +julia> c = find_zero(f, 1) 0.7390851332151607 -julia> x, f(x) +julia> c, f(c) (0.7390851332151607, 0.0) ``` @@ -244,10 +244,10 @@ For the polynomial $f(x) = x^3 - 2x - 5$, an initial guess of $2$ seems reasonab ```jldoctest roots julia> f(x) = x^3 - 2x - 5; -julia> x = find_zero(f, 2) +julia> c = find_zero(f, 2) 2.0945514815423265 -julia> f(x), sign(f(prevfloat(x)) * f(x)), sign(f(x) * f(nextfloat(x))) +julia> f(c), sign(f(prevfloat(c)) * f(c)), sign(f(c) * f(nextfloat(c))) (-8.881784197001252e-16, 1.0, -1.0) ``` @@ -255,10 +255,10 @@ julia> f(x), sign(f(prevfloat(x)) * f(x)), sign(f(x) * f(nextfloat(x))) For even more precision, `BigFloat` numbers can be used ```jldoctest roots -julia> x = find_zero(sin, big(3)) +julia> c = find_zero(sin, big(3)) 3.141592653589793238462643383279502884197169399375105820974944592307816406286198 -julia> x, sin(x), x - pi +julia> c, sin(c), c - pi (3.141592653589793238462643383279502884197169399375105820974944592307816406286198, 1.096917440979352076742130626395698021050758236508687951179005716992142688513354e-77, 0.0) ``` @@ -280,10 +280,10 @@ method after the initial starting point: ```jldoctest roots julia> f(x) = 2x - exp(-x); -julia> x = find_zero(f, 1, Order1()) +julia> c = find_zero(f, 1, Order1()) 0.3517337112491958 -julia> x, f(x) +julia> c, f(c) (0.3517337112491958, -1.1102230246251565e-16) ``` @@ -293,20 +293,20 @@ Similarly, ```jldoctest roots julia> f(x) = (x + 3) * (x - 1)^2; -julia> x = find_zero(f, -2, Order2()) +julia> c = find_zero(f, -2, Order2()) -3.0 -julia> x, f(x) +julia> c, f(c) (-3.0, 0.0) ``` ```jldoctest roots -julia> x = find_zero(f, 2, Order8()) +julia> c = find_zero(f, 2, Order8()) 1.0000000131073141 -julia> x, f(x) +julia> c, f(c) (1.0000000131073141, 6.87206736323862e-16) ``` @@ -318,9 +318,6 @@ longer super-linear. This is the case here, where `Order2` uses $51$ function calls, `Order8` uses $42$, and `Order0` takes $80$. The `Roots.Order2B` method is useful when a multiplicity is expected; on this problem it takes ``17`` function calls. -To investigate an algorithm and its convergence a `Roots.Tracks` object can be used to store the intermediate values. - - For some functions, adjusting the default tolerances may be necessary to achieve convergence. The tolerances include `atol` and `rtol`, which are used to check if $f(x_n) \approx 0$; @@ -328,6 +325,8 @@ used to check if $f(x_n) \approx 0$; number of iterations in the algorithm. +To investigate an algorithm and its convergence a `Roots.Tracks` object can be used to store the intermediate values. + ## Classical methods @@ -342,10 +341,10 @@ julia> f(x) = x^3 - 2x - 5; julia> fp(x) = 3x^2 - 2; -julia> x = Roots.find_zero((f, fp), 2, Roots.Newton()) +julia> c = Roots.find_zero((f, fp), 2, Roots.Newton()) 2.0945514815423265 -julia> x, f(x) +julia> c, f(c) (2.0945514815423265, -8.881784197001252e-16) ``` @@ -358,10 +357,10 @@ quadratic convergence. It uses the second derivative as well: ```jldoctest roots julia> fpp(x) = 6x; -julia> x = Roots.find_zero((f, fp, fpp), 2, Roots.Halley()) +julia> c = Roots.find_zero((f, fp, fpp), 2, Roots.Halley()) 2.0945514815423265 -julia> x, f(x), sign(f(prevfloat(x)) * f(nextfloat(x))) +julia> c, f(c), sign(f(prevfloat(c)) * f(nextfloat(c))) (2.0945514815423265, -8.881784197001252e-16, -1.0) ``` @@ -474,11 +473,11 @@ julia> f(x) = tan(x); g(x) = x/(B*(Λ*x^2 - 1)); julia> h(x) = f(x) - g(x) h (generic function with 1 method) -julia> x = find_zero(h, (k*pi, (k + 1/2)*pi)); x, h(x) +julia> c = find_zero(h, (k*pi, (k + 1/2)*pi)); c, h(c) (9.530477156207574, 8.326672684688674e-16) ``` -As of version 1.9 of `Julia`, an extension is provided so that when `SymPy` is loaded, an equation can be used to specify the left and right hand sides, as with: +As of version 1.9 of `Julia`, an extension is provided so that when `SymPy` is loaded, an equation of a single symbolic variable can be used to specify the left and right hand sides, as with: ``` using SymPy @@ -1039,20 +1038,14 @@ produce zeros as exact as possible, and the fact that the error in function evaluation, $f'(x)|x|\epsilon$, is not typically on the scale of `1e-8`, leads to a desire for more precision, if available. -In `Roots`, the faster algorithms use a check on both the size of -`f(xn)` and the size of the difference between the last two `xn` values. The check on `f(xn)` -is done with a tight tolerance, as is the check on $x_n \approx -x_{n-1}$. If the function values get close to zero, an -approximate zero is declared. Further, if the $x$ values get close to each other -*and* the function value is close to zero with a *relaxed* tolerance, -then an approximate zero is declared. In practice this seems to work -reasonably well. The relaxed tolerance uses the cube root of the -absolute and relative tolerances. - - - - - +In `Roots`, the faster algorithms (and a few bracketing intervals) use a check on both the size the +size of the difference between the last two `xn` values and the size +of the residual `f(xn)`. The check on `f(xn)` is done with a tight +tolerance, as is the check on $x_n \approx x_{n-1}$. If the function +values get close to zero, an approximate zero is declared. Further, if +the $x$ values get close to each other *and* the function value is +close to zero with a *relaxed* tolerance, then an approximate zero is +declared. In practice this seems to work reasonably well. ## Searching for all zeros in an interval @@ -1208,11 +1201,15 @@ To add a solver the minimum needed is a type to declare the solver and an `updat The [Wikipedia](https://en.wikipedia.org/wiki/Brent%27s_method) page for Brent's method suggest a modern improvement, Chandrapatla's method, described [here](https://www.google.com/books/edition/Computational_Physics/cC-8BAAAQBAJ?hl=en&gbpv=1&pg=PA95&printsec=frontcover). That description is mostly followed below and in the package implementation `Roots.Chandrapatla`. -To implement Chandrapatla's algorithm we first define a type to indicate the method and a state object which records the values ``x_n``, ``x_{n-1}``, and ``x_{n-2}``, needed for the inverse quadratic step. +To implement Chandrapatla's algorithm we first define a type to indicate the method: ```julia julia> struct Chandrapatla <: Roots.AbstractBracketingMethod end +``` + +For this method, the default state object isn't sufficient, as the algorithm tracks three values, not the default two. Here we define a state object for this type: +```julia julia> struct ChandrapatlaState{T,S} <: Roots.AbstractUnivariateZeroState{T,S} xn1::T xn0::T @@ -1227,7 +1224,7 @@ end An `init_state` method can be used by some methods to add more detail to the basic state object. Here it starts the old value, `c`, off as `a` as a means to ensure an initial bisection step. ```julia -julia> function init_state(::Chandrapatla, F, x₀, x₁, fx₀, fx₁) +julia> function Roots.init_state(::Chandrapatla, F, x₀, x₁, fx₀, fx₁) a, b, fa, fb = x₁, x₀, fx₁, fx₀ c, fc = a, fa ChandrapatlaState(b, a, c, fb, fa, fc) @@ -1250,7 +1247,7 @@ julia> function Roots.update_state(::Chandrapatla, F, o, options, l=Roots.NullTr ϕ² = ϕ^2 Δ = (ϕ² < ξ) && (1 - 2ϕ + ϕ² < 1 - ξ) # Chandrapatla's inequality to determine next step - xₜ = Δ ? Roots.inverse_quadratic_step(a, b, c, fa, fb, fc) : a + (b-a)/2 + xₜ = Δ ? Roots.inverse_quadratic_step(a, b, c, fa, fb, fc) : a/2 + b/2 fₜ = F(xₜ) incfn(l) diff --git a/src/Bracketing/alefeld_potra_shi.jl b/src/Bracketing/alefeld_potra_shi.jl index 61497815..7b634a7a 100644 --- a/src/Bracketing/alefeld_potra_shi.jl +++ b/src/Bracketing/alefeld_potra_shi.jl @@ -21,7 +21,7 @@ This implementation deviates slightly from the printed algorithm, as it may use These algorithms do not check the size of `f` for convergence, so the `atol` or `rtol` are not utilized. """ =# -abstract type AbstractAlefeldPotraShi <: AbstractBracketingMethod end +abstract type AbstractAlefeldPotraShi <: AbstractNonStrictBracketingMethod end initial_fncalls(::AbstractAlefeldPotraShi) = 3 # worst case assuming fx₀, fx₁,fc must be computed @@ -134,11 +134,9 @@ function update_state( ā, b̄, d, fā, fb̄, fd = bracket(a, b, x, fa, fb, fx) - if ( - (b̄ - ā) <= tolₑ(ā, b̄, fā, fb̄, atol, rtol) || + if ((b̄ - ā) <= tolₑ(ā, b̄, fā, fb̄, atol, rtol) || iszero(fx) || # exact zero - !isbracket(fā, fb̄) - ) # catch non bracket?, issue #453 + !isbracket(fā, fb̄)) # catch non bracket?, issue #453 @reset o.xn0 = ā @reset o.xn1 = b̄ @reset o.fxn0 = fā diff --git a/src/Bracketing/bracketing.jl b/src/Bracketing/bracketing.jl index 6a2bfe2e..dbf83d64 100644 --- a/src/Bracketing/bracketing.jl +++ b/src/Bracketing/bracketing.jl @@ -54,6 +54,20 @@ function default_tolerances( (xatol, xrtol, atol, rtol, maxevals, strict) end +# Non strict has check on non convergence with relaxed tolerances +function default_tolerances( + ::AbstractNonStrictBracketingMethod, + ::AbstractUnivariateZeroState{T,S}, +) where {T,S} + xatol = eps(real(T))^3 * oneunit(real(T)) + xrtol = eps(real(T)) # unitless + atol = zero(oneunit(real(S))) + rtol = zero(one(real(S))) + maxevals = 60 + strict = false + (xatol, xrtol, atol, rtol, maxevals, strict) +end + ## -------------------------------------------------- const bracketing_error = """The interval [a,b] is not a bracketing interval. diff --git a/src/Bracketing/brent.jl b/src/Bracketing/brent.jl index 78892076..9d3b2dee 100644 --- a/src/Bracketing/brent.jl +++ b/src/Bracketing/brent.jl @@ -7,7 +7,7 @@ This method uses a choice of inverse quadratic interpolation or a secant step, falling back on bisection if necessary. """ -struct Brent <: AbstractBracketingMethod end +struct Brent <: AbstractNonStrictBracketingMethod end struct BrentState{T,S} <: AbstractUnivariateZeroState{T,S} xn1::T diff --git a/src/Bracketing/chandrapatlu.jl b/src/Bracketing/chandrapatlu.jl index c506396d..c6a77f06 100644 --- a/src/Bracketing/chandrapatlu.jl +++ b/src/Bracketing/chandrapatlu.jl @@ -10,7 +10,7 @@ Chandrapatla's algorithm chooses between an inverse quadratic step or a bisectio """ -struct Chandrapatla <: AbstractBracketingMethod end +struct Chandrapatla <: AbstractNonStrictBracketingMethod end struct ChandrapatlaState{T,S} <: AbstractUnivariateZeroState{T,S} xn1::T diff --git a/src/Bracketing/false_position.jl b/src/Bracketing/false_position.jl index bfacc31a..daa69522 100644 --- a/src/Bracketing/false_position.jl +++ b/src/Bracketing/false_position.jl @@ -1,4 +1,4 @@ -struct FalsePosition{R} <: AbstractSecantMethod end +struct FalsePosition{R} <: AbstractRegulaFalsiMethod end """ @@ -30,27 +30,13 @@ Examples ``` find_zero(x -> x^5 - x - 1, (-2, 2), FalsePosition()) ``` + +Note: The alternative `RegulaFalsi` method provides most of this and seems to suffer less from numerical issues. + """ FalsePosition FalsePosition(x=:anderson_bjork) = FalsePosition{x}() -# 12 is tough; needs more evaluations -function default_tolerances( - ::FalsePosition{12}, - ::AbstractUnivariateZeroState{T,S}, -) where {T,S} - xatol = eps(real(T)) * oneunit(real(T)) - xrtol = eps(real(T)) # unitless - atol = 4 * eps(real(float(S))) * oneunit(real(S)) - rtol = 4 * eps(real(float(S))) * one(real(S)) - maxiters = 250 - strict = false - (xatol, xrtol, atol, rtol, maxiters, strict) -end - -init_state(M::FalsePosition, F, x₀, x₁, fx₀, fx₁) = - init_state(Bisection(), F, x₀, x₁, fx₀, fx₁) - function update_state( method::FalsePosition, fs, @@ -62,7 +48,7 @@ function update_state( fa, fb = o.fxn0, o.fxn1 lambda = fb / (fb - fa) - ϵ = √eps(T) / 100 # some engineering to avoid short moves; still fails on some + ϵ = √eps(T) * abs(b-a) # some engineering to avoid short moves; still fails on some ϵ ≤ lambda ≤ 1 - ϵ || (lambda = 1 / 2) x::T = b - lambda * (b - a) fx::S = fs(x) @@ -82,33 +68,29 @@ function update_state( return (o, false) end -# the 12 reduction factors offered by Galdino -# In RootsTesting.jl, we can see :12 has many more failures. -galdino = Dict{Union{Int,Symbol},Function}( - :1 => (fa, fb, fx) -> fa * fb / (fb + fx), - :2 => (fa, fb, fx) -> (fa - fb) / 2, - :3 => (fa, fb, fx) -> (fa - fx) / (2 + fx / fb), - :4 => (fa, fb, fx) -> (fa - fx) / (1 + fx / fb)^2, - :5 => (fa, fb, fx) -> (fa - fx) / (3*one(fa)/2 + fx / fb)^2, - :6 => (fa, fb, fx) -> (fa - fx) / (2 + fx / fb)^2, - :7 => (fa, fb, fx) -> (fa + fx) / (2 + fx / fb)^2, - :8 => (fa, fb, fx) -> fa / 2, - :9 => (fa, fb, fx) -> fa / (1 + fx / fb)^2, - :10 => (fa, fb, fx) -> (fa - fx) / 4, - :11 => (fa, fb, fx) -> fx * fa / (fb + fx), - :12 => (fa, fb, fx) -> (fa * (1 - fx / fb > 0 ? 1 - fx / fb : one(fa)/2)), -) - -# give common names -for (nm, i) in [(:pegasus, 1), (:illinois, 8), (:anderson_bjork, 12)] - galdino[nm] = galdino[i] -end - -# from Chris Elrod; https://raw.githubusercontent.com/chriselrod/AsymptoticPosteriors.jl/master/src/false_position.jl -@generated function galdino_reduction(methods::FalsePosition{R}, fa, fb, fx) where {R} - f = galdino[R] - quote - $Expr(:meta, :inline) - $f(fa, fb, fx) - end -end +galdino_reduction(::FalsePosition{1}, fa, fb, fx) = fa * fb / (fb + fx) +galdino_reduction(::FalsePosition{2}, fa, fb, fx) = (fa - fb) / 2 +galdino_reduction(::FalsePosition{3}, fa, fb, fx) = (fa - fx) / (2 + fx / fb) +galdino_reduction(::FalsePosition{4}, fa, fb, fx) = (fa - fx) / (1 + fx / fb)^2 +galdino_reduction(::FalsePosition{5}, fa, fb, fx) = (fa - fx) / (3*one(fa)/2 + fx / fb)^2 +galdino_reduction(::FalsePosition{6}, fa, fb, fx) = (fa - fx) / (2 + fx / fb)^2 +galdino_reduction(::FalsePosition{7}, fa, fb, fx) = (fa + fx) / (2 + fx / fb)^2 +galdino_reduction(::FalsePosition{8}, fa, fb, fx) = fa / 2 +galdino_reduction(::FalsePosition{9}, fa, fb, fx) = fa / (1 + fx / fb)^2 +galdino_reduction(::FalsePosition{10}, fa, fb, fx) = (fa - fx) / 4 +galdino_reduction(::FalsePosition{11}, fa, fb, fx) = fx * fa / (fb + fx) +galdino_reduction(::FalsePosition{12}, fa, fb, fx) = + (fa * (1 - fx / fb > 0 ? 1 - fx / fb : one(fa)/2)) + +galdino_reduction(::FalsePosition{:pegasus}, fa, fb, fx) = + galdino_reduction(FalsePosition(1), fa, fb, fx) +galdino_reduction(::FalsePosition{:Pegasus}, fa, fb, fx) = + galdino_reduction(FalsePosition(1), fa, fb, fx) +galdino_reduction(::FalsePosition{:illinois}, fa, fb, fx) = + galdino_reduction(FalsePosition(8), fa, fb, fx) +galdino_reduction(::FalsePosition{:Illinois}, fa, fb, fx) = + galdino_reduction(FalsePosition(8), fa, fb, fx) +galdino_reduction(::FalsePosition{:anderson_bjork}, fa, fb, fx) = + galdino_reduction(FalsePosition(12), fa, fb, fx) +galdino_reduction(::FalsePosition{:AndersonBjork}, fa, fb, fx) = + galdino_reduction(FalsePosition(12), fa, fb, fx) diff --git a/src/Bracketing/itp.jl b/src/Bracketing/itp.jl index 550e0b00..f2237273 100644 --- a/src/Bracketing/itp.jl +++ b/src/Bracketing/itp.jl @@ -21,7 +21,7 @@ Suggested on by `@TheLateKronos`, who supplied the original version of the code. """ -struct ITP{T,S} <: AbstractBracketingMethod +struct ITP{T,S} <: AbstractNonStrictBracketingMethod κ₁::T κ₂::S n₀::Int diff --git a/src/Bracketing/modAB.jl b/src/Bracketing/modAB.jl index 4d99694b..576ba754 100644 --- a/src/Bracketing/modAB.jl +++ b/src/Bracketing/modAB.jl @@ -48,7 +48,7 @@ julia> F = Cnt(f); x = find_zero(F, (a,b), Roots.ModAB()); (x, F.cnt) * That the value found in the example by `ModAB` and `Bisection` agree only through the first 5 decimal points is due to the `ModAB` algorithm stopping on small `f(xₙ)` values, as `Bisection` iterates up to the last floating point bit unless it finds an exact numeric zero. -## Reference +## References N Ganchovski and A Traykov 2023 IOP Conf. Ser.: Mater. Sci. Eng. 1276 012010 @@ -56,6 +56,10 @@ DOI 10.1088/1757-899X/1276/1/012010 [https://iopscience.iop.org/article/10.1088/1757-899X/1276/1/012010](https://iopscience.iop.org/article/10.1088/1757-899X/1276/1/012010) +Improvements to the Modified Anderson–Björck (modAB) Root-Finding Algorithm +by Nedelcho Ganchovski, Oscar Smith,Christopher Rackauckas, Lachezar Tomov, and Alexander Traykov +[https://doi.org/10.3390/a19050332](https://doi.org/10.3390/a19050332) + A new high order method of regula falsi type for computing a root of an equation; Ned Anderson & Åke Björck. [https://link.springer.com/article/10.1007/BF01951936](https://link.springer.com/article/10.1007/BF01951936). diff --git a/src/Bracketing/regula_falsi.jl b/src/Bracketing/regula_falsi.jl new file mode 100644 index 00000000..3845f972 --- /dev/null +++ b/src/Bracketing/regula_falsi.jl @@ -0,0 +1,134 @@ +struct RegulaFalsi{R} <: AbstractRegulaFalsiMethod end + +""" + RegulaFalsi{M} + ScalingFactorRegulaFalsi(M::Symbol=:AndersonBjork) + SCRF{M} + +Implements several different *scaling factor* *regula falsi* methods following Diaz and Perez. The *regula falsi* method uses the secant method to find the next value with a bracketing interval providing the history. This is a linearly convergent method. By scaling, or adjusting, the ``f(x_i)`` values, superlinear convergence can be achieved. + +The scaling factor refers to the `f(xᵢ)` values being scaled by `γ` which is dynamically computed based on `ξ = f(xᵢ₊₁) / f(xᵢ)` and `ζ = - f(xᵢ₊₁) / f(xᵢ₋₁)`. + +The value for `M` is one of `Roots._regula_falsi_names` or `(:classic, :Illinois, :GIll01, :Pegasus, :AndersonBjork, :AB_GIll01, :Ford3, :Ford4)`. The default is :AndersonBjork. + +# Examples + +``` +find_zero(x -> x^5 - x - 1, (-2, 2), RegulaFalsi()) # default :AndersonBjork +find_zero(x -> x^5 - x - 1, (-2, 2), RegulaFalsi(:Illinois)) +``` + +New scaling factors can be introduced by defining a method for `Fᵧ(::ScalingFactorRegulaFalsi{M}, ξ::T, ζ::T) where T` and then calling with `ScalingFactorRegulaFalsi{M}` for `M` a symbol. + +Convergence rates from the reference paper are: + +1. `:classic`---Classic *regula falsi*: p = 1 (linear). +2. `:Illinois`---Illinois method: p = 1.442. +3. `:Pegasus`---Pegasus method: p = 1.642. +4. `:AndersonBjork`---A&B method: 1.681 < p < 1.710. +5. `:Ford4`---Ford fourth method: p = 1.681. + + +# References: + +*A common framework for modified Regula Falsi methods and new methods of this kind* by Julio M. Fernández-Díaz, César O. Menéndez-Pérez [url](https://digibuo.uniovi.es/dspace/bitstream/handle/10651/66581/1-s2.0-S0378475422004335-main.pdf) + +!!! note + Compared to similar methods in `FalsePosition` the parameterization used here has fewer issues with floating point differences and should be preferred. + +""" +RegulaFalsi +RegulaFalsi(x=:AndersonBjork) = RegulaFalsi{x}() +const ScalingFactorRegulaFalsi{R} = RegulaFalsi{R} +const SCRF{R} = RegulaFalsi{R} + +# for testing, this might be helpful +_regula_falsi_names = + (:classic, :Illinois, :GIll01, :Pegasus, :AndersonBjork, :AB_GIll01, :Ford3, :Ford4) + +# factors +Fᵧ(::RegulaFalsi{:classic}, ξ::T, ζ::T) where {T} = one(T) +Fᵧ(::RegulaFalsi{:Illinois}, ξ::T, ζ::T) where {T} = one(T)/2 +Fᵧ(::RegulaFalsi{:GIll01}, ξ::T, ζ::T) where {T} = one(T)/10 # Generalized Illinois method with γ = 0.1 +Fᵧ(::RegulaFalsi{:Pegasus}, ξ::T, ζ::T) where {T} = one(T) / (one(T) + ξ) +Fᵧ(::RegulaFalsi{:AndersonBjork}, ξ::T, ζ::T) where {T} = ξ < 1 ? (one(T) - ξ) : one(T)/2 +Fᵧ(::RegulaFalsi{:AB_GIll01}, ξ::T, ζ::T) where {T} = max(one(T) - ξ, one(T)/10) # AB + Generalized Illinois method with γ = 0.1, though paper conditions on no. of steps +Fᵧ(::RegulaFalsi{:Ford3}, ξ::T, ζ::T) where {T} = (one(T) - ξ + ζ) / (one(T) - ζ) +Fᵧ(::RegulaFalsi{:Ford4}, ξ::T, ζ::T) where {T} = one(T) - ξ + ζ + +# take one step so the :right/:left is set up +function init_state( + M::AbstractRegulaFalsiMethod, + F, + x₀::T, + x₁::T, + fx₀::S, + fx₁::S, +) where {T,S} + assert_bracket(fx₀, fx₁) + + c::T = (x₀ * fx₁ - x₁ * fx₀) / (fx₁ - fx₀) + fc::S = F(c) + + if sign(fx₀) * sign(fc) < 0 + a, fa = x₀, fx₀ + else + a, fa = x₁, fx₁ + end + UnivariateZeroState(c, a, fc, fa) +end + +initial_fncalls(M::AbstractRegulaFalsiMethod) = 3 + +function default_tolerances( + ::AbstractRegulaFalsiMethod, + ::AbstractUnivariateZeroState{T,S}, +) where {T,S} + xatol = eps(T)^3 * oneunit(T) + xrtol = eps(T) * one(T) # unitless + atol = 0 * oneunit(S) + rtol = 0 * one(S) + maxiters = 250 + strict = false + (xatol, xrtol, atol, rtol, maxiters, strict) +end + +function update_state( + M::RegulaFalsi, + fs, + o::AbstractUnivariateZeroState{T,S}, + options, + l=NullTracks(), +) where {T,S} + xₙ₋₁::T, xₙ::T = o.xn0, o.xn1 + f̂xₙ₋₁::S, f̂xₙ::S = o.fxn0, o.fxn1 + + c::T = (xₙ₋₁ * f̂xₙ - xₙ * f̂xₙ₋₁) / (f̂xₙ - f̂xₙ₋₁) + + ϵ = abs(xₙ - xₙ₋₁) * sqrt(eps(T)) # some engineering to avoid short moves + if abs(c - xₙ) ≤ ϵ || abs(c - xₙ₋₁) ≤ ϵ + c = xₙ₋₁/2 + xₙ/2 + end + + fc::S = fs(c) + incfn(l) + + iszero(fc) && return (_set(o, (c, fc)), true) + + xₙ₊₁, f̂xₙ₊₁ = c, fc + + if sign(f̂xₙ) * sign(fc) > 0 # stuck on same side as last time, scale other side + ξ::S = fc / f̂xₙ # \xi + ζ::S = - fc / f̂xₙ₋₁ # \zeta + γ::S = Fᵧ(M, ξ, ζ) + γ = max(zero(ξ), min(one(ξ), γ)) + + xₙ = xₙ₋₁ + f̂xₙ = γ * f̂xₙ₋₁ + iszero(f̂xₙ) && (f̂xₙ = f̂xₙ₋₁/2) # Some engineering to avoid early termination by f̂(xₙ) == 0.0 + end + + o = _set(o, (xₙ₊₁, f̂xₙ₊₁), (xₙ, f̂xₙ)) + + return (o, false) +end diff --git a/src/Bracketing/ridders.jl b/src/Bracketing/ridders.jl index 5d230968..ad7d8a33 100644 --- a/src/Bracketing/ridders.jl +++ b/src/Bracketing/ridders.jl @@ -23,7 +23,7 @@ true `f=F', g=F''/2, h=F'''/6`, suggesting converence at rate `≈ 1.839...`. It uses two function evaluations per step, so its order of convergence is `≈ 1.225...`. """ -struct Ridders <: AbstractBracketingMethod end +struct Ridders <: AbstractNonStrictBracketingMethod end function update_state( M::Ridders, diff --git a/src/DerivativeFree/sidi.jl b/src/DerivativeFree/sidi.jl index 0edb0c71..7b7f1565 100644 --- a/src/DerivativeFree/sidi.jl +++ b/src/DerivativeFree/sidi.jl @@ -24,6 +24,7 @@ find_zero(sin, 3, Roots.Sidi(2)) """ struct Sidi{k} <: AbstractSecantMethod end Sidi(k::Int) = Sidi{k}() +initial_fncalls(::Sidi{K}) where {K} = K + 1 struct SidiState{T,S} <: AbstractUnivariateZeroState{T,S} xn1::T @@ -37,7 +38,7 @@ end function init_state(M::Sidi{k}, F::Callable_Function, x) where {k} x₀, x₁ = x₀x₁(x) fx₀, xs, fs = _init_sidi(F, (x₀, x₁), k) - state = SidiState(xs[k], xs[k + 1], fx₀, fs[1], xs, fs) + state = SidiState(xs[k + 1], xs[k], fs[1], fx₀, xs, fs) end function update_state( @@ -103,6 +104,7 @@ function _init_sidi(f, x, k) xⱼ = xs[j] = xⱼ₋₁ - fs[1] / pk′ end Δ = f(xⱼ) + j == k && (fx₀ = Δ) # initial value for x₁ for i in 2:j Δ₀ = fs[i - 1] fs[i - 1] = Δ diff --git a/src/Roots.jl b/src/Roots.jl index f623b14c..77c0ec82 100644 --- a/src/Roots.jl +++ b/src/Roots.jl @@ -23,8 +23,6 @@ import CommonSolve import CommonSolve: solve, solve!, init using Accessors -export fzero, fzeros, secant_method - export find_zero, find_zeros, ZeroProblem, @@ -60,6 +58,7 @@ include("Bracketing/brent.jl") include("Bracketing/ridders.jl") include("Bracketing/itp.jl") include("Bracketing/chandrapatlu.jl") +include("Bracketing/regula_falsi.jl") include("Bracketing/false_position.jl") include("Bracketing/modAB.jl") @@ -80,11 +79,16 @@ include("Derivative/thukralb.jl") include("Derivative/lith.jl") include("find_zeros.jl") -include("simple.jl") -include("alternative_interfaces.jl") if !isdefined(Base, :get_extension) include("../ext/RootsChainRulesCoreExt.jl") end +## ---- deprecated +include("deprecated/simple.jl") +export secant_method +include("deprecated/alternative_interfaces.jl") +export fzero, fzeros + + end diff --git a/src/abstract_types.jl b/src/abstract_types.jl index c8cf2279..e16c4588 100644 --- a/src/abstract_types.jl +++ b/src/abstract_types.jl @@ -4,6 +4,8 @@ Base.broadcastable(method::AbstractUnivariateZeroMethod) = Ref(method) abstract type AbstractBracketingMethod <: AbstractUnivariateZeroMethod end abstract type AbstractBisectionMethod <: AbstractBracketingMethod end +abstract type AbstractRegulaFalsiMethod <: AbstractBracketingMethod end +abstract type AbstractNonStrictBracketingMethod <: AbstractBracketingMethod end abstract type AbstractNonBracketingMethod <: AbstractUnivariateZeroMethod end abstract type AbstractSecantMethod <: AbstractNonBracketingMethod end diff --git a/src/convergence.jl b/src/convergence.jl index e0b2bb05..c70f74e0 100644 --- a/src/convergence.jl +++ b/src/convergence.jl @@ -87,10 +87,13 @@ end ## test f ≈ 0 not f == 0 function is_approx_zero_f( - ::AbstractUnivariateZeroMethod, + ::M, state::AbstractUnivariateZeroState, options::O, -) where {O<:AbstractUnivariateZeroOptions} +) where { + M<:Union{AbstractUnivariateZeroMethod,AbstractRegulaFalsiMethod}, + O<:AbstractUnivariateZeroOptions, +} ab, afb = abs(state.xn1), abs(state.fxn1) ϵₐ, ϵᵣ = options.abstol, options.reltol Δ = max(_unitless(ϵₐ), _unitless(ab) * ϵᵣ) @@ -214,6 +217,8 @@ If algorithm hasn't converged this returns `(:not_converged, false)`. If algorithm has stopped or converged, return flag and `true`. Flags are: +* :`exact_zero` if `|f(xn1)| == 0` + * `:x_converged` if `xn1 ≈ xn`, typically with non-zero tolerances specified. * `:f_converged` if `|f(xn1)| < max(atol, |xn1|*rtol)` @@ -341,29 +346,53 @@ end function decide_convergence( M::AbstractBracketingMethod, F, - state::AbstractUnivariateZeroState, + state::AbstractUnivariateZeroState{T,S}, options, val, -) - a, b = state.xn0, state.xn1 - fa, fb = state.fxn0, state.fxn1 +) where {T,S} + b, a = xs = state.xn1, state.xn0 + fb, fa = fxs = state.fxn1, state.fxn0 + m, i = findmin(abs, fxs) + α = xs[i] + + if val == :exact_zero + return α + elseif val == :x_converged + # when exact closeness is *possible* + # get as close as possible with one extra function call + u, v, fu, fv = a < b ? (a, b, fa, fb) : (b, a, fb, fa) + u₊ = nextfloat(float(u)) + u₊₊ = nextfloat(u₊) + if v == u₊₊ + fu₊ = first(F(u₊)) + m, i = findmin(abs, (fu, fu₊, fv)) + return (u, u₊, v)[i] + else + return α + end + elseif val == :not_converged + if !options.strict + # check relaxed convergence on Δx, f(x) + # Δ ≤ 2^8 * max(δ + |x|⋅ϵ) + # |f(x)| ≤ 16*min(16δ + |x|⋅ϵ) -- 0 by default + mx = maximum(abs, xs) + + Δ = abs(xs[1] - xs[2]) + ϵ = 256 * max(options.xabstol, mx * sqrt(options.xreltol)) + _unitless(Δ) ≤ _unitless(ϵ) && return α + + atol = max(eps(oneunit(real(S))), options.abstol) + rtol = max(eps(one(real(S))), options.xreltol) + δ = 16 * min(16*atol, maximum(abs, xs) * rtol) + _unitless(m) ≤ _unitless(δ) && return α + end + return nan(T) * state.xn1 + end - iszero(fa) && return a - iszero(fb) && return b isnan(fa) && return a isnan(fb) && return b + iszero(fa) && return a + iszero(fb) && return b - # get as close as possible with one extra function call - # when exact closeness is possible - a₊₊ = nextfloat(nextfloat(float(a))) - if b == a₊₊ - c = nextfloat(float(a)) - fc = first(F(c)) - m = minimum(abs, (fa, fb, fc)) - abs(fc) == m && return c - abs(fa) == m && return a - return b - end - - return (abs(fa) < abs(fb)) ? a : b + return α end diff --git a/src/alternative_interfaces.jl b/src/deprecated/alternative_interfaces.jl similarity index 53% rename from src/alternative_interfaces.jl rename to src/deprecated/alternative_interfaces.jl index 6a0c6962..06dafce3 100644 --- a/src/alternative_interfaces.jl +++ b/src/deprecated/alternative_interfaces.jl @@ -1,90 +1,6 @@ ## Some **legacy** alternative interfaces. - +## These should be deprecated! ## several named interfaces to methods -## newton, halley, quadratic_inverse, superhalley, chebyshevlike -""" - Roots.newton(f, fp, x0; kwargs...) - -Implementation of Newton's method: `xᵢ₊₁ = xᵢ - f(xᵢ)/f'(xᵢ)`. - -Arguments: - -* `f::Function` -- function to find zero of - -* `fp::Function` -- the derivative of `f`. - -* `x0::Number` -- initial guess. For Newton's method this may be complex. - -With the `ForwardDiff` package derivatives may be computed automatically. For example, defining -`D(f) = x -> ForwardDiff.derivative(f, float(x))` allows `D(f)` to be used for the first derivative. - -Keyword arguments are passed to `find_zero` using the `Roots.Newton()` method. - -See also `Roots.newton((f,fp), x0)` and `Roots.newton(fΔf, x0)` for simpler implementations. - -""" -newton(f, fp, x0; kwargs...) = find_zero((f, fp), x0, Newton(); kwargs...) - -## -------------------------------------------------- -#= -""" - Roots.halley(f, fp, fpp, x0; kwargs...) - -Implementation of Halley's method (cf `?Roots.Halley()`). - -Arguments: - -* `f::Function` -- function to find zero of - -* `fp::Function` -- derivative of `f`. - -* `fpp:Function` -- second derivative of `f`. - -* `x0::Number` -- initial guess - -With the `ForwardDiff` package derivatives may be computed automatically. For example, defining -`D(f) = x -> ForwardDiff.derivative(f, float(x))` allows `D(f)` and `D(D(f))` to be used for the first and second -derivatives, respectively. - -Keyword arguments are passed to `find_zero` using the `Roots.Halley()` method. - -""" -=# -halley(f, fp, fpp, x0; kwargs...) = find_zero((f, fp, fpp), x0, Halley(); kwargs...) - -#= -""" - Roots.quadratic_inverse(f, fp, fpp, x0; kwargs...) - -Implementation of the quadratic inverse method (cf `?Roots.QuadraticInverse()`). - -Arguments: - -* `f::Function` -- function to find zero of - -* `fp::Function` -- derivative of `f`. - -* `fpp:Function` -- second derivative of `f`. - -* `x0::Number` -- initial guess - -With the `ForwardDiff` package derivatives may be computed automatically. For example, defining -`D(f) = x -> ForwardDiff.derivative(f, float(x))` allows `D(f)` and `D(D(f))` to be used for the first and second -derivatives, respectively. - -Keyword arguments are passed to `find_zero` using the `Roots.QuadraticInverse()` method. - -""" -=# -quadratic_inverse(f, fp, fpp, x0; kwargs...) = - find_zero((f, fp, fpp), x0, QuadraticInverse(); kwargs...) - -superhalley(f, fp, fpp, x0; kwargs...) = - find_zero((f, fp, fpp), x0, SuperHalley(); kwargs...) - -chebyshev_like(f, fp, fpp, x0; kwargs...) = - find_zero((f, fp, fpp), x0, ChebyshevLike(); kwargs...) - ## -------------------------------------------------- ## MATLAB interface to find_zero @@ -152,10 +68,25 @@ fzero(sin, cos, 3) # use Newton's method This has the advantage of making the first use of the function `f` faster, but subsequent uses slower. """ -function fzero(f, x0::Number; kwargs...) - x = float(x0) - isinf(x) && throw(ConvergenceFailed("An initial value must be finite")) - derivative_free(f, x; kwargs...) +function fzero(f, x0; order=nothing, verbose=false, tracks=NullTracks(), kwargs...) + #Base.depwarn("`fzero(f, x0)` is deprecated; use `find_zero(f, x0)` instead.", :fzero) + + M = haskey(_method_lookup, order) ? _method_lookup[order] : + isa(x0, Tuple) ? Bisection() : + isa(x0, Number) ? Order0() : + Order0() + if verbose && isa(tracks, NullTracks) + tracks = Tracks() + end + @noinline α = find_zero(FnWrapper(f), x0, M; tracks, kwargs...) + verbose && display(tracks) + α +end + + +function fzero(f, a::Number, b::Number, args...; kwargs...) + #Base.depwarn("`fzero(f, a, b)` is deprecated; use `find_zero(f, (a,b))` instead.", :fzero) + fzero(f, (a, b), args...; kwargs...) end function fzero( @@ -166,6 +97,7 @@ function fzero( tracks=NullTracks(), kwargs..., ) + #Base.depwarn("`fzero(f, x0, M)` is deprecated; use `find_zero(f, x0, M)` instead.", :fzero) tracks = (verbose && isa(tracks, NullTracks)) ? Tracks() : tracks α = find_zero(FnWrapper(f), x0, M; tracks, kwargs...) verbose && display(tracks) @@ -181,39 +113,13 @@ function fzero( tracks=NullTracks(), kwargs..., ) + #Base.depwarn("`fzero(f, x0, M, N)` is deprecated; use `find_zero(f, x0, M, N)` instead.", :fzero) tracks = (verbose && isa(tracks, NullTracks)) ? Tracks() : tracks a = find_zero(FnWrapper(f), x0, M, N; tracks, kwargs...) verbose && display(tracks) a end -function fzero( - f, - bracket::Tuple{T,S}; - verbose=false, - tracks=NullTracks(), - kwargs..., -) where {T<:Number,S<:Number} - d = Dict(kwargs...) - tracks = (verbose && isa(tracks, NullTracks)) ? Tracks() : tracks - if haskey(d, :order) - val = find_zero(FnWrapper(f), bracket, _method_lookup[d[:order]]; tracks, kwargs...) - else - val = find_zero(FnWrapper(f), bracket, Bisection(); tracks, kwargs...) - end - verbose && display(tracks) - val -end - -fzero(f, a::Number, b::Number, args...; kwargs...) = fzero(f, (a, b), args...; kwargs...) - -function fzero(f, x; verbose=false, tracks=NullTracks(), kwargs...) - tracks = (verbose && isa(tracks, NullTracks)) ? Tracks() : tracks - α = find_zero(FnWrapper(f), x; kwargs...) - verbose && display(tracks) - α -end - function fzero( f::Function, fp::Function, @@ -222,8 +128,10 @@ function fzero( tracks=NullTracks(), kwargs..., ) + #Base.depwarn("`fzero(f, fp, x0)` is deprecated; use `find_zero((f, fp), x0, Roots.Newton())` instead.", :fzero) + M = Newton() tracks = (verbose && isa(tracks, NullTracks)) ? Tracks() : tracks - α = find_zero((f, fp), x0, Newton(); tracks, kwargs...) + α = find_zero((f, fp), x0, M; tracks, kwargs...) verbose && display(tracks) α end @@ -260,38 +168,7 @@ _method_lookup = Dict( "16" => Order16(), ) -@noinline function derivative_free( - f, - x0; - verbose::Bool=false, - tracks=NullTracks(), - order=0, - kwargs..., -) - if haskey(_method_lookup, order) - M = _method_lookup[order] - else - throw(ArgumentError("Invalid order specified. See ?fzero.")) - end - - # d = (kv[1] == :ftol ? :atol=>kv[2] : - # kv[1] == :ftolrel ? :rtol=>kv[2] : - # kv[1] == :xtol ? :xatol=>kv[2] : - # kv[1] == :xtolrel ? xrtol=>kv[2] : - # kv[1] => kv[1] for kv in kwargs) - - d = Dict(kwargs...) - for (o, n) in ((:ftol, :atol), (:ftolrel, :rtol), (:xtol, :xatol), (:xtolrel, :xrtol)) - if haskey(d, o) - d[n] = d[o] - end - end - tracks = (verbose && isa(tracks, NullTracks)) ? Tracks() : tracks - a = find_zero(FnWrapper(f), x0, M; tracks, d...) - verbose && display(tracks) - a -end ## fzeros """ @@ -303,6 +180,7 @@ Searches for all zeros of `f` within an interval `(a,b)`. Assumes neither `a` or Compatibility interface for [`find_zeros`](@ref). """ function fzeros(f, a::Number, b::Number; kwargs...) + #Base.depwarn("`fzeros(f, a, b)` is deprecated; use `find_zeros(f, (a,b)) instead.", :fzeros) find_zeros(FnWrapper(f), float(a), float(b); kwargs...) end fzeros(f, ab; kwargs...) = fzeros(f, _extrema(ab)...; kwargs...) diff --git a/src/simple.jl b/src/deprecated/simple.jl similarity index 61% rename from src/simple.jl rename to src/deprecated/simple.jl index 30a8421b..5ffcc2e5 100644 --- a/src/simple.jl +++ b/src/deprecated/simple.jl @@ -10,88 +10,6 @@ # `Roots.dfree(f, xs)` (Order0) more robust secant method # -## Bisection -## -## Essentially from Jason Merrill https://gist.github.com/jwmerrill/9012954 -## cf. http://squishythinking.com/2014/02/22/bisecting-floats/ -## This also borrows a trick from https://discourse.julialang.org/t/simple-and-fast-bisection/14886 -## where we keep x1 so that y1 is negative, and x2 so that y2 is positive -## this allows the use of signbit over y1*y2 < 0 which avoid < and a multiplication -## this has a small, but noticeable impact on performance. -""" - bisection(f, a, b; [xatol, xrtol]) - -Performs bisection method to find a zero of a continuous -function. - -It is assumed that `(a,b)` is a bracket, that is, the function has -different signs at `a` and `b`. The interval `(a,b)` is converted to floating point -and shrunk when `a` or `b` is infinite. The function `f` may be infinite for -the typical case. If `f` is not continuous, the algorithm may find -jumping points over the x axis, not just zeros. - - -If non-trivial tolerances are specified, the process will terminate -when the bracket `(a,b)` satisfies `isapprox(a, b, atol=xatol, -rtol=xrtol)`. For zero tolerances, the default, for `Float64`, `Float32`, -or `Float16` values, the process will terminate at a value `x` with -`f(x)=0` or `f(x)*f(prevfloat(x)) < 0 ` or `f(x) * f(nextfloat(x)) < -0`. For other number types, the `Roots.A42` method is used. - -""" -function bisection(f, a::Number, b::Number; xatol=nothing, xrtol=nothing) - x1, x2 = adjust_bracket(float.((a, b))) - T = eltype(x1) - - atol = xatol === nothing ? zero(T) : abs(xatol) - rtol = xrtol === nothing ? zero(one(T)) : abs(xrtol) - CT = iszero(atol) && iszero(rtol) ? Val(:exact) : Val(:inexact) - - x1, x2 = float(x1), float(x2) - y1, y2 = f(x1), f(x2) - - _unitless(y1 * y2) >= 0 && error("the interval provided does not bracket a root") - - if isneg(y2) - x1, x2, y1, y2 = x2, x1, y2, y1 - end - - xm = Roots._middle(x1, x2) # for possibly mixed sign x1, x2 - ym = f(xm) - - while true - if has_converged(CT, x1, x2, xm, ym, atol, rtol) - return xm - end - - if isneg(ym) - x1, y1 = xm, ym - else - x2, y2 = xm, ym - end - - xm = Roots.__middle(x1, x2) - ym = f(xm) - end -end - -# -0.0 not returned by __middle, so isneg true on [-Inf, 0.0) -@inline isneg(x::T) where {T<:AbstractFloat} = signbit(x) -@inline isneg(x) = _unitless(x) < 0 - -@inline function has_converged(::Val{:exact}, x1, x2, m, ym, atol, rtol) - iszero(ym) && return true - isnan(ym) && return true - x1 != m && m != x2 && return false - return true -end - -@inline function has_converged(::Val{:inexact}, x1, x2, m, ym, atol, rtol) - iszero(ym) && return true - isnan(ym) && return true - val = abs(x1 - x2) <= atol + max(abs(x1), abs(x2)) * rtol - return val -end #= """ @@ -107,13 +25,14 @@ Not exported """ =# function a42(f, ab; atol=nothing, rtol=nothing, λ=0.7, μ=0.5) + #Base.depwarn("`a42(f, ab)` is deprecated; use `find_zero(f, ab, Roots.AlefeldPotraShi())` instead.", :a42) a, b = adjust_bracket(ab) δ₀ = b - a fa, fb = f(a), f(b) assert_bracket(fa, fb) tols = ( - λ = λ, + λ = λ, atol = isnothing(atol) ? zero(one(a)) : atol, rtol = isnothing(rtol) ? eps(one(a)) : rtol, ) @@ -181,6 +100,97 @@ function a42(f, ab; atol=nothing, rtol=nothing, λ=0.7, μ=0.5) end end +""" + modab(f, left, right, args...; kwargs...) + +Implementation of "Modified Anderson-Bjork’s method for solving non-linear equations in structural mechanics" by N Ganchovski and A Traykov. Code contributed by @Proektsoft-EOOD in issue [#487](https://github.com/JuliaMath/Roots.jl/issues/487) + +This is a modified Anderson-Bjork bracketing algorithm which takes the least number of iterations over a wide-ranging test suite. +""" +function modab(f, left::Real, right::Real, target::Real=0.0; precision::Float64=1e-14, maxIter::Int=200) + #Base.depwarn("`modab(f, left, right)` is deprecated; use `find_zero(f, (left, right), Roots.ModAB())` instead.", :a42) + x1, x2 = min(left, right), max(left, right) + y1 = f(x1) - target + abs(y1) <= precision && return x1 + y2 = f(x2) - target + abs(y2) <= precision && return x2 + eps1 = precision * 1e-3 + eps2 = precision * (x2 - x1) + if abs(target) >= 1 + eps1 *= abs(target) + else + eps1 = 0 + end + side = 0 + bisection = true + C = 16 # safety factor for threshold corresponding to 4 iterations = 2^4 + threshold = x2 - x1 # Threshold to fall back to bisection if AB fails to shrink the interval enough + # calculate k on each bisection step with account for local function properties and symmetry + for i in 1:maxIter + local x3, y3 + if bisection + x3 = (x1 + x2) / 2 + y3 = f(x3) - target # Function value at midpoint + ym = (y1 + y2) / 2 # Ordinate of chord at midpoint + r = 1 - abs(ym / (y2 - y1)) # Symmetry factor + k = r * r # Deviation factor + # Check if the function is close enough to linear + if abs(ym - y3) < k * (abs(y3) + abs(ym)) + bisection = false + threshold = (x2 - x1) * C + end + else + x3 = (x1 * y2 - y1 * x2) / (y2 - y1) + if x3 <= x1 + x3 = x1 + y3 = y1 + elseif x3 >= x2 + x3 = x2 + y3 = y2 + else + y3 = f(x3) - target + end + threshold /= 2 + end + + if abs(y3) <= eps1 || abs(x2 - x1) <= eps2 # Convergence check + return x3 + end + + if sign(y1) == sign(y3) + if side == 1 + m = 1 - y3 / y1 + if m <= 0 + y2 /= 2 + else + y2 *= m + end + elseif !bisection + side = 1 + end + x1, y1 = x3, y3 + else + if side == -1 + m = 1 - y3 / y2 + if m <= 0 + y1 /= 2 + else + y1 *= m + end + elseif !bisection + side = -1 + end + x2, y2 = x3, y3 + end + if x2 - x1 > threshold # in case AB failed to shrink the interval enough + bisection = true + side = 0 + end + end + return NaN +end + +## --- non bracketing """ secant_method(f, xs; [atol=0.0, rtol=8eps(), maxevals=1000]) @@ -193,7 +203,7 @@ given by `b - fb/m` where `m` is the slope of the secant line between The initial values can be specified as a pair of 2, as in `(x₀, x₁)` or `[x₀, x₁]`, or as a single value, `x₁` in which case a value of `x₀` is chosen. -The algorithm returns m when `abs(fm) <= max(atol, abs(m) * rtol)`. +The algorithm returns `m` when `abs(fm) <= max(atol, abs(m) * rtol)`. If this doesn't occur before `maxevals` steps or the algorithm encounters an issue, a value of `NaN` is returned. If too many steps are taken, the current value is checked to see if there is a sign @@ -224,6 +234,7 @@ function secant_method( rtol=8eps(one(float(real(first(xs))))), maxevals=100, ) + #Base.depwarn("`secant_method(f, xs)` is deprecated; use `find_zero(f, xs, Secant())` instead.", :secant_method) if length(xs) == 1 # secant needs x0, x1; only x0 given a = float(xs[1]) @@ -238,6 +249,7 @@ function secant_method( end function secant(f, a::T, b::T, atol=zero(T), rtol=8eps(T), maxevals=100) where {T} + #Base.depwarn("`secant(f, a, b)` is deprecated; use `find_zero(f, (a,b), Secant())` instead.", :secant) nan = (0a) / (0a) cnt = 0 @@ -314,6 +326,7 @@ function muller( xrtol=nothing, maxevals=300, ) where {T} + #Base.depwarn("`muller(f, x)` is deprecated; use `find_zero(f, x, Roots.Muller())` instead.", :muller) @assert old ≠ older ≠ oldest ≠ old # we want q to be non-degenerate xᵢ₋₂, xᵢ₋₁, xᵢ = oldest, older, old fxᵢ₋₂, fxᵢ₋₁ = f(xᵢ₋₂), f(xᵢ₋₁) @@ -409,6 +422,7 @@ If the convergence fails, will return a `ConvergenceFailed` error. =# newton(f::Tuple, x0; kwargs...) = newton(TupleWrapper(f[1], f[2]), x0; kwargs...) function newton(f, x0; xatol=nothing, xrtol=nothing, maxevals=100) + #Base.depwarn("`newton(f, x0)` is deprecated; use `find_zero(f, x0, Roots.Newton())` instead.", :newton) x = float(x0) T = typeof(x) atol = xatol !== nothing ? xatol : oneunit(T) * (eps(one(T)))^(4 / 5) @@ -430,126 +444,98 @@ function newton(f, x0; xatol=nothing, xrtol=nothing, maxevals=100) throw(ConvergenceFailed("No convergence")) end +## newton, halley, quadratic_inverse, superhalley, chebyshevlike +""" + Roots.newton(f, fp, x0; kwargs...) -## -------------------------------------------------- +Implementation of Newton's method: `xᵢ₊₁ = xᵢ - f(xᵢ)/f'(xᵢ)`. + +Arguments: + +* `f::Function` -- function to find zero of + +* `fp::Function` -- the derivative of `f`. + +* `x0::Number` -- initial guess. For Newton's method this may be complex. + +With the `ForwardDiff` package derivatives may be computed automatically. For example, defining +`D(f) = x -> ForwardDiff.derivative(f, float(x))` allows `D(f)` to be used for the first derivative. + +Keyword arguments are passed to `find_zero` using the `Roots.Newton()` method. + +See also `Roots.newton((f,fp), x0)` and `Roots.newton(fΔf, x0)` for simpler implementations. -## This is basically Order0(), but with different, default, tolerances employed -## It takes more function calls, but works harder to find exact zeros -## where exact means either iszero(fx), adjacent floats have sign change, or -## abs(fxn) <= 8 eps(xn) """ - dfree(f, xs) +function newton(f, fp, x0; kwargs...) + #Base.depwarn("`newton(f,fp, x0)` is deprecated; use `find_zero((f,fp), x0, Roots.Newton())` instead.", :newton) + find_zero((f, fp), x0, Newton(); kwargs...) +end -A more robust secant method implementation +## -------------------------------------------------- +#= +""" + Roots.halley(f, fp, fpp, x0; kwargs...) -Solve for `f(x) = 0` using an algorithm from *Personal Calculator Has Key -to Solve Any Equation f(x) = 0*, the SOLVE button from the -[HP-34C](http://www.hpl.hp.com/hpjournal/pdfs/IssuePDFs/1979-12.pdf). +Implementation of Halley's method (cf `?Roots.Halley()`). -This is also implemented as the `Order0` method for `find_zero`. +Arguments: -The initial values can be specified as a pair of two values, as in -`(a,b)` or `[a,b]`, or as a single value, in which case a value of `b` -is computed, possibly from `fb`. The basic idea is to follow the -secant method to convergence unless: +* `f::Function` -- function to find zero of -* a bracket is found, in which case `AlefeldPotraShi` is used; +* `fp::Function` -- derivative of `f`. -* the secant method is not converging, in which case a few steps of a - quadratic method are used to see if that improves matters. +* `fpp:Function` -- second derivative of `f`. -Convergence occurs when `f(m) == 0`, there is a sign change between -`m` and an adjacent floating point value, or `f(m) <= 2^3*eps(m)`. +* `x0::Number` -- initial guess -A value of `NaN` is returned if the algorithm takes too many steps -before identifying a zero. +With the `ForwardDiff` package derivatives may be computed automatically. For example, defining +`D(f) = x -> ForwardDiff.derivative(f, float(x))` allows `D(f)` and `D(D(f))` to be used for the first and second +derivatives, respectively. -# Examples +Keyword arguments are passed to `find_zero` using the `Roots.Halley()` method. -```julia -Roots.dfree(x -> x^5 - x - 1, 1.0) -``` +""" +=# +function halley(f, fp, fpp, x0; kwargs...) + #Base.depwarn("`halley(f,fp,fpp, x0)` is deprecated; use `find_zero((f, fp, fpp), x0, Roots.Halley())` instead.", :halley) + find_zero((f, fp, fpp), x0, Halley(); kwargs...) +end +#= """ -function dfree(f, xs) - if length(xs) == 1 - a = float(xs[1]) - fa = f(a) + Roots.quadratic_inverse(f, fp, fpp, x0; kwargs...) - h = eps(one(a))^(1 / 3) - da = h * oneunit(a) + abs(a) * h^2 # adjust for if eps(a) > h - b = float(a + da) - fb = f(b) - else - a, b = promote(float(xs[1]), float(xs[2])) - fa, fb = f(a), f(b) - end +Implementation of the quadratic inverse method (cf `?Roots.QuadraticInverse()`). - nan = (0 * a) / (0 * a) # try to preserve type - cnt, MAXCNT = 0, 5 * ceil(Int, -log(eps(one(a)))) # must be higher for BigFloat - MAXQUAD = 3 +Arguments: - if abs(fa) > abs(fb) - a, fa, b, fb = b, fb, a, fa - end +* `f::Function` -- function to find zero of - # we keep a, b, fa, fb, gamma, fgamma - quad_ctr = 0 - while !iszero(fb) - cnt += 1 +* `fp::Function` -- derivative of `f`. - if sign(fa) * sign(fb) < 0 - return solve(ZeroProblem(f, (a, b))) # faster than bisection(f, a, b) - end +* `fpp:Function` -- second derivative of `f`. - # take a secant step - gamma = float(b - (b - a) * fb / (fb - fa)) - # modify if gamma is too small or too big - if iszero(abs(gamma - b)) - gamma = b + 1 / 1000 * abs(b - a) # too small - elseif abs(gamma - b) >= 100 * abs(b - a) - gamma = b + sign(gamma - b) * 100 * abs(b - a) ## too big - end - fgamma = f(gamma) +* `x0::Number` -- initial guess - # change sign - if sign(fgamma) * sign(fb) < 0 - return solve(ZeroProblem(f, (gamma, b))) # faster than bisection(f, gamma, b) - end +With the `ForwardDiff` package derivatives may be computed automatically. For example, defining +`D(f) = x -> ForwardDiff.derivative(f, float(x))` allows `D(f)` and `D(D(f))` to be used for the first and second +derivatives, respectively. - # decreasing - if abs(fgamma) < abs(fb) - a, fa, b, fb = b, fb, gamma, fgamma - quad_ctr = 0 - cnt < MAXCNT && continue - end +Keyword arguments are passed to `find_zero` using the `Roots.QuadraticInverse()` method. - gamma = float(quad_vertex(a, fa, b, fb, gamma, fgamma)) - fgamma = f(gamma) - # decreasing now? - if abs(fgamma) < abs(fb) - a, fa, b, fb = b, fb, gamma, fgamma - quad_ctr = 0 - cnt < MAXCNT && continue - end +""" +=# +function quadratic_inverse(f, fp, fpp, x0; kwargs...) + #Base.depwarn("`quadratic_inverse(f,fp,fpp, x0)` is deprecated; use `find_zero((f, fp, fpp), x0, Roots.QuadraticInverse())` instead.", :quadratic_inverse) + find_zero((f, fp, fpp), x0, QuadraticInverse(); kwargs...) +end - quad_ctr += 1 - if (quad_ctr > MAXQUAD) || (cnt > MAXCNT) || iszero(gamma - b) || isnan(gamma) - bprev, bnext = prevfloat(b), nextfloat(b) - fbprev, fbnext = f(bprev), f(bnext) - sign(fb) * sign(fbprev) < 0 && return b - sign(fb) * sign(fbnext) < 0 && return b - for (u, fu) in ((b, fb), (bprev, fbprev), (bnext, fbnext)) - abs(fu) / oneunit(fu) <= 2^3 * eps(u / oneunit(u)) && return u - end - return nan # Failed. - end +function superhalley(f, fp, fpp, x0; kwargs...) + #Base.depwarn("`superhalley(f,fp,fpp, x0)` is deprecated; use `find_zero((f, fp, fpp), x0, Roots.SuperHalley())` instead.", :superhalley) + find_zero((f, fp, fpp), x0, SuperHalley(); kwargs...) +end - if abs(fgamma) < abs(fb) - b, fb, a, fa = gamma, fgamma, b, fb - else - a, fa = gamma, fgamma - end - end - b +function chebyshev_like(f, fp, fpp, x0; kwargs...) + #Base.depwarn("`chebyshev_like(f,fp,fpp, x0)` is deprecated; use `find_zero((f, fp, fpp), x0, Roots.ChebyshevLike())` instead.", :chebyshev_like) + find_zero((f, fp, fpp), x0, ChebyshevLike(); kwargs...) end diff --git a/src/find_zeros.jl b/src/find_zeros.jl index b686117d..7473dafe 100644 --- a/src/find_zeros.jl +++ b/src/find_zeros.jl @@ -2,6 +2,215 @@ # Algorithm due to @djsegal in https://github.com/JuliaMath/Roots.jl/pull/113 +## ----- +## Bisection +## +## Essentially from Jason Merrill https://gist.github.com/jwmerrill/9012954 +## cf. http://squishythinking.com/2014/02/22/bisecting-floats/ +## This also borrows a trick from https://discourse.julialang.org/t/simple-and-fast-bisection/14886 +## where we keep x1 so that y1 is negative, and x2 so that y2 is positive +## this allows the use of signbit over y1*y2 < 0 which avoid < and a multiplication +## this has a small, but noticeable impact on performance. +""" + bisection(f, a, b; [xatol, xrtol]) + +Performs bisection method to find a zero of a continuous +function. + +It is assumed that `(a,b)` is a bracket, that is, the function has +different signs at `a` and `b`. The interval `(a,b)` is converted to floating point +and shrunk when `a` or `b` is infinite. The function `f` may be infinite for +the typical case. If `f` is not continuous, the algorithm may find +jumping points over the x axis, not just zeros. + + +If non-trivial tolerances are specified, the process will terminate +when the bracket `(a,b)` satisfies `isapprox(a, b, atol=xatol, +rtol=xrtol)`. For zero tolerances, the default, for `Float64`, `Float32`, +or `Float16` values, the process will terminate at a value `x` with +`f(x)=0` or `f(x)*f(prevfloat(x)) < 0 ` or `f(x) * f(nextfloat(x)) < +0`. For other number types, the `Roots.A42` method is used. + +""" +function bisection(f, a::Number, b::Number; xatol=nothing, xrtol=nothing) + x1, x2 = adjust_bracket(float.((a, b))) + T = typeof(x1) + + atol = xatol === nothing ? zero(T) : abs(xatol) + rtol = xrtol === nothing ? zero(one(T)) : abs(xrtol) + CT = iszero(atol) && iszero(rtol) ? Val(:exact) : Val(:inexact) + + x1, x2 = float(x1), float(x2) + y1, y2 = f(x1), f(x2) + + _unitless(y1 * y2) >= 0 && error("the interval provided does not bracket a root") + + if isneg(y2) + x1, x2, y1, y2 = x2, x1, y2, y1 + end + + xm = Roots._middle(x1, x2) # for possibly mixed sign x1, x2 + ym = f(xm) + + while true + if has_converged(CT, x1, x2, xm, ym, atol, rtol) + return xm + end + + if isneg(ym) + x1, y1 = xm, ym + else + x2, y2 = xm, ym + end + + xm = Roots.__middle(x1, x2) + ym = f(xm) + end +end + +# -0.0 not returned by __middle, so isneg true on [-Inf, 0.0) +@inline isneg(x::T) where {T<:AbstractFloat} = signbit(x) +@inline isneg(x) = _unitless(x) < 0 + +@inline function has_converged(::Val{:exact}, x1, x2, m, ym, atol, rtol) + iszero(ym) && return true + isnan(ym) && return true + x1 != m && m != x2 && return false + return true +end + +@inline function has_converged(::Val{:inexact}, x1, x2, m, ym, atol, rtol) + iszero(ym) && return true + isnan(ym) && return true + val = abs(x1 - x2) <= atol + max(abs(x1), abs(x2)) * rtol + return val +end + +## ---- +## -------------------------------------------------- + +## This is basically Order0(), but with different, default, tolerances employed +## It takes more function calls, but works harder to find exact zeros +## where exact means either iszero(fx), adjacent floats have sign change, or +## abs(fxn) <= 8 eps(xn) +""" + dfree(f, xs) + +A more robust secant method implementation + +Solve for `f(x) = 0` using an algorithm from *Personal Calculator Has Key +to Solve Any Equation f(x) = 0*, the SOLVE button from the +[HP-34C](http://www.hpl.hp.com/hpjournal/pdfs/IssuePDFs/1979-12.pdf). + +This is also implemented as the `Order0` method for `find_zero`. + +The initial values can be specified as a pair of two values, as in +`(a,b)` or `[a,b]`, or as a single value, in which case a value of `b` +is computed, possibly from `fb`. The basic idea is to follow the +secant method to convergence unless: + +* a bracket is found, in which case `AlefeldPotraShi` is used; + +* the secant method is not converging, in which case a few steps of a + quadratic method are used to see if that improves matters. + +Convergence occurs when `f(m) == 0`, there is a sign change between +`m` and an adjacent floating point value, or `f(m) <= 2^3*eps(m)`. + +A value of `NaN` is returned if the algorithm takes too many steps +before identifying a zero. + +# Examples + +```julia +Roots.dfree(x -> x^5 - x - 1, 1.0) +``` + +""" +function dfree(f, xs) + if length(xs) == 1 + a = float(xs[1]) + fa = f(a) + + h = eps(one(a))^(1 / 3) + da = h * oneunit(a) + abs(a) * h^2 # adjust for if eps(a) > h + b = float(a + da) + fb = f(b) + else + a, b = promote(float(xs[1]), float(xs[2])) + fa, fb = f(a), f(b) + end + + nan = (0 * a) / (0 * a) # try to preserve type + cnt, MAXCNT = 0, 5 * ceil(Int, -log(eps(one(a)))) # must be higher for BigFloat + MAXQUAD = 3 + + if abs(fa) > abs(fb) + a, fa, b, fb = b, fb, a, fa + end + + # we keep a, b, fa, fb, gamma, fgamma + quad_ctr = 0 + while !iszero(fb) + cnt += 1 + + if sign(fa) * sign(fb) < 0 + return solve(ZeroProblem(f, (a, b))) # faster than bisection(f, a, b) + end + + # take a secant step + gamma = float(b - (b - a) * fb / (fb - fa)) + # modify if gamma is too small or too big + if iszero(abs(gamma - b)) + gamma = b + 1 / 1000 * abs(b - a) # too small + elseif abs(gamma - b) >= 100 * abs(b - a) + gamma = b + sign(gamma - b) * 100 * abs(b - a) ## too big + end + fgamma = f(gamma) + + # change sign + if sign(fgamma) * sign(fb) < 0 + return solve(ZeroProblem(f, (gamma, b))) # faster than bisection(f, gamma, b) + end + + # decreasing + if abs(fgamma) < abs(fb) + a, fa, b, fb = b, fb, gamma, fgamma + quad_ctr = 0 + cnt < MAXCNT && continue + end + + gamma = float(quad_vertex(a, fa, b, fb, gamma, fgamma)) + fgamma = f(gamma) + # decreasing now? + if abs(fgamma) < abs(fb) + a, fa, b, fb = b, fb, gamma, fgamma + quad_ctr = 0 + cnt < MAXCNT && continue + end + + quad_ctr += 1 + if (quad_ctr > MAXQUAD) || (cnt > MAXCNT) || iszero(gamma - b) || isnan(gamma) + bprev, bnext = prevfloat(b), nextfloat(b) + fbprev, fbnext = f(bprev), f(bnext) + sign(fb) * sign(fbprev) < 0 && return b + sign(fb) * sign(fbnext) < 0 && return b + for (u, fu) in ((b, fb), (bprev, fbprev), (bnext, fbnext)) + abs(fu) / oneunit(fu) <= 2^3 * eps(u / oneunit(u)) && return u + end + return nan # Failed. + end + + if abs(fgamma) < abs(fb) + b, fb, a, fa = gamma, fgamma, b, fb + else + a, fa = gamma, fgamma + end + end + b +end + + # A naive approach to find zeros: split (a,b) by n points, look into each for a zero # * k is oversampling rate for bisection. (It is relatively cheap to check for a bracket so we # oversample our intervals looking for brackets diff --git a/src/trace.jl b/src/trace.jl index a2d2f852..667a53b7 100644 --- a/src/trace.jl +++ b/src/trace.jl @@ -32,13 +32,13 @@ julia> using Roots julia> f(x) = x^2-2 f (generic function with 1 method) -julia> tracker = Roots.Tracks() +julia> tracks = Roots.Tracks() Algorithm has not been run -julia> find_zero(f, (0, 2), Roots.Secant(), tracks=tracker) ≈ √2 +julia> find_zero(f, (0, 2), Roots.Secant(); tracks) ≈ √2 # keyword match true -julia> tracker +julia> tracks Results of univariate zero finding: * Converged to: 1.4142135623730947 @@ -58,13 +58,13 @@ x₇ = 1.4142114384748701 fx₇ = -6.0072868388605372e-06 x₈ = 1.4142135626888697 fx₈ = 8.9314555751229818e-10 x₉ = 1.4142135623730947 fx₉ = -8.8817841970012523e-16 -julia> tracker = Roots.Tracks() +julia> tracks = Roots.Tracks() Algorithm has not been run -julia> find_zero(sin, (3, 4), Roots.A42(), tracks=tracker) ≈ π +julia> find_zero(sin, (3, 4), Roots.A42(); tracks) ≈ π true -julia> tracker +julia> tracks Results of univariate zero finding: * Converged to: 3.141592653589793 diff --git a/test/runtests.jl b/test/runtests.jl index ec671556..4bc8e8e2 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -23,20 +23,23 @@ avg(x) = sum(x) / length(x) include("./test_find_zero.jl") include("./test_bracketing.jl") include("./test_derivative_free.jl") -include("./test_simple.jl") include("./test_find_zeros.jl") -include("./test_fzero.jl") include("./test_newton.jl") include("./test_chain_rules.jl") -include("./test_simple.jl") include("./test_composable.jl") VERSION >= v"1.6.0" && include("./test_allocations.jl") VERSION >= v"1.9.0" && include("./test_extensions.jl") VERSION >= v"1.12.0" && include("./test_jet.jl") +@testset "Aqua" begin + Aqua.test_all(Roots) +end #include("./runbenchmarks.jl") #include("./benchmark-bracketing-methods.jl") #include("./test_derivative_free_interactive.jl") -Aqua.test_all(Roots) + +# deprecated +include("./test_fzero.jl") +include("./test_simple.jl") diff --git a/test/test_allocations.jl b/test/test_allocations.jl index 5f2b8472..57454fa4 100644 --- a/test/test_allocations.jl +++ b/test/test_allocations.jl @@ -20,6 +20,10 @@ import BenchmarkTools Roots.Brent(), Roots.Ridders(), Roots.ITP(), + Roots.RegulaFalsi(:AndersonBjork), + Roots.RegulaFalsi(:Illinois), + Roots.RegulaFalsi(:Ford3), + Roots.FalsePosition(), ) # not FalsePosition() Ns = (Roots.Newton(), Roots.Halley(), Roots.Schroder()) for M in Ms @@ -102,11 +106,3 @@ import BenchmarkTools @test @inferred solve(ZeroProblem(sin, 3.0f0)) ≈ pi end end - -@testset "simple: zero allocations" begin - @test BenchmarkTools.@ballocated(Roots.bisection(sin, 3, 4)) == 0 - @test BenchmarkTools.@ballocated(Roots.secant_method(sin, 3)) == 0 - @test BenchmarkTools.@ballocated(Roots.muller(sin, 2.9, 3.0, 3.1)) == 0 - @test BenchmarkTools.@ballocated(Roots.newton((sin, cos), 3)) == 0 - @test BenchmarkTools.@ballocated(Roots.dfree(sin, 3)) == 0 -end diff --git a/test/test_bracketing.jl b/test/test_bracketing.jl index 746ac578..a0fb987e 100644 --- a/test/test_bracketing.jl +++ b/test/test_bracketing.jl @@ -278,44 +278,35 @@ end # end @testset "bracketing methods" begin - - ## Test for failures, ideally all of these would be 0 - ## test for residual, ideally small - ## test for evaluation counts, ideally not so low for these problems - - ## exact_bracket - Ms = [ - Roots.Brent(), + Ms = ( Roots.A42(), Roots.AlefeldPotraShi(), Roots.Chandrapatla(), - Roots.ITP(), - Roots.Ridders(), Roots.Bisection(), Roots.ModAB(), - ] - results = [run_tests((f, b) -> find_zero(f, b, M), name="$M") for M in Ms] - maxfailures = maximum(length(result.failures) for result in results) - maxresidual = maximum(result.maxresidual for result in results) - cnts = [result.evalcount for result in results] - - @test maxfailures == 0 - @test maxresidual <= 5e-13 - @test avg(cnts) <= 4700 + Roots.Brent(), + Roots.Ridders(), + Roots.ITP(), + [Roots.FalsePosition(i) for i in 1:12]..., + [ + Roots.RegulaFalsi(m) for + m in (:Illinois, :Pegasus, :AndersonBjork, :Ford3, :Ford4) + ]..., + ) - ## False position has larger residuals - ## Fn #13 fails on numbers 2 and 4 until maxsteps is increased; 100 works - Ms = [Roots.FalsePosition(i) for i in 1:12] + ## Test for failures, ideally all of these would be 0 + ## test for residual, ideally small + ## test for evaluation counts, ideally not so low for these problems results = [run_tests((f, b) -> find_zero(f, b, M), name="$M") for M in Ms] - maxfailures = maximum(length(result.failures) for result in results) - maxresidual = maximum(result.maxresidual for result in results) + failures = [length(result.failures) for result in results] + residuals = [result.maxresidual for result in results] cnts = [result.evalcount for result in results] - @test maxfailures <= 1 - @test maxresidual <= 1e-5 - @test avg(cnts) <= 3000 + @test maximum(failures) == 0 + @test maximum(residuals) <= 5e-14 + @test avg(cnts) <= 4000 - ## issue 412 check for bracket + ## issue 412 check for bracket in bracketing methods for M in Ms @test_throws ArgumentError find_zero(x -> x - 1, (-3, 0), M) @test_throws ArgumentError find_zero(x -> 1 + x^2, (10, 20), M) @@ -350,12 +341,25 @@ end ((x, n=20) -> x^2 + sin(x / n) - 1 / 4, [0, 1]), ] - for (fn_, ab) in galadino_probs - for M in (FalsePosition(i) for i in 1:12) + for (no, (fn_, ab)) in enumerate(galadino_probs) + for i in vcat(1, 3:12) # issue with 2 not converging + M = FalsePosition(i) g = Cnt(fn_) x0_ = find_zero(g, ab, M) @test abs(fn_(x0_)) <= 1e-7 - @test g.cnt <= 50 + @test g.cnt <= 60 + end + end + for (fn_, ab) in galadino_probs + Ms = [ + Roots.RegulaFalsi(m) for + m in (:Illinois, :Pegasus, :AndersonBjork, :Ford3, :Ford4) + ] + for M in Ms + g = Cnt(fn_) + x0_ = find_zero(g, ab, M) + @test abs(fn_(x0_)) <= 1e-14 + @test g.cnt <= 60 end end end @@ -377,7 +381,7 @@ end @test find_zero(f, (-1, 1), M, p=eps()) ≈ eps() atol = 2eps() end - @test iszero(@inferred(find_zero(f, (-1, 1), Roots.Bisection()))) + @test iszero(@inferred(find_zero(f, (-1, 1), Roots.Bisection()))) # too fiddly to get zero, might get eps(0.0) # XXX changes with relaxed tolerance (adding non-zero xatol) #@test_throws Roots.ConvergenceFailed find_zero(f, (-1, 1), Roots.A42()) #@test_throws Roots.ConvergenceFailed find_zero(f, (-1, 1), Roots.AlefeldPotraShi()) @@ -434,8 +438,4 @@ end x = find_zero(f, (xleft, xright)) @test abs(f(x)) <= 2eps(BigFloat) - # simple a42() - m = run_tests(Roots.a42) - VERSION >= v"1.6" && @test isempty(m.failures) - @test m.evalcount <= 3000 # paper says 2884, this has 2877 end diff --git a/test/test_composable.jl b/test/test_composable.jl index 326e0e3a..0aca15c6 100644 --- a/test/test_composable.jl +++ b/test/test_composable.jl @@ -91,9 +91,10 @@ using Measurements c = measurement(-10.0, 0.1) f(x) = a * x^2 + b * x + c x₀ = (measurement(-3.0, 0.1), measurement(0.0, 0.1)) - for M in (A42(), AlefeldPotraShi(), Bisection()) + for M in (Secant(), Bisection()) @test find_zero(f, x₀, M) ≈ -2.0 end + @test find_zero(f, measurement(0.0, 0.1), Secant()) ≈ -2.0 @test find_zero(f, measurement(0.0, 0.1)) ≈ -2.0 end end diff --git a/test/test_derivative_free_interactive.jl b/test/test_derivative_free_interactive.jl index 77615834..6df3af8e 100644 --- a/test/test_derivative_free_interactive.jl +++ b/test/test_derivative_free_interactive.jl @@ -4,7 +4,7 @@ ### Benchmarking tests ## We have ## -## * visualize_diagonostics(which): to see summaries of the methods +## * visualize_diagnostics(which): to see summaries of the methods ## over the different functions with which in (:summary, :counts, ## :residuals) ## @@ -102,7 +102,7 @@ function vvta1(vs, T) end ## Return Dict of arrays -function create_diagonostics(Ms, Fs, nms) +function create_diagnostics(Ms, Fs, nms) @assert length(nms) == length(Ms) out = Array{Any}(undef, length(Ms), length(Fs)) @@ -367,7 +367,7 @@ end ## Main interface for interactive use fname = joinpath(@__DIR__, "derivative_free_diagnostics.json") elide_ascii(x, n=12) = length(x) > n ? x[1:(n - 3)] * "..." * x[(end - 1):end] : x -function create_diagonostics() +function create_diagnostics() meths = [ Order0(), Order1(), @@ -388,12 +388,12 @@ function create_diagonostics() Ms = [(f, b) -> find_zero(f, b, M) for M in meths] # F(f,b), name nms = elide_ascii.([replace(string(M), r"^Roots." => "") for M in meths]) Fs = known_functions - create_diagonostics(Ms, Fs, nms) + create_diagnostics(Ms, Fs, nms) end ## write out current diagnostic test function write_out() - println("Creating diagonostics to save") - write_out(fname, create_diagonostics()) + println("Creating diagnostics to save") + write_out(fname, create_diagnostics()) end ## visualize state @@ -404,7 +404,7 @@ Show diagnostics summary `which` is one of `(:all, :summary, :counts, :residuals)` """ -visualize_diagnostics(which=:summary) = visualize_diagnostics(create_diagonostics(), which) +visualize_diagnostics(which=:summary) = visualize_diagnostics(create_diagnostics(), which) ## identify regressions from currently saved state """ @@ -420,7 +420,7 @@ function identify_regressions() return String[] # empty end - Dnew = create_diagonostics() + Dnew = create_diagnostics() Dold = read_in(fname) out = identify_regressions(Dnew, Dold) @@ -430,7 +430,7 @@ end ## tests for newton, halley import ForwardDiff: derivative D(f, n=1) = n > 1 ? D(D(f), n - 1) : x -> derivative(f, float(x)) -function derivative_based_diagonostics() +function derivative_based_diagnostics() Ms = ( (f, b) -> Roots.find_zero((f, D(f)), b, Roots.Newton()), (f, b) -> Roots.find_zero((f, D(f), D(f, 2)), b, Roots.Halley()), @@ -439,7 +439,7 @@ function derivative_based_diagonostics() ) nms = ("Newton", "Halley", "Schroder", "Order5") Fs = known_functions - create_diagonostics(Ms, Fs, nms) + create_diagnostics(Ms, Fs, nms) end ## Order of convergence diff --git a/test/test_find_zero.jl b/test/test_find_zero.jl index d57597f7..39dda506 100644 --- a/test/test_find_zero.jl +++ b/test/test_find_zero.jl @@ -429,9 +429,6 @@ end @test @inferred(find_zero(x -> x + 1, (-1, 1))) == -1 end - ## issue #178 passing through method - @test fzero(sin, 3, 4, Roots.Brent()) ≈ π - ## issue #188 with A42 f = let a = 0.18 x -> x * (1 - x^2) / ((x^2 + a^2) * (1 + a^2 * x^2)) diff --git a/test/test_find_zeros.jl b/test/test_find_zeros.jl index a8b87b4f..81b0ef32 100644 --- a/test/test_find_zeros.jl +++ b/test/test_find_zeros.jl @@ -2,6 +2,7 @@ using Roots using Test +using BenchmarkTools mutable struct CallableFunction f @@ -136,3 +137,26 @@ end @test eltype(rts) == T end end + +@testset "dfree/bisection" begin + + # dfree + fpoly = x -> x^5 - x - 1 + xrt = Roots.dfree(fpoly, 1.0) + @test abs(fpoly(xrt)) <= 1e-14 + + @test BenchmarkTools.@ballocated(Roots.dfree(sin, 3)) == 0 + + # bisection + # bisection + xrt = Roots.bisection(sin, 3.0, 4.0) + @test isapprox(xrt, pi) + + xrt = Roots.bisection(sin, 3.0, 4.0, xatol=1e-3) + @test abs(sin(xrt)) >= 1e-7 # not too close + + xrt = Roots.bisection(sin, big(3.0), big(4.0)) + @test isapprox(xrt, pi) + @test BenchmarkTools.@ballocated(Roots.bisection(sin, 3, 4)) == 0 + +end diff --git a/test/test_fzero.jl b/test/test_fzero.jl index 4d1d0c21..778235f7 100644 --- a/test/test_fzero.jl +++ b/test/test_fzero.jl @@ -67,4 +67,8 @@ import Roots.fzero ### issue with fzeros and roots near 'b' @test 0 < maximum(fzeros(x -> sin(x) - 1 / 1000 * x, 0, pi)) < pi + + ## issue #178 passing through method + @test fzero(sin, 3, 4, Roots.Brent()) ≈ π + end diff --git a/test/test_newton.jl b/test/test_newton.jl index 89d31253..aae5fe8a 100644 --- a/test/test_newton.jl +++ b/test/test_newton.jl @@ -1,19 +1,7 @@ using Test -import Roots.newton, - Roots.halley, Roots.superhalley, Roots.quadratic_inverse, Roots.chebyshev_like +using Roots @testset "Test Newton, Halley, Schroder methods" begin - @test abs(newton(sin, cos, 0.5) - 0.0) <= 100 * eps(1.0) - @test newton(cos, x -> -sin(x), 1.0) ≈ pi / 2 - @test newton(x -> x^2 - 2x - 1, x -> 2x - 2, 3.0) ≈ 2.414213562373095 - @test abs(newton(x -> exp(x) - cos(x), x -> exp(x) + sin(x), 3.0) - 0.0) <= 1e-14 - @test halley(x -> x^2 - 2x - 1, x -> 2x - 2, x -> 2, 3.0) ≈ 2.414213562373095 - @test quadratic_inverse(x -> x^2 - 2x - 1, x -> 2x - 2, x -> 2, 3.0) ≈ 2.414213562373095 - @test superhalley(x -> x^2 - 2x - 1, x -> 2x - 2, x -> 2, 3.0) ≈ 2.414213562373095 - @test chebyshev_like(x -> x^2 - 2x - 1, x -> 2x - 2, x -> 2, 3.0) ≈ 2.414213562373095 - a = halley(x -> exp(x) - cos(x), x -> exp(x) + sin(x), x -> exp(x) + cos(x), 3.0) - @test abs(a - 0.0) <= 1e-14 - # find_zero calls @test find_zero((x -> x^2 - 2x - 1, x -> 2x - 2), 3.0, Roots.Newton()) ≈ 2.414213562373095 @@ -31,18 +19,8 @@ import Roots.newton, @test find_zero((x -> x^2 - 2x - 1, x -> 2x - 2, x -> 2), 3.0, Roots.Schroder()) ≈ 2.414213562373095 - ## test with Complex input - - @test real(Roots.newton(x -> x^3 - 1, x -> 3x^2, 1 + im)) ≈ 1.0 - @test real(Roots.newton(x -> x^3 - 1, x -> 3x^2, 1 + 10im)) ≈ (-1 / 2) - - ## Issue #143 test with new interface - Roots.newton(sin, cos, 3.0) ≈ π # uses find_zero - Roots.newton((sin, cos), 3.0) ≈ π # uses simple - fdf = x -> (sin(x), sin(x) / cos(x)) # (f, f/f') @test Roots.find_zero(fdf, 3.0, Roots.Newton()) ≈ π # uses find_zero - Roots.newton(fdf, 3.0) ≈ π # uses simple fdfdf = x -> (sin(x), sin(x) / cos(x), -cos(x) / sin(x), sin(x) / cos(x)) # (f, f/f', f'/f'', f''/f''') @test Roots.find_zero(fdfdf, 3.0, Roots.Halley()) ≈ π @@ -73,7 +51,7 @@ import Roots.newton, Roots.Bisection(), ) ≈ sqrt(2) - @test_throws Roots.ConvergenceFailed Roots.newton((x -> x^2 + 1, x -> 2x), 0) + end @testset "bracketed Halley" begin diff --git a/test/test_simple.jl b/test/test_simple.jl index 84a759f4..8db8ee3f 100644 --- a/test/test_simple.jl +++ b/test/test_simple.jl @@ -1,18 +1,9 @@ using Test using BenchmarkTools +import Roots: newton, halley, superhalley, quadratic_inverse, chebyshev_like @testset "simpler implementations" begin - # bisection - xrt = Roots.bisection(sin, 3.0, 4.0) - @test isapprox(xrt, pi) - - xrt = Roots.bisection(sin, 3.0, 4.0, xatol=1e-3) - @test abs(sin(xrt)) >= 1e-7 # not to0 close - - xrt = Roots.bisection(sin, big(3.0), big(4.0)) - @test isapprox(xrt, pi) - # secant_method fpoly = x -> x^5 - x - 1 xrt = Roots.secant_method(fpoly, 1.0) @@ -37,13 +28,49 @@ using BenchmarkTools @test Roots.muller(expoly, -0.7 - 0.5im) ≈ -1.0 - # dfree - fpoly = x -> x^5 - x - 1 - xrt = Roots.dfree(fpoly, 1.0) - @test abs(fpoly(xrt)) <= 1e-14 - - # newton + # newton/halley/superhalley/quadratic_inverse/chebyshev_like @test Roots.newton((sin, cos), 3.0) ≈ pi u = Roots.newton(x -> (sin(x), sin(x) / cos(x)), 3.0, xatol=1e-10, xrtol=1e-10) @test abs(u - pi) <= 1e-8 + + @test abs(newton(sin, cos, 0.5) - 0.0) <= 100 * eps(1.0) + @test newton(cos, x -> -sin(x), 1.0) ≈ pi / 2 + @test newton(x -> x^2 - 2x - 1, x -> 2x - 2, 3.0) ≈ 2.414213562373095 + @test abs(newton(x -> exp(x) - cos(x), x -> exp(x) + sin(x), 3.0) - 0.0) <= 1e-14 + @test halley(x -> x^2 - 2x - 1, x -> 2x - 2, x -> 2, 3.0) ≈ 2.414213562373095 + @test quadratic_inverse(x -> x^2 - 2x - 1, x -> 2x - 2, x -> 2, 3.0) ≈ 2.414213562373095 + @test superhalley(x -> x^2 - 2x - 1, x -> 2x - 2, x -> 2, 3.0) ≈ 2.414213562373095 + @test chebyshev_like(x -> x^2 - 2x - 1, x -> 2x - 2, x -> 2, 3.0) ≈ 2.414213562373095 + a = halley(x -> exp(x) - cos(x), x -> exp(x) + sin(x), x -> exp(x) + cos(x), 3.0) + @test abs(a - 0.0) <= 1e-14 + @test_throws Roots.ConvergenceFailed Roots.newton((x -> x^2 + 1, x -> 2x), 0) + + ## test with Complex input + @test real(Roots.newton(x -> x^3 - 1, x -> 3x^2, 1 + im)) ≈ 1.0 + @test real(Roots.newton(x -> x^3 - 1, x -> 3x^2, 1 + 10im)) ≈ (-1 / 2) + + ## Issue #143 test with new interface + @test Roots.newton(sin, cos, 3.0) ≈ π # uses find_zero + @test Roots.newton((sin, cos), 3.0) ≈ π # uses simple + + fdf = x -> (sin(x), sin(x) / cos(x)) # (f, f/f') + @test Roots.find_zero(fdf, 3.0, Roots.Newton()) ≈ π # uses find_zero + @test Roots.newton(fdf, 3.0) ≈ π # uses simple + + + + + # a42 + # simple a42() + #m = run_tests(Roots.a42) # in test_bracketing + #VERSION >= v"1.6" && @test isempty(m.failures) + #@test m.evalcount <= 3000 # paper says 2884, this has 2877 + +end + + +@testset "simple: zero allocations" begin + @test BenchmarkTools.@ballocated(Roots.secant_method(sin, 3)) == 0 + @test BenchmarkTools.@ballocated(Roots.muller(sin, 2.9, 3.0, 3.1)) == 0 + @test BenchmarkTools.@ballocated(Roots.newton((sin, cos), 3)) == 0 end