JuliaIntervals / TaylorModels.jl

Rigorous function approximation using Taylor models in Julia
Other
63 stars 14 forks source link

AssertionError in validated_integ if rhs has a sine term #109

Closed mforets closed 3 years ago

mforets commented 3 years ago
using Test
using TaylorModels
using TaylorModels: validated_integ

function pendulum!(dx, x, p, t)
    aux = 2 * sin(x[1])
    dx[1] = x[2]
    dx[2] = aux + 8*x[3]
    dx[3] = zero(x[1])
    return nothing
end

# Initial conditions
tini, tend = 0.0, 0.03
q0 = [1.2, 0.2, -0.78]
δq0 = IntervalBox(0 .. 0, 3)
X0 = IntervalBox(q0 .+ δq0)

# Parameters
abstol = 1e-20
orderQ = 3
orderT = 10
ξ = set_variables("ξ", order=2*orderQ, numvars=length(q0))

normalized_box = IntervalBox(-1 .. 1, Val(orderQ))

tTM, qv, qTM = validated_integ(pendulum!, X0, tini, tend, orderQ, orderT, abstol)

AssertionError: get_numvars() == length(vals)

Stacktrace:
  [1] _evaluate(a::TaylorN{Interval{Float64}}, vals::Tuple{Interval{Float64}})
    @ TaylorSeries ~/.julia/packages/TaylorSeries/tveWm/src/evaluate.jl:229
  [2] _evaluate
    @ ~/.julia/packages/TaylorSeries/tveWm/src/evaluate.jl:240 [inlined]
  [3] evaluate(a::TaylorN{Interval{Float64}}, vals::Interval{Float64}; sorting::Bool)
    @ TaylorSeries ~/.julia/packages/TaylorSeries/tveWm/src/evaluate.jl:223
  [4] evaluate
    @ ~/.julia/packages/TaylorSeries/tveWm/src/evaluate.jl:223 [inlined]
  [5] _broadcast_getindex_evalf
    @ ./broadcast.jl:648 [inlined]
  [6] _broadcast_getindex
    @ ./broadcast.jl:621 [inlined]
  [7] getindex
    @ ./broadcast.jl:575 [inlined]
  [8] macro expansion
    @ ./broadcast.jl:984 [inlined]
  [9] macro expansion
    @ ./simdloop.jl:77 [inlined]
 [10] copyto!
    @ ./broadcast.jl:983 [inlined]
 [11] copyto!
    @ ./broadcast.jl:936 [inlined]
 [12] copy
    @ ./broadcast.jl:908 [inlined]
 [13] materialize
    @ ./broadcast.jl:883 [inlined]
 [14] broadcasted
    @ ~/.julia/packages/IntervalArithmetic/dhnJB/src/multidim/arithmetic.jl:29 [inlined]
 [15] picard_remainder!(f!::typeof(pendulum!), t::Taylor1{Float64}, x::Vector{Taylor1{TaylorN{Float64}}}, dx::Vector{Taylor1{TaylorN{Float64}}}, xxI::Vector{Taylor1{TaylorN{Interval{Float64}}}}, dxxI::Vector{Taylor1{TaylorN{Interval{Float64}}}}, δI::IntervalBox{3, Float64}, δt::Interval{Float64}, Δx::IntervalBox{3, Float64}, Δ0::IntervalBox{3, Float64}, params::Nothing)
    @ TaylorModels ~/.julia/dev/TaylorModels/src/validatedODEs.jl:108
 [16] remainder_taylorstep!(f!::Function, t::Taylor1{Float64}, x::Vector{Taylor1{TaylorN{Float64}}}, dx::Vector{Taylor1{TaylorN{Float64}}}, xI::Vector{Taylor1{Interval{Float64}}}, dxI::Vector{Taylor1{Interval{Float64}}}, δI::IntervalBox{3, Float64}, δt::Interval{Float64}, params::Nothing)
    @ TaylorModels ~/.julia/dev/TaylorModels/src/validatedODEs.jl:32
 [17] validated_step!(f!::typeof(pendulum!), t::Taylor1{Float64}, x::Vector{Taylor1{TaylorN{Float64}}}, dx::Vector{Taylor1{TaylorN{Float64}}}, xaux::Vector{Taylor1{TaylorN{Float64}}}, tI::Taylor1{Float64}, xI::Vector{Taylor1{Interval{Float64}}}, dxI::Vector{Taylor1{Interval{Float64}}}, xauxI::Vector{Taylor1{Interval{Float64}}}, t0::Float64, tmax::Float64, sign_tstep::Int64, xTMN::Vector{TaylorModelN{3, Float64, Float64}}, xv::Vector{IntervalBox{3, Float64}}, rem::Vector{Interval{Float64}}, zbox::IntervalBox{3, Float64}, symIbox::IntervalBox{3, Float64}, nsteps::Int64, orderT::Int64, abstol::Float64, params::Nothing, parse_eqs::Bool, check_property::TaylorModels.var"#36#38")
    @ TaylorModels ~/.julia/dev/TaylorModels/src/validatedODEs.jl:361
 [18] validated_integ(f!::Function, X0::IntervalBox{3, Float64}, t0::Float64, tmax::Float64, orderQ::Int64, orderT::Int64, abstol::Float64, params::Nothing; maxsteps::Int64, parse_eqs::Bool, check_property::Function)
    @ TaylorModels ~/.julia/dev/TaylorModels/src/validatedODEs.jl:533
 [19] validated_integ (repeats 2 times)
    @ ~/.julia/dev/TaylorModels/src/validatedODEs.jl:481 [inlined]
 [20] top-level scope
    @ In[15]:27
 [21] eval
    @ ./boot.jl:360 [inlined]
 [22] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
    @ Base ./loading.jl:1094

Let me remark that for aux = 2 there is no error.

lbenet commented 3 years ago

Thanks for reporting! As far as I have checked, this is related with some changes in TaylorSeries.jl, related to evaluate rather than sin. Let me dig into it...