How to Return the Roots from IntervalRootFinding.jl in pi terms?

Hi all,

I want to know how to make the roots computed by IntervalRootFinding returned in \pi terms instead of decimals / float value / numeric term ?

using SymPy, IntervalArithmetic, IntervalRootFinding
@syms x

v1(x) = log(1.5 + sin(x))
v2(x) = cos(x)/(1.5 + sin(x))

vd = diff((v1(x)), (x))
vdd = simplify(diff((vd), (x)))
V1 = IntervalRootFinding.roots(v2, -5..5)

println("Computing with Julia")
println("f'(x) = ", vd)
println("f''(x) = ", vdd)

println("a.  absolute extreme points: ")
V1

This might be related to the difference between symbolic expressions and Float64 numbers. Floating point numbers are just an approximation to most “real” values that we might be interested in. That’s usually fine, because the approximation is very good, but if we want the “true” mathematical value, like in your case, the distinction is important.

It looks like IntervalRootFinding just finds the roots of normal Julia functions using floats. So even if the result is actually \pi/2 we just get the closest float that it could find and we can’t really tell that the solution is \pi / 2.

However, symbolic libraries like SymPy are exactly made for this kind of scenario and provides some tools to find roots (it’s usually just called solve since we want to solve an equation f(x) = 0.

In this example:

julia> f1(x) = cos(x)/(1.5 + sin(x));
julia> solve(f1(x), x)      # solving the equation f1(x) == 0 for x
2-element Vector{Sym}:
 1.57079632679490
 4.71238898038469

Hmm, it seems like it still didn’t work. The result are now symbolic numbers (Vector{Sym} not Vector{Float64}), but that has a different, slightly subtle reason: The 1.5 is a Julia floating point number and will be converted to a SymPy symbol when calling the function.

julia> Sym(1.5)
1.50000000000000

Although it is the right value, SymPy has some difficulty realizing that it is actually 3/2. I’m not exactly sure why this matters here (because the numerator cos(x) is zero regardless of the denominator), but the following would work:

julia> f2(x) = cos(x)/(Sym(3) / Sym(2) + sin(x)); # Define the fraction explicitly as symbolic value 3 / 2

julia> solve(f2(x), x)
2-element Vector{Sym}:
   pi/2
 3*pi/2

Note that IntervalRootFinding doesn’t work now anymore, since it doesn’t know how to use Sym values.

Hope that helps!

2 Likes

Just to add to this excellent answer, using a rational (3//2) instead of the float (1.5 or 3/2) should be enough, as rationals promote to symbolic equivalents.

2 Likes

Very convenient, thanks for the hint!

It helps a lot thanks for this amazing guidance.

1 Like