BoundsError for nonlinear ODE

[SebastianGuadalupe]

using ReachabilityAnalysis

@taylorize function benchmark7!(du, u, p, t)
    
    du[1] = u[3]^3 - u[2] + u[4]
    du[2] = u[3]
    du[3] = 2.0
    du[4] = u[4]
    return du
end

# define the initial-value problem
X₀ = Hyperrectangle(low=[0.35, 0.45, 0.25, -0.01], high=[0.45, 0.55, 0.35, 0.01])

prob = @ivp(x' = benchmark7!(x), dim: 4, x(0) ∈ X₀)

# solve it
sol = solve(prob, T=2.0)

Stacktrace:
 [1] getindex at ./array.jl:787 [inlined]
 [2] getindex at /home/sguadalupe/.julia/packages/TaylorSeries/XWMlP/src/auxiliary.jl:80 [inlined]
 [3] jetcoeffs!(::typeof(benchmark7!), ::Taylor1{Float64}, ::Array{Taylor1{TaylorN{Float64}},1}, ::Array{Taylor1{TaylorN{Float64}},1}, ::Array{Taylor1{TaylorN{Float64}},1}, ::Nothing) at /home/sguadalupe/.julia/packages/TaylorIntegration/KHxcX/src/explicitode.jl:82
 [4] __jetcoeffs!(::Val{false}, ::Function, ::Taylor1{Float64}, ::Array{Taylor1{TaylorN{Float64}},1}, ::Array{Taylor1{TaylorN{Float64}},1}, ::Array{Taylor1{TaylorN{Float64}},1}, ::Nothing) at /home/sguadalupe/.julia/packages/TaylorIntegration/KHxcX/src/explicitode.jl:101
 [5] taylorstep!(::Function, ::Taylor1{Float64}, ::Array{Taylor1{TaylorN{Float64}},1}, ::Array{Taylor1{TaylorN{Float64}},1}, ::Array{Taylor1{TaylorN{Float64}},1}, ::Float64, ::Nothing, ::Bool) at /home/sguadalupe/.julia/packages/TaylorIntegration/KHxcX/src/explicitode.jl:237
 [6] validated_step!(::typeof(benchmark7!), ::Taylor1{Float64}, ::Array{Taylor1{TaylorN{Float64}},1}, ::Array{Taylor1{TaylorN{Float64}},1}, ::Array{Taylor1{TaylorN{Float64}},1}, ::Taylor1{Float64}, ::Array{Taylor1{IntervalArithmetic.Interval{Float64}},1}, ::Array{Taylor1{IntervalArithmetic.Interval{Float64}},1}, ::Array{Taylor1{IntervalArithmetic.Interval{Float64}},1}, ::Float64, ::Float64, ::Int64, ::Array{TaylorModels.TaylorModelN{4,Float64,Float64},1}, ::Array{IntervalBox{4,Float64},1}, ::Array{IntervalArithmetic.Interval{Float64},1}, ::IntervalBox{4,Float64}, ::IntervalBox{4,Float64}, ::Int64, ::Int64, ::Float64, ::Nothing, ::Bool, ::ReachabilityAnalysis.var"#111#113") at /home/sguadalupe/.julia/dev/ReachabilityAnalysis/src/Algorithms/TMJets/validated_integ.jl:365
 [7] validated_integ!(::Array{TaylorModelReachSet{Float64},1}, ::Function, ::Array{Float64,1}, ::IntervalBox{4,Float64}, ::Float64, ::Float64, ::Int64, ::Int64, ::Float64, ::Int64, ::Nothing; parse_eqs::Bool, check_property::Function) at /home/sguadalupe/.julia/dev/ReachabilityAnalysis/src/Algorithms/TMJets/validated_integ.jl:492
 [8] validated_integ! at /home/sguadalupe/.julia/dev/ReachabilityAnalysis/src/Algorithms/TMJets/validated_integ.jl:431 [inlined] (repeats 2 times)
 [9] post(::TMJets{Float64}, ::InitialValueProblem{BlackBoxContinuousSystem{typeof(benchmark7!)},Hyperrectangle{Float64,Array{Float64,1},Array{Float64,1}}}, ::IntervalArithmetic.Interval{Float64}; kwargs::Base.Iterators.Pairs{Symbol,Float64,Tuple{Symbol},NamedTuple{(:T,),Tuple{Float64}}}) at /home/sguadalupe/.julia/dev/ReachabilityAnalysis/src/Algorithms/TMJets/post.jl:27
 [10] solve(::InitialValueProblem{BlackBoxContinuousSystem{typeof(benchmark7!)},Hyperrectangle{Float64,Array{Float64,1},Array{Float64,1}}}; kwargs::Base.Iterators.Pairs{Symbol,Float64,Tuple{Symbol},NamedTuple{(:T,),Tuple{Float64}}}) at /home/sguadalupe/.julia/dev/ReachabilityAnalysis/src/Continuous/solve.jl:34
 [11] top-level scope at In[21]:16

[mforets]

The workaround that i know is to define du[3] = 2.0 * one(u[3]). I have marked this issue as external as it should probably be considered in TaylorIntegration. cc'ing @lbenet who may have a better answer.

@SebastianGuadalupe can you try to see if this problem also arises if you use the DifferentialEquations interface directly? (cf. https://perezhz.github.io/TaylorIntegration.jl/stable/common/)

[lbenet]

Thanks for reporting. I agree with the workaround of @mforets.

I would also suggest using the following:

@taylorize function benchmark7!(du, u, p, t)
    local two = 2.0+zero(u[1])
    du[1] = u[3]^3 - u[2] + u[4]
    du[2] = u[3]
    du[3] = two
    du[4] = u[4]
    return du
end

[lbenet]

can you try to see if this problem also arises if you use the DifferentialEquations interface directly?

I also think the prolem will still be present

[mforets]

Alright, thanks. Should we report this in TaylorIntegration.jl or it is known?

I just uploaded some documentation for this library related to this (and some other common issues), see #180. Comments welcome :)

[lbenet]

Yes, please, report to TaylorIntegration.jl, and if possible, add further comments in the the docs as well.

[mforets]

👍 @SebastianGuadalupe i'm looking at you 😄

[mforets]

xref: PerezHz/TaylorIntegration.jl#104