Fix evaluation in picard_remainder! (#110)

* Fix evaluation in picard_remainder!

* Two minor details

... in validated integration methods

* Add test: integrate pendulum with constant torque

* Bump patch versoin

Project.toml

@@ -1,7 +1,7 @@
 name = "TaylorModels"
 uuid = "314ce334-5f6e-57ae-acf6-00b6e903104a"
 repo = "https://github.com/JuliaIntervals/TaylorModels.jl.git"
-version = "0.3.9"
+version = "0.3.10"
 
 [deps]
 IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"

src/validatedODEs.jl

@@ -105,7 +105,7 @@ function picard_remainder!(f!::Function, t::Taylor1{T},
     f!(dxxI, xxI, params, t)
 
     # Picard iteration, considering only the bound of `f` and the last coeff of f
-    Δdx = IntervalBox( evaluate.( (dxxI - dx)(δt), δI ) )
+    Δdx = IntervalBox( evaluate.( (dxxI - dx)(δt), (δI,) ) )
     Δ = Δ0 + Δdx * δt
     return Δ
 end
@@ -540,14 +540,7 @@ function validated_integ(f!, X0, t0::T, tmax::T, orderQ::Int, orderT::Int, absto
     @inbounds xv[1] = evaluate(xTMN, S)
 
     # Determine if specialized jetcoeffs! method exists (built by @taylorize)
-    parse_eqs = parse_eqs && (length(methods(TaylorIntegration.jetcoeffs!)) > 2)
-    if parse_eqs
-        try
-            TaylorIntegration.jetcoeffs!(Val(f!), t, x, dx, params)
-        catch
-            parse_eqs = false
-        end
-    end
+    parse_eqs = TaylorIntegration._determine_parsing!(parse_eqs, f!, t, x, dx, params)
 
     # Integration
     nsteps = 1
@@ -752,10 +745,9 @@ function validated_integ2(f!, X0, t0::T, tf::T, orderQ::Int, orderT::Int,
         δt = sign_tstep * δt
 
         @inbounds for i in eachindex(x)
-            dom = sign_tstep > 0 ? Interval(0, δt) : Interval(δt, 0)
-            x0 = sign_tstep > 0 ? dom.lo : dom.hi
-            Δ = zero(Interval{Float64})
-            xTM1[i] = TaylorModel1(deepcopy(x[i]), Δ, x0, dom)
+            dom = sign_tstep * Interval{T}(0, sign_tstep*δt)
+            Δ = zero(dom)
+            xTM1[i] = TaylorModel1(deepcopy(x[i]), Δ, zI, dom)
         end
 
         # to reuse the previous TaylorModel and save some allocations

test/validated_integ.jl

@@ -235,4 +235,52 @@ interval_rand(X::IntervalBox) = interval_rand.(X)
             end
         end
     end
+
+    @testset "Pendulum with constant torque" begin
+        @taylorize function pendulum!(dx, x, p, t)
+            si = sin(x[1])
+            aux = 2 *  si
+            dx[1] = x[2]
+            dx[2] = aux + 8*x[3]
+            dx[3] = zero(x[1])
+            nothing
+        end
+        # Conseerved quantity
+        ene_pendulum(x) = x[2]^2/2 + 2 * cos(x[1]) - 8 * x[3]
+
+        # Initial conditions
+        tini, tend = 0.0, 12.0
+        q0 = [1.1, 0.1, 0.0]
+        δq0 = IntervalBox(-0.1 .. 0.1, -0.1 .. 0.1, 0..0)
+        X0 = IntervalBox(q0 .+ δq0)
+        ene0 = ene_pendulum(X0)
+
+        # Parameters
+        abstol = 1e-20
+        orderQ = 3
+        orderT = 10
+        ξ = set_variables("ξ", order=2*orderQ, numvars=length(q0))
+
+        tTM, qv, qTM = validated_integ(pendulum!, X0, tini, tend, orderQ, orderT, abstol,
+            maxsteps=1800);
+        @test all(ene0 .⊆ ene_pendulum.(qv))
+
+        # tTM, qv, qTM = validated_integ2(pendulum!, X0, tini, tend, orderQ, orderT, abstol,
+        #     maxsteps=1800);
+        # @test all(ene0 .⊆ ene_pendulum.(qv))
+
+        # Initial conditions 2
+        q0 = [1.1, 0.1, 0.0]
+        δq0 = IntervalBox(-0.1 .. 0.1, -0.1 .. 0.1, -0.01 .. 0.01)
+        X0 = IntervalBox(q0 .+ δq0)
+        ene0 = ene_pendulum(X0)
+
+        tTM, qv, qTM = validated_integ(pendulum!, X0, tini, tend, orderQ, orderT, abstol,
+            maxsteps=1800);
+        @test all(ene0 .⊆ ene_pendulum.(qv))
+
+        # tTM, qv, qTM = validated_integ2(pendulum!, X0, tini, tend, orderQ, orderT, abstol,
+        #     maxsteps=1800);
+        # @test all(ene0 .⊆ ene_pendulum.(qv))
+    end
 end