IntervalRootFinding and Legendre polynomials: error

General Usage
1

When I try to find the roots of the 26th degree Legendre polynomial with IntervalRootFinding.jl, I get an error I can’t decode.

julia> using IntervalArithmetic

julia> using LegendrePolynomials

julia> using IntervalRootFinding

julia> f(x)=(x-2)*(x+3)
f (generic function with 1 method)

julia> roots(f,-4..3)
2-element Vector{Root{Interval{Float64}}}:
 Root([1.99999, 2.00001], :unique)
 Root([-3.00001, -2.99999], :unique)

So far, so good.

But if I define:

julia> g(x) = Pl(x,26)
g (generic function with 1 method)

julia> roots(g,-1..1)
ERROR: TypeError: non-boolean (Interval{Float64}) used in boolean context
Stacktrace:
  [1] _abs_deriv
    @ C:\Users\amgough\.julia\packages\DiffRules\Hgr7t\src\rules.jl:72 [inlined]
  [2] abs
    @ C:\Users\amgough\.julia\packages\ForwardDiff\QOqCN\src\dual.jl:203 [inlined]
  [3] checkdomain(x::ForwardDiff.Dual{ForwardDiff.Tag{typeof(g), Interval{Float64}}, Interval{Float64}, 1})
    @ LegendrePolynomials C:\Users\amgough\.julia\packages\LegendrePolynomials\jbIWy\src\LegendrePolynomials.jl:13
  [4] LegendrePolynomialIterator
    @ C:\Users\amgough\.julia\packages\LegendrePolynomials\jbIWy\src\LegendrePolynomials.jl:92 [inlined]
  [5] LegendrePolynomialIterator
    @ C:\Users\amgough\.julia\packages\LegendrePolynomials\jbIWy\src\LegendrePolynomials.jl:96 [inlined]
  [6] Pl(x::ForwardDiff.Dual{ForwardDiff.Tag{typeof(g), Interval{Float64}}, Interval{Float64}, 1}, l::Int64)
    @ LegendrePolynomials C:\Users\amgough\.julia\packages\LegendrePolynomials\jbIWy\src\LegendrePolynomials.jl:149
  [7] g
    @ .\REPL[326]:1 [inlined]
  [8] derivative
    @ C:\Users\amgough\.julia\packages\ForwardDiff\QOqCN\src\derivative.jl:14 [inlined]
  [9] #41
    @ C:\Users\amgough\.julia\packages\IntervalRootFinding\8d8pn\src\roots.jl:162 [inlined]
 [10] � (f::typeof(g), f′::IntervalRootFinding.var"#41#42"{typeof(g)}, X::Interval{Float64}, α::Float64)
    @ IntervalRootFinding C:\Users\amgough\.julia\packages\IntervalRootFinding\8d8pn\src\contractors.jl:14
 [11] #35
    @ C:\Users\amgough\.julia\packages\IntervalRootFinding\8d8pn\src\contractors.jl:115 [inlined]
 [12] determine_region_status(op::IntervalRootFinding.var"#35#36"{Newton{typeof(g), IntervalRootFinding.var"#41#42"{typeof(g)}}, F
loat64}, f::typeof(g), R::Root{Interval{Float64}})
    @ IntervalRootFinding C:\Users\amgough\.julia\packages\IntervalRootFinding\8d8pn\src\contractors.jl:153
 [13] Newton
    @ C:\Users\amgough\.julia\packages\IntervalRootFinding\8d8pn\src\contractors.jl:116 [inlined]
 [14] Newton
    @ C:\Users\amgough\.julia\packages\IntervalRootFinding\8d8pn\src\contractors.jl:115 [inlined]
 [15] process(search::DepthFirstSearch{Interval{Float64}, Newton{typeof(g), IntervalRootFinding.var"#41#42"{typeof(g)}}, Float64},
 r::Root{Interval{Float64}})
    @ IntervalRootFinding C:\Users\amgough\.julia\packages\IntervalRootFinding\8d8pn\src\roots.jl:56
 [16] iterate(search::DepthFirstSearch{Interval{Float64}, Newton{typeof(g), IntervalRootFinding.var"#41#42"{typeof(g)}}, Float64},
 wt::IntervalRootFinding.BBTree{Root{Interval{Float64}}})
    @ IntervalRootFinding C:\Users\amgough\.julia\packages\IntervalRootFinding\8d8pn\src\branch_and_bound.jl:307
 [17] iterate
    @ C:\Users\amgough\.julia\packages\IntervalRootFinding\8d8pn\src\branch_and_bound.jl:303 [inlined]
 [18] branch_and_prune(r::Root{Interval{Float64}}, contractor::Newton{typeof(g), IntervalRootFinding.var"#41#42"{typeof(g)}}, sear
ch::Type{DepthFirstSearch}, tol::Float64)
    @ IntervalRootFinding C:\Users\amgough\.julia\packages\IntervalRootFinding\8d8pn\src\roots.jl:86
 [19] _roots(f::Function, deriv::Function, r::Root{Interval{Float64}}, contractor::Type{Newton}, strategy::Type{DepthFirstSearch},
 tol::Float64)
    @ IntervalRootFinding C:\Users\amgough\.julia\packages\IntervalRootFinding\8d8pn\src\roots.jl:169
 [20] _roots
    @ C:\Users\amgough\.julia\packages\IntervalRootFinding\8d8pn\src\roots.jl:163 [inlined]
 [21] _roots
    @ C:\Users\amgough\.julia\packages\IntervalRootFinding\8d8pn\src\roots.jl:192 [inlined]
 [22] roots (repeats 2 times)
    @ C:\Users\amgough\.julia\packages\IntervalRootFinding\8d8pn\src\roots.jl:123 [inlined]
 [23] top-level scope
    @ REPL[327]:1

