# IntervalRootFinding and Legendre polynomials: error

**URL:** <https://discourse.julialang.org/t/intervalrootfinding-and-legendre-polynomials-error/65541>\
**Category:** General Usage\
**Tags:** roots\
**Created:** [July 30, 2021, 6:30am UTC](https://discourse.julialang.org/t/intervalrootfinding-and-legendre-polynomials-error/65541 "2021-07-30T06:30:47Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![blackeneth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/blackeneth/32/10353_2.png) [@blackeneth](https://discourse.julialang.org/u/blackeneth)\
**Post date:** [July 30, 2021, 6:30am UTC](https://discourse.julialang.org/t/intervalrootfinding-and-legendre-polynomials-error/65541/1 "2021-07-30T06:30:48Z")

</div>

When I try to find the roots of the 26th degree Legendre polynomial with [IntervalRootFinding.jl](https://github.com/JuliaIntervals/IntervalRootFinding.jl/issues), I get an error I can’t decode.

```julia
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
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
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
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

```

---

<div class="post-metadata">

**Author:** ![bvrb](https://avatars.discourse-cdn.com/v4/letter/b/76d3ee/32.png) [@bvrb](https://discourse.julialang.org/u/bvrb)\
**Post date:** [July 31, 2021, 9:54am UTC](https://discourse.julialang.org/t/intervalrootfinding-and-legendre-polynomials-error/65541/2 "2021-07-31T09:54:15Z")

</div>

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
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](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
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
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
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)
 
 ...

```