JuliaIntervals/ValidatedNumerics.jl

Only finds one root of tan

Opened this issue · 6 comments

dpsanders commented 
julia> find_roots(tan, -50, 50)
1-element Array{ValidatedNumerics.RootFinding.Root{Float64},1}:
 ValidatedNumerics.RootFinding.Root{Float64}([3.14159, 3.1416],:unique)
dpsanders commented

And yet

julia> find_roots(tan, 30..40)
1-element Array{ValidatedNumerics.RootFinding.Root{Float64},1}:
 ValidatedNumerics.RootFinding.Root{Float64}([34.5575, 34.5576],:unique)

Which is wrongcorrect:

julia> tan(34.5575)
-1.9189487730536465e-5

julia> 34.5575/pi
10.999993891796345
dpsanders commented

This seems to be because the Newton method behaves badly for this function, presumably due to the singularities.

dpsanders commented

This can be captured with decorations:

julia> displaymode(decorations=true)

julia> x = @decorated(30, 40)
[30, 40]_com

julia> tan(x)
[-∞, ∞]_trv

julia> ForwardDiff.derivative(tan, x)
[1, ∞]_trv
lbenet commented

I think that there are two issues here: First, the theorem requires the function needs to be twice continuously differentiable in the interval, which I think is useful if the interval contains no singularities.

Secondly, I think that (part of) the odd behavior is related to guarded_mid, which is somehow defined in relative terms:

julia> RootFinding.guarded_mid(-5pi..5pi)
0.9817477042468106

julia> RootFinding.guarded_mid(-50pi..50pi)
9.817477042468113

In both cases, the shifted point is actually too far from zero.

The consequence is that the Newton operator N then yields something different:

julia> RootFinding.N(tan,-50pi..50pi, 1.0..Inf)
[9.40326, 9.81748]

julia> RootFinding.N(tan,-5pi..5pi, 1.0..Inf)
[-0.514859, 0.981748]
dpsanders commented

Somehow we need to encode in tan where the discontinuities are. ApproxFun does this "automatically", I believe by treating tan as sin/cos and knowing (or finding!) the roots of cos.

dpsanders commented

Note that the following works. It is currently unable to recognise the singularity EDIT as a singularity, but it correctly reports that it exists and is not known to be a root; and the roots are correctly identified.

julia v0.5> newton(x->sin(x)/cos(x), -10..10, tolerance=1e-3)
13-element Array{ValidatedNumerics.RootFinding.Root{Float64},1}:
 ValidatedNumerics.RootFinding.Root{Float64}([-9.42485, -9.42471],:unique)
 ValidatedNumerics.RootFinding.Root{Float64}([-7.85412, -7.85346],:unknown)
 ValidatedNumerics.RootFinding.Root{Float64}([-6.28319, -6.28318],:unique)
 ValidatedNumerics.RootFinding.Root{Float64}([-4.7128, -4.71214],:unknown)
 ValidatedNumerics.RootFinding.Root{Float64}([-3.1416, -3.14159],:unique)
 ValidatedNumerics.RootFinding.Root{Float64}([-1.57082, -1.57016],:unknown)
 ValidatedNumerics.RootFinding.Root{Float64}([-1.45505e-06, 2.08322e-06],:unique)
 ValidatedNumerics.RootFinding.Root{Float64}([1.57028, 1.57086],:unknown)
 ValidatedNumerics.RootFinding.Root{Float64}([3.14159, 3.1416],:unique)
 ValidatedNumerics.RootFinding.Root{Float64}([4.71225, 4.71283],:unknown)
 ValidatedNumerics.RootFinding.Root{Float64}([6.28318, 6.28319],:unique)
 ValidatedNumerics.RootFinding.Root{Float64}([7.85366, 7.85424],:unknown)
 ValidatedNumerics.RootFinding.Root{Float64}([9.42477, 9.42479],:unique)