Question: Which package is throwing the error?

Maybe it doesn’t matter. For if I define the 26th Legendre polynomial explicitly, roots() takes forever or never finishes:

julia> h(x) = (61989816618513x^26 - 395033145117975x^24 + 1112542327066950x^22 - 1822675727322450x^20 + 1923935489951475x^18 - 136
9126185872445x^16 + 667866432132900x^14 - 222622144044300x^12 + 49638721307175x^10 - 7091245901025x^8 + 601681470390x^6 - 26466926
850x^4 + 456326325x^2 - 1300075)/8388608
h (generic function with 1 method)

julia> roots(h,-1..1)

It’s been a few hours now, and no result.

Thanks for your help.

Environment info:

Julia Version 1.6.2
Commit 1b93d53fc4 (2021-07-14 15:36 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Core(TM) i7-4960X CPU @ 3.60GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-11.0.1 (ORCJIT, ivybridge)
Environment:
  JULIA.EXECUTABLEPATH = "C:\Users\amgough\Julia\bin\julia.exe"
  JULIA_GLPK_LIBRARY_PATH = C:\Program Files\glpk
  JULIA_HOME = C:\Users\amgough\Julia-1.6.2\bin\julia.exe
  JULIA_NUM_THREADS = 4
  JULIA_EDITOR = code
2

I believe the problem lies in the fact that IntervalRootFinding.jl uses the Newton method by default, which also needs to evaluate the derivative of the function. If you don’t provide it explicitly, this will be calculated by ForwardDiff. However, as the stacktrace indicates an error occurs during this step when evaluating the derivative of g(x), apparently connected to a call to abs somewhere in the function.

I get the same error using a simpler example involving the absolute value:

julia> using IntervalArithmetic, IntervalRootFinding

julia> roots(x -> abs(x)-0.5, -1..1)
ERROR: TypeError: non-boolean (Interval{Float64}) used in boolean context
Stacktrace:
  [1] _abs_deriv
    @ ~/.julia/packages/DiffRules/Hgr7t/src/rules.jl:72 [inlined]

...

I am not too familiar with the internals of ForwardDiff and DiffRules, but I think the following issue might be related: https://github.com/JuliaIntervals/IntervalRootFinding.jl/issues/98.

To circumvent this, it is probably easiest to explicitly provide the derivative yourself. For your Legendre polynomial example, this is conveniently already provided by LegendrePolynomials.jl (dnPl(x, l, n) gives the n th derivative of the Legendre polynomial of degree l). Using a second order polynomial as an example:

julia> using LegendrePolynomials

julia> roots(x -> Pl(x, 2), x -> dnPl(x, 2, 1), -1..1)
2-element Vector{Root{Interval{Float64}}}:
 Root([0.57735, 0.577351], :unique)
 Root([-0.577351, -0.57735], :unique)

This also works (in the sense that no error is thrown) for higher orders, although that is of course a much harder problem and at some point quite a few intervals with status “:unknown” appear in the output.

Alternatively, since the Bisection method doesn’t need a derivative, the following will also not error:

julia> roots(x -> Pl(x, 2), -1..1, Bisection)
2-element Vector{Root{Interval{Float64}}}:
 Root([-0.577351, -0.57735], :unknown)
 Root([0.57735, 0.577351], :unknown)

For Pl(x, 26) it does indeed seem to take quite a long time with the Newton method, switching to the Krawczyk method terminates much faster for me (output truncated):

julia> roots(x -> Pl(x, 26), x -> dnPl(x, 26, 1), -1..1, Krawczyk)
211-element Vector{Root{Interval{Float64}}}:
 Root([0.995885, 0.995886], :unknown)
 Root([-0.978357, -0.978356], :unknown)
 Root([0.947159, 0.94716], :unknown)
 Root([-0.947159, -0.947158], :unknown)
 Root([0.902637, 0.902638], :unknown)
 Root([0.902612, 0.902613], :unknown)
 Root([0.845445, 0.845446], :unknown)
 Root([-0.845446, -0.845445], :unknown)
 Root([-0.978378, -0.978377], :unknown)
 Root([-0.292005, -0.292004], :unique)
 
 ...
1 Like