rework to save bracket

[jverzani]

Introduces find_bracket to return a tuple containing a zero, an enclosing bracket, and a flat indicating if the value is an exact zero.

To do so, required adding a field in the struct to store the approximate zero instead of reusing xn1, which contains one part of the bracket.

[ChrisRackauckas]

@kanav99

[jverzani]

@kanav99

Should I merge this? I certainly can, but wanted to make sure the interface is sufficient before doing so.

[kanav99]

Sorry @jverzani this one slipped my mind. We have developed a workaround for the time being which seems to be stable. What we did right now is instead of solving for θ (which is in [0,1]) we solve for 1-θ. Thus we get an upper limit for 1-θ i.e. lower limit of θ. But this can only be done using Bisection.

If this PR makes sure that xn0 is just less than the root and xn1 is just greater than the root then this works for us!

[jverzani]

There are two possibilities due to the stopping rules employed: one with a bracket satisfying the condition or one where the algorithm terminates with an exact zero is found. In the latter case, the bracket may not satisfy the above. The return value includes a flag to determine which occurred. Is that a satisfactory design for you?

[kanav99]

I think that works for us. If we do prevfloat(root) in case of exact solution, we will get the previous sign right?

[jverzani]

No, that's the worry. For example, you might see this:

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

julia> Roots.find_bracket(f, (-2, 2))
(xstar = 0.0, bracket = (0.0, 2.0), exact = true)

Here, the first step is to consider the midpoint 0.0 which is an exact zero. But at zero, there is no sign change.

[kanav99]

I don't think that would be issue, isn't it @ChrisRackauckas ? This issue is not because of floating point error or inconsistent lower or upper limit but a property of the function. If Chris agrees too, this interface is good enough for us.

[ChrisRackauckas]

The true problem is that f(x) = 1e-12 and then f(prevfloat(x)) = 1e-12 still because of precision stuff, so f(prevfloat(prevfloat(x)) = -1e-12 is what we really needed. So we need the bracket and can just take the left, and that would solve it. The issue is that we've been having to assume the bracket right prevfloat once is equal to the left, which isn't always true. If we have a left bracket that always works, 👍

[jverzani]

Hi @ChrisRackauckas and @kanav99 . After your comments, I squeezed the tolerances for the Alefeld-Potra-Shi cases and the following is satisfied for the test cases: if l is the return value from find_bracket then either:

• l.exact is true (in which case the bracket may be quite large, as the algorithm can terminate early)

• abs(l.bracket[2] - l.bracket[1]) <= 2eps(maximum(abs.(l.bracket))). For BisectionExact, this bound is smaller (the bracket consists of adjacent floting point values), but this isn't the case for the faster A42 method.

In the latter case, l.xstar can be any of the 4 possible values in the bracket, so starting at l.xstar, you might need to step back 3 floats to get the sign change.

Hopefully that will suit for your use case.

[ChrisRackauckas]

In the latter case, l.xstar can be any of the 4 possible values in the bracket, so starting at l.xstar, you might need to step back 3 floats to get the sign change.

Yeah, having a known fact that it's at most 4 steps fixes this! Before we just put to run it like 5 steps otherwise through an error, so the larger tolerance likely was why ewe hit the error.

[kanav99]

@jverzani In case of exact = false what is the suggested root? I @showed a result of the find_bracket function and output was:

(xstar = 0.9746722903038366, bracket = (0.9746722903038367, 0.9746722903038368), exact = false)

Note that xstar < bracket[1] < bracket[2]. Is this expected behaviour?

[jverzani]

Well, that isn't expected. Do you have the function available that I can look at?

[kanav99]

It is a private code someone sent me. I don't think I can send that to you. I will try to make a MWE

[jverzani]

Thanks. If you can, that would be great. If it proves too much work, let me try to find the issue without an example.

…

On Fri, Mar 13, 2020 at 12:14 PM Kanav Gupta ***@***.***> wrote: It is a private code someone sent me. I don't think I can send that to you. I will try to make a MWE — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#173 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AADG6TGUKPHPON57UB3JHTLRHJLURANCNFSM4KR4RA2Q> .

-- John Verzani Department of Mathematics College of Staten Island, CUNY verzani@math.csi.cuny.edu

[jverzani]

I think I found the culprit here. I made a PR #174. If you can check, it would be most appreciated.

…

On Fri, Mar 13, 2020 at 12:21 PM John Verzani ***@***.***> wrote: Thanks. If you can, that would be great. If it proves too much work, let me try to find the issue without an example. On Fri, Mar 13, 2020 at 12:14 PM Kanav Gupta ***@***.***> wrote: > It is a private code someone sent me. I don't think I can send that to > you. I will try to make a MWE > > — > You are receiving this because you were mentioned. > Reply to this email directly, view it on GitHub > <#173 (comment)>, > or unsubscribe > <https://github.com/notifications/unsubscribe-auth/AADG6TGUKPHPON57UB3JHTLRHJLURANCNFSM4KR4RA2Q> > . > -- John Verzani Department of Mathematics College of Staten Island, CUNY ***@***.***

-- John Verzani Department of Mathematics College of Staten Island, CUNY verzani@math.csi.cuny.edu