Interpolation problems for SimpleATsit5 ?

[Datseris]

This is the MWE I could create. Seriously!

using SimpleDiffEq, OrdinaryDiffEq, Roots, StaticArrays
const ROOTS_ALG = Roots.A42()
const rootkw = (xrtol = 1e-6, atol = 1e-6)

function loop(u, p, t)
    σ = p[1]; ρ = p[2]; β = p[3]
    du1 = σ*(u[2]-u[1])
    du2 = u[1]*(ρ-u[3]) - u[2]
    du3 = u[1]*u[2] - β*u[3]
    return SVector{3}(du1, du2, du3)
end
prob = ODEProblem{false}(loop, SVector{3}([0, 10.0, 0]), (0.0, Inf),  [10, 28, 8/3])

# Minimal poincare surface of section function (can't get smaller)
function poincare(prob, alg, tfinal = 2000.0)
    integ = init(prob, alg, atol = 1e-6, rtol = 1e-6, save_everystep = false)
    plane(u) = u[2]  # crosses section when 2nd variable crosses 0
    planet(t) = plane(integ(t))

    psos = Vector{SVector{3,Float64}}()
    side = plane(integ.u)

    while integ.t < tfinal
        while side < 0
            integ.t > tfinal && break
            step!(integ)
            side = plane(integ.u)
        end
        while side ≥ 0
            integ.t > tfinal && break
            step!(integ)
            side = plane(integ.u)
        end
        integ.t > tfinal && break

        # After the above two loops I am now guaranteed to have `t` in negative
        # side and `tprev` in positive side. Thus some time inbetween is in
        # exactly side = 0, which means on the PSOS (for the lorentz system)
        tcross = find_zero(planet, (integ.tprev, integ.t), ROOTS_ALG; rootkw...)
        ucross = integ(tcross)
        push!(psos, ucross)
    end
    return psos
end

psos = poincare(prob, Tsit5())
psos2 = poincare(prob, SimpleATsit5())

what you will notice is that psos returns fine (and what is expected i.e. second variable close to zero). However, psos2 errors:

ERROR: ArgumentError: The interval [a,b] is not a bracketing interval.
You need f(a) and f(b) to have different signs (f(a) * f(b) < 0).
Consider a different bracket or try fzero(f, c) with an initial guess c.

I believe this has to be some problem in the interpolation of SimpleATsit5 since I use identical code for both algorithms.

Any ideas?

[Datseris]

and of course this is not a problem of Roots. Adding the line

        @assert planet(integ.t) * planet(integ.tprev) < 0

before the find_zero call spits out the assertion error. I guess this means that interpolation is not working perfectly exactly at the boundaries?

[Datseris]

@ChrisRackauckas oh yeah! I think I really found the problem!!! Check this out: For the SimpleTsit:

julia> integ(integ.t)
3-element SArray{Tuple{3},Float64,1,3}:
 11.071218851365108
 22.9963922538833
  9.966343862342143

julia> integ.u
3-element SArray{Tuple{3},Float64,1,3}:
 10.987883478548856
 22.95174900928437
  9.51383768714711

while for Vern9:

julia> integ(integ.t)
3-element SArray{Tuple{3},Float64,1,3}:
 0.00763628449871269
 9.99924424819131
 2.9169291953420574e-6

julia> integ.u
3-element SArray{Tuple{3},Float64,1,3}:
 0.007636284498717652
 9.99924424819131
 2.9169291953472726e-6

[Datseris]

This was closed although not tested. For me, even though the current test suite of the master branch passes the test, the above code snippet still fails.

Intuitively what the code snippet assumes should be true: the boolean functions that pass using integ.u should also pass using integ(integ.t) but they don't.

You can add the statement @assert planett(integ.t) < 0 after the while loops, which is what the code assumes to be true. (It isn't true otherwise Roots wouldn't fail)