JuliaIntervals / IntervalRootFinding.jl

Library for finding the roots of a function using interval arithmetic
https://juliaintervals.github.io/IntervalRootFinding.jl/
Other
125 stars 26 forks source link

Add root-status "root(s)" ? #164

Closed KronosTheLate closed 3 years ago

KronosTheLate commented 3 years ago

In finding roots, I feel there is a huge difference between no roots and [1, many] roots. These two cases are however currently represented by the status :unknown.

In three cases I have seen (out of not that many unknown roots), a simple check was available to rule out the possibility of no roots (where the interval was given to the function, and the returned interval contained zero. A code example is provided at the bottom.

Would it not be smart to compute such a test, and change the status from :unknown to no root (leading to no return value, requiring no status), or e.g. root(s)?

Below is a code example (excuse the units) of returned roots of unknown status, that seem like they could easily be "upgraded" to status root(s) by the check in the end. A plot is included for convenience. Kudos for smooth integration with Unitful BTW ^_^

using Unitful, IntervalArithmetic, IntervalRootFinding
const U = Unitful
k_B     = U.k
mp      = U.mp
h_planc = U.h

T = 1e3u"K"
R = 13.6u"eV"   # Ground state energy, hydrogen atom
β = 1/(k_B*T)
n = 1e20u"m^-3"
n_Q = /((2*pi*mp*k_B*T)^(3/2), h_planc^3) |> u"m^-3"

# Equation to solve: y^2/(1-y) = A

A = exp(-β*R)*n_Q/n

g(y) = A - frac(y^2, 1-y)

root = roots(g, -1..3)
root[1].interval |> g
root[2].interval |> g

image

image

lucaferranti commented 3 years ago

if the interval returned by the function contains zero, unfortunately you cannot be sure that the function will contain a zero, because the 0 being included might be due to overestimation, like in the following example:

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

julia> f(3.1..3.5)
[-0.890001, 2.95001]

0 is in f(3.1..3.5) but you can easily see that the function does not have a zero in there.

I cannot reproduce your code because of a "tau undefined error" when I try to define n_Q, but why do you think the first interval should be upgrated to root(s)? I think in this case the range [-inf, inf] is because of the verical asymptote, but the interval does not actually contain a root?

lucaferranti commented 3 years ago

For the second interval, I think it's unknown because the root seems to be close to a local minimum. This means that the derivative of the function evaluated at that interval will contain a zero and thus Newton method cannot deduce the unicity. There is an explanation article about interval newton method in the organization webpage here

KronosTheLate commented 3 years ago

I cannot reproduce your code because of a "tau undefined error" when I try to define n_Q, but why do you think the first interval should be upgrated to root(s)?

Sorry about that. I like using Tau instead of 2*pi, I have changed the example now.

I think in this case the range [-inf, inf] is because of the verical asymptote, but the interval does not actually contain a root?

I have not concidered that posibility to be honest. When I plot it (where Plots.jl arbitrarily chooses some points and draws a connecting line), the line crosses the x-axis, and in my simple understanding of this domain I swiftly conclude that a zero is present. I do not know the math well enough to see why there is not nessecarily a root somewhere in the asymptote, but I accept that there is not nessecarily one there.

I had not though of overestimation. I was inspired by an example in which one could use the logic I did in my example to rule out a root, but one can clearly not use the same logic to prove a root. My bad.

You example is very much illuminating. I think you have also completely killed the premise of this issue, and so, I concider it closed. Thanks for the intelligible reply ^_^

(Aand you wrote a second comment:)

For the second interval, I think it's unknown because the root seems to be close to a local minimum (or even be the local minimum itself). This means that the derivative of the function evaluated at that interval will contain a zero and thus Newton method cannot deduce the unicity. There is an explanation article about interval newton method in the organization webpage here

That makes a lot of sense! I should check out the possible contractors in the docs, and try another one then I suppose ^_^

KronosTheLate commented 3 years ago

For the second interval, I think it's unknown because the root seems to be close to a local minimum. This means that the derivative of the function evaluated at that interval will contain a zero and thus Newton method cannot deduce the unicity. There is an explanation article about interval newton method in the organization webpage here

That is a great and visual resource, thanks!

I tried the two other methods, but I ended up getting the exact same roots. It feels like that should not be the case, but below you can see my computations: image