Convergence logic sometimes fails to notice that no progress is being made

[ikirill]

Convergence logic sometimes fails to notice that no progress is being made. I don't have any idea why this particular test function triggers that behaviour.

module M1

import Roots
import ForwardDiff

function check_rootfinding(which; args...)
    g, T = 1.62850, 14.60000
    α, t1, tf = 0.00347, 40.91375, 131.86573
    y, ya, yf = 0.0, 9000.0, 10000.0
    vy = sqrt(2g*(ya-y))
    θ0, θ1 = atan(α*tf), atan(α*(tf-t1))
    nf = 0
    I_sintan(x) = tan(x)/2cos(x) - atanh(tan(x/2))
    I_sintan(x, y) = I_sintan(y) - I_sintan(x)
    function lhs(θ)
        nf += 1
        tRem = (vy - T/α*(sec(θ1) - sec(θ))) / g
        val = -yf + y + vy*tRem - 0.5g*tRem^2 - T/α^2*I_sintan(θ, θ1)
        isa(θ, AbstractFloat) && @printf "[% 3d]: θ=% .18e val=% .18e\n" nf θ val
        val
    end
    ans = Roots.find_zero(lhs, [atan(α*tf), atan(α*(tf-t1))], which; args...)
    @show nf
    ans_d = ForwardDiff.derivative(x -> lhs(x), ans)
    @show ans
    @show ans_d
    abs(lhs(ans) / ans_d)
end

end

Notice that it keeps calling the function at exactly the same two points starting from iteration 10, making no progress:

julia> M1.check_rootfinding(Roots.FalsePosition())
[  1]: θ= 3.057096351246670896e-01 val=-1.000000000000000000e+03
[  2]: θ= 3.057096351246670896e-01 val=-1.000000000000000000e+03
[  3]: θ= 4.291346551666863629e-01 val= 8.995313325134375191e+03
[  4]: θ= 3.180579243705532466e-01 val= 5.404732458036085063e+02
[  5]: θ= 3.137255414570386813e-01 val= 5.773644115446813885e+00
[  6]: θ= 3.136790321558284855e-01 val=-2.980464060328813503e-03
[  7]: θ= 3.136790561524200882e-01 val= 1.918206180562265217e-07
[  8]: θ= 3.136790561508757680e-01 val= 1.068656274583190680e-11
[  9]: θ= 3.136790561508756570e-01 val=-2.296474121976643801e-11
[ 10]: θ= 3.136790561508757125e-01 val=-2.296474121976643801e-11
[ 11]: θ= 3.136790561508757680e-01 val= 1.068656274583190680e-11
[ 12]: θ= 3.136790561508757680e-01 val= 1.068656274583190680e-11
[ 13]: θ= 3.136790561508757680e-01 val= 1.068656274583190680e-11
[ 14]: θ= 3.136790561508757125e-01 val=-2.296474121976643801e-11
[ 15]: θ= 3.136790561508757680e-01 val= 1.068656274583190680e-11
[ 16]: θ= 3.136790561508757680e-01 val= 1.068656274583190680e-11
[ 17]: θ= 3.136790561508757680e-01 val= 1.068656274583190680e-11
[ 18]: θ= 3.136790561508757125e-01 val=-2.296474121976643801e-11
[ 19]: θ= 3.136790561508757680e-01 val= 1.068656274583190680e-11
[ 20]: θ= 3.136790561508757680e-01 val= 1.068656274583190680e-11
[ 21]: θ= 3.136790561508757680e-01 val= 1.068656274583190680e-11
[ 22]: θ= 3.136790561508757125e-01 val=-2.296474121976643801e-11
[ 23]: θ= 3.136790561508757680e-01 val= 1.068656274583190680e-11
[ 24]: θ= 3.136790561508757680e-01 val= 1.068656274583190680e-11
[ 25]: θ= 3.136790561508757680e-01 val= 1.068656274583190680e-11
[ 26]: θ= 3.136790561508757125e-01 val=-2.296474121976643801e-11
[ 27]: θ= 3.136790561508757680e-01 val= 1.068656274583190680e-11
[ 28]: θ= 3.136790561508757680e-01 val= 1.068656274583190680e-11
[ 29]: θ= 3.136790561508757680e-01 val= 1.068656274583190680e-11
[ 30]: θ= 3.136790561508757125e-01 val=-2.296474121976643801e-11
[ 31]: θ= 3.136790561508757680e-01 val= 1.068656274583190680e-11
[ 32]: θ= 3.136790561508757680e-01 val= 1.068656274583190680e-11
[ 33]: θ= 3.136790561508757680e-01 val= 1.068656274583190680e-11
[ 34]: θ= 3.136790561508757125e-01 val=-2.296474121976643801e-11
[ 35]: θ= 3.136790561508757680e-01 val= 1.068656274583190680e-11
[ 36]: θ= 3.136790561508757680e-01 val= 1.068656274583190680e-11
[ 37]: θ= 3.136790561508757680e-01 val= 1.068656274583190680e-11
[ 38]: θ= 3.136790561508757125e-01 val=-2.296474121976643801e-11
[ 39]: θ= 3.136790561508757680e-01 val= 1.068656274583190680e-11
[ 40]: θ= 3.136790561508757680e-01 val= 1.068656274583190680e-11
[ 41]: θ= 3.136790561508757680e-01 val= 1.068656274583190680e-11
[ 42]: θ= 3.136790561508757125e-01 val=-2.296474121976643801e-11
[ 43]: θ= 3.136790561508757680e-01 val= 1.068656274583190680e-11
[ 44]: θ= 3.136790561508757680e-01 val= 1.068656274583190680e-11
nf = 44
ans = 0.31367905615087577
ans_d = 124211.6309131711
[ 46]: θ= 3.136790561508757680e-01 val= 1.068656274583190680e-11
8.603512140744889e-17

This is with Pkg.checkout("Roots"):

julia> Pkg.status("Roots")
 - Roots                         0.4.1+             master

Here's a plot of that function, it doesn't look pathological: https://imgur.com/a/KWmTl

[jverzani]

This case is an interesting problem. The algorithm is set to stop when the f(x) values get close to 0, but not when subsequent x values get close (as xabstol=zero and not even eps()). It seems to suggest relaxing the bounds a tad is a good idea. (check_rootfinding(Roots.FalsePosition(), xabstol=eps()) takes 8 steps.

I'll make this change.