Merge fde2df2 into 32297dd

src/TaylorModels.jl

@@ -35,7 +35,7 @@ export TaylorModel1, RTaylorModel1, TaylorModelN
 
 export remainder, polynomial, domain, expansion_point,
     rpa, fp_rpa, bound_remainder,
-    validated_integ
+    validated_integ, validated_integ2
 
 export linear_dominated_bounder, quadratic_fast_bounder
 

src/integration.jl

@@ -17,7 +17,15 @@ function integrate(a::TaylorModel1{T,S}, c0::T) where {T,S}
 
     return TaylorModel1( integ_pol, Δ, a.x0, a.dom )
 end
+function integrate(a::TaylorModel1{T, S}, c0, cc0) where {T <: TaylorN, S}
+    integ_pol = integrate(a.pol, c0)
+    δ = a.dom - a.x0
+    Δ = bound_integration(a, δ, cc0)
+    return TaylorModel1(integ_pol, Δ, a.x0, a.dom)
+end
+
 integrate(a::TaylorModel1{T,S}) where {T,S} = integrate(a, zero(T))
+integrate(a::TaylorModel1{T, S}, δI) where {T <: TaylorN, S} = integrate(a, zero(T), δI)
 
 function integrate(a::RTaylorModel1{T,S}, c0::T) where {T,S}
     order = get_order(a)
@@ -79,6 +87,13 @@ function bound_integration(a::TaylorModel1{T,S}, δ) where {T,S}
     Δ = δ * (remainder(a) + getcoeff(polynomial(a), order) * aux)
     return Δ
 end
+function bound_integration(a::TaylorModel1{T, S}, δ, δI) where {T <: TaylorN, S}
+    order = get_order(a)
+    aux = δ^order / (order+1)
+    Δ = δ * (remainder(a) + getcoeff(polynomial(a), order)(δI) * aux)
+    return Δ
+end
+
 function bound_integration(a::Vector{TaylorModel1{T,S}}, δ) where {T,S}
     order = get_order(a[1])
     aux = δ^order / (order+1)

src/validatedODEs.jl

@@ -497,3 +497,206 @@ function validated_integ(f!, qq0::AbstractArray{T,1}, δq0::IntervalBox{N,T},
 
     return view(tv,1:nsteps), view(xv,1:nsteps), view(xTM1v, :, 1:nsteps)
 end
+
+"""
+    picard(dx, x0, box)
+
+Computes the picard (integral) operator for the initial condition `x0`.
+`dx` must be the rhs of the differential equation.
+"""
+function picard(dx, x0, box)
+    ∫f = integrate(dx, 0., box)
+    pol = ∫f.pol + x0.pol # picard operator
+    return TaylorModel1(deepcopy(pol), ∫f.rem + x0.rem, ∫f.x0, ∫f.dom)
+end
+
+function _picard(dx, x0, box)
+    ∫f = integrate(dx, 0., box)
+    pol = ∫f.pol + x0.pol
+    Δk = ∫f.rem
+    return pol, Δk
+end
+
+"""
+    verify_contraction(f!, dx, xTM1K, params, t, x0, box)
+
+Verifies contraction property for a given set of `f!`, `dx`, `xTM1K`,
+`params`, `t`, `x0`, `box`. It asserts the contention remainder(∫f!(dx, xTM1K, params, t)) ⊆ remainder(xTM1K).
+"""
+function verify_contraction(f!, dx, xTM1K, params, t, x0, box)
+    f!(dx, xTM1K, params, t)
+    @inbounds for i in eachindex(dx)
+        pii, Δ = _picard(dx[i], x0[i], box)
+        @assert (Δ + x0[i].rem) ⊆ xTM1K[i].rem
+    end
+end
+
+function _bound_integration(a::Taylor1, Δr, δ, δI)
+    order = get_order(a)
+    aux = δ^order / (order+1)
+    Δ = δ * (Δr + getcoeff(a, order)(δI) * aux)
+    return Δ
+end
+
+function _validate_step!(xTM1K, f!, dx, x0, params, t, box, dof; ε=1e-10, δ=1e-5,
+                      maxsteps=20, extrasteps=50)
+    E = [interval(0) for i in eachindex(x0)]
+    E′ = [interval(0) for _ in 1:dof]
+    polv = polynomial.(xTM1K)
+    x0v = expansion_point.(xTM1K)
+    domv = domain.(xTM1K)
+    low_ratiov = Array{Float64, 1}(undef, dof)
+    hi_ratiov = Array{Float64, 1}(undef, dof)
+    nsteps = 0
+
+    while nsteps < maxsteps
+        nsteps += 1
+        f!(dx, xTM1K, params, t)
+        @inbounds for i in eachindex(dx)
+            pᵢ₁, Δ = _picard(dx[i], x0[i], box)
+            E′[i] = Δ + x0[i].rem
+        end
+
+        # Only inflates the required component
+        @inbounds for i in eachindex(dx)
+            if E′[i] ⊆ E[i] || E′[i] == E[i] == Interval(0)
+                xTM1K[i] = TaylorModel1(polv[i], E′[i], x0v[i], domv[i])
+            else
+                εi = (1 - ε) .. (1 + ε)
+                δi = -δ .. δ
+                E[i] = E′[i] * εi + δi
+                xTM1K[i] = TaylorModel1(polv[i], E[i], x0v[i], domv[i])
+            end
+        end
+
+        all(E′ .⊆ E) && break
+
+    end
+
+    @inbounds for i in eachindex(dx)
+        low_ratiov[i] = E′[i].lo / E[i].lo
+        hi_ratiov[i] = E′[i].hi / E[i].hi
+    end
+
+    # Contract further the remainders if the last contraction improves more than 5%
+    for ind = 1:extrasteps
+        all(low_ratiov .> 0.95) && all(hi_ratiov .> 0.95) && break
+        f!(dx, xTM1K, params, t)
+        @inbounds for i in eachindex(dx)
+            E[i] = E′[i]
+            _, Δ = _picard(dx[i], x0[i], box)
+            E′[i] = Δ + x0[i].rem
+            xTM1K[i] = TaylorModel1(polv[i], E′[i], x0v[i], domv[i])
+            low_ratiov[i] = E′[i].lo / E[i].lo
+            hi_ratiov[i] = E′[i].hi / E[i].hi
+        end
+    end
+
+
+    verify_contraction(f!, dx, xTM1K, params, t, x0, box)
+
+    if nsteps == maxsteps
+        @warn ("Maximum number of validate steps reached.")
+    end
+
+    return nothing
+end
+
+function validated_integ2(f!, qq0, δq0::IntervalBox{N, T}, t0, tf, orderQ, orderT,
+                         abstol, params=nothing; parse_eqs=true, maxsteps=500,
+                         validatesteps=30, ε=1e-10, δ=1e-3,
+                         absorb_steps=3) where {N, T}
+    dof = N
+    @assert N == get_numvars()
+    zI = zero(Interval{T})
+    zbox = IntervalBox(zI, Val(N))
+    symIbox = IntervalBox(Interval{T}(-1 .. 1), Val(N))
+    q0 = IntervalBox(qq0)
+    t = t0 + Taylor1(orderT)
+
+    tv = Array{T}(undef, maxsteps+1)
+    xv = Array{IntervalBox{N, T}}(undef, maxsteps+1)
+    xTM1v = Array{TaylorModel1{TaylorN{T}, T}}(undef, dof, maxsteps+1)
+    x = Array{Taylor1{TaylorN{T}}}(undef, dof)
+    dx = Array{Taylor1{TaylorN{T}}}(undef, dof)
+    xaux = Array{Taylor1{TaylorN{T}}}(undef, dof)
+    xTMN = Array{TaylorModelN{N, T, T}}(undef, dof)
+    dxTM1 = Array{TaylorModel1{TaylorN{T}, T}}(undef, dof)
+    xTM1 = Array{TaylorModel1{TaylorN{T}, T}}(undef, dof)
+    xTM1K = Array{TaylorModel1{TaylorN{T}, T}}(undef, dof)
+
+    rem = Array{Interval{T}}(undef, dof)
+
+    @inbounds for i in eachindex(x)
+        qaux = normalize_taylor(qq0[i] + TaylorN(i, order=orderQ), δq0, true)
+        x[i] = Taylor1(qaux, orderT)
+        dx[i] = x[i]
+        xTMN[i] = TaylorModelN(qaux, zI, zbox, symIbox)
+        rem[i] = zI
+        xTM1v[i, 1] = TaylorModel1(deepcopy(x[i]), zI, zI, zI)
+    end
+
+    sign_tstep = copysign(1, tf - t0)
+    nsteps = 1
+    @inbounds tv[1] = t0
+    @inbounds xv[1] = evaluate(xTMN, symIbox)
+
+    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
+
+
+    while t0 * sign_tstep < tf * sign_tstep
+        δt = TaylorIntegration.taylorstep!(f!, t, x, dx, xaux, abstol, params, parse_eqs)
+        f!(dx, x, params, t)
+
+        δt = min(δt, sign_tstep*(tf-t0))
+        δt = sign_tstep * δt
+
+        @inbounds for i in eachindex(x)
+            dom = sign_tstep > 0 ? 0 .. δt : δt .. 0
+            x0 = sign_tstep > 0 ? dom.lo : dom.hi
+            Δ = zero(Interval{Float64})
+            xTM1[i] = TaylorModel1(deepcopy(x[i]), Δ, x0, dom)
+        end
+
+        # to reuse the previous TaylorModel and save some allocations
+        _validate_step!(xTM1, f!, dxTM1, xTMN, params, t, symIbox, dof,
+                        maxsteps=validatesteps, ε=ε, δ=δ)
+        t0 += δt
+        nsteps += 1
+
+        @inbounds t[0] = t0
+        @inbounds tv[nsteps] = t0
+        xv[nsteps] = evaluate(xTMN, symIbox)
+
+        @inbounds for i in eachindex(x)
+            aux_pol = evaluate(xTM1[i], Interval(δt))
+            rem[i] = remainder(xTM1[i])
+            xTMN[i] = fp_rpa(TaylorModelN(deepcopy(aux_pol), 0 .. 0, zbox, symIbox))
+            # temporal solution
+            j = 0
+            while (j < absorb_steps) && (mag(rem[i]) > 1.0e-10)
+                j += 1
+                xTMN[i] = absorb_remainder(xTMN[i])
+                rem[i] = remainder(xTMN[i])
+            end
+            x[i] = Taylor1(xTMN[i].pol, orderT)
+            xTM1v[i, nsteps] = xTM1[i]
+        end
+
+        if nsteps > maxsteps
+            @warn("""
+            Maximum number of integration steps reached; exiting.
+            """)
+            break
+        end
+    end
+
+    return view(tv, 1:nsteps), view(xv, 1:nsteps), view(xTM1v, :, 1:nsteps)
+end

test/validated_integ.jl

@@ -33,7 +33,7 @@ interval_rand(X::IntervalBox) = interval_rand.(X)
         ξ = set_variables("ξₓ ξᵥ", order=2*orderQ, numvars=length(q0))
         normalized_box = IntervalBox(-1 .. 1, Val(orderQ))
 
-        tTM, qv, qTM = validated_integ(falling_ball!, q0, δq0,
+        tTM, qv, qTM = validated_integ2(falling_ball!, q0, δq0,
             tini, tend, orderQ, orderT, abstol)
 
         @test length(qv) == length(qTM[1, :]) == length(tTM)
@@ -47,6 +47,52 @@ interval_rand(X::IntervalBox) = interval_rand.(X)
                 @test all(exactsol(tTM[n-1]+δt, tini, q0 .+ q0ξ) .∈ q)
             end
         end
+
+        tTM, qv, qTM = validated_integ2(falling_ball!, q0, δq0, tini, tend,
+                                        orderQ, orderT, abstol)
+
+        @test length(qv) == length(qTM[1, :]) == length(tTM)
+        @test length(tTM) < 501
+
+        for n = 2:length(tTM)
+            for it = 1:10
+                δt = interval_rand(domain(qTM[1, n]))
+                q0ξ = interval_rand(δq0)
+                q = evaluate.(evaluate(qTM[:, n], δt), (normalized_box,))
+                @test all(exactsol(tTM[n-1] + δt, tini, q0 .+ q0ξ) .∈ q)
+            end
+        end
+
+        @taylorize function x_square(dx, x, p, t)
+            dx[1] = x[1]^2
+            nothing
+        end
+
+        exactsol(t, x0) = x0[1] / (1 - x0[1] * t)
+
+        tini, tend = 0., 0.45
+        abstol = 1e-15
+        orderQ = 5
+        orderT = 20
+        q0 = [2.]
+        δq0 = IntervalBox(-0.1 .. 0.1, Val(1))
+        ξ = set_variables("ξₓ", numvars=1, order=2*orderQ)
+        # tTM, qv, qTM = validated_integ(x_square, q0, δq0, tini, tend, orderQ, orderT, abstol)
+        tTM, qv, qTM = validated_integ2(x_square, q0, δq0, tini, tend, orderQ, orderT, abstol)
+
+        @test length(qv) == length(qTM[1, :]) == length(tTM)
+        @test length(tTM) < 501
+        normalized_box = IntervalBox(-1 .. 1, Val(1))
+
+        for n = 2:length(tTM)
+            for it = 1:10
+                δt = interval_rand(domain(qTM[1, n]))
+                q0ξ = interval_rand(δq0)
+                q = evaluate.(evaluate.(qTM[:, n], δt), (normalized_box,))
+                @test all(exactsol(tTM[n-1]+δt, q0 .+ q0ξ) .∈ q)
+            end
+        end
+
     end
 
     @testset "Backward integration" begin
@@ -69,7 +115,7 @@ interval_rand(X::IntervalBox) = interval_rand.(X)
         ξ = set_variables("ξₓ ξᵥ", order=2*orderQ, numvars=length(q0))
         normalized_box = IntervalBox(-1 .. 1, Val(orderQ))
 
-        tTM, qv, qTM = validated_integ(falling_ball!, q0, δq0,
+        tTM, qv, qTM = validated_integ2(falling_ball!, q0, δq0,
             tini, tend, orderQ, orderT, abstol)
 
         @test length(qv) == length(qTM[1, :]) == length(tTM)