Add methods to output local Taylor polynomials for dense output (#134)

* Add methods to output local Taylor polynomial sols (dense output)

* Add DiffEqBase to extras section of Project.toml

* Update docs.yml

* Fix broken test

* Update packages compatibility and bump patch version

.github/workflows/docs.yml

@@ -25,4 +25,5 @@ jobs:
         env:
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # For authentication with GitHub Actions token
           DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # For authentication with SSH deploy key
+          GKSwstype: nul
         run: julia --project=docs/ docs/make.jl

Project.toml

@@ -1,10 +1,11 @@
 name = "TaylorIntegration"
 uuid = "92b13dbe-c966-51a2-8445-caca9f8a7d42"
 repo = "https://github.com/PerezHz/TaylorIntegration.jl.git"
-version = "0.8.10"
+version = "0.8.11"
 
 [deps]
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
+Elliptic = "b305315f-e792-5b7a-8f41-49f472929428"
 Espresso = "6912e4f1-e036-58b0-9138-08d1e6358ea9"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -20,15 +21,16 @@ TaylorSeries = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea"
 DiffEqBase = "6"
 Elliptic = "0.5, 1.0"
 Espresso = "0.6.1"
-OrdinaryDiffEq = "5"
-RecursiveArrayTools = "2.2"
-Reexport = "0.2, 1.0"
-Requires = "0.5.2, 1.0"
-StaticArrays = "0.12.5, 1.0"
-TaylorSeries = "0.10.9, 0.11"
+OrdinaryDiffEq = "5, 6"
+RecursiveArrayTools = "2"
+Reexport = "0.2, 1"
+Requires = "0.5.2, 1"
+StaticArrays = "0.12.5, 1"
+TaylorSeries = "0.12"
 julia = "1"
 
 [extras]
+DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 Elliptic = "b305315f-e792-5b7a-8f41-49f472929428"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -38,4 +40,5 @@ TaylorSeries = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["TaylorSeries", "LinearAlgebra", "Test", "Elliptic", "InteractiveUtils", "Pkg", "StaticArrays"]
+test = ["TaylorSeries", "LinearAlgebra", "Test", "Elliptic", "InteractiveUtils",
+    "Pkg", "DiffEqBase", "StaticArrays"]

src/explicitode.jl

@@ -286,6 +286,7 @@ end
 # taylorinteg
 @doc doc"""
     taylorinteg(f, x0, t0, tmax, order, abstol, params[=nothing]; kwargs... )
+    taylorinteg(f, x0, t0, tmax, order, abstol, Val(true), params[=nothing]; kwargs... )
 
 General-purpose Taylor integrator for the explicit ODE ``\dot{x}=f(x, p, t)``.
 The initial condition are specified by `x0` at time `t0`, and any parameters
@@ -299,8 +300,9 @@ convention of `DifferentialEquations.jl` to define this function, i.e.,
 
 It returns a vector with the values of time (independent variable),
 and a vector with the computed values of
-the dependent variable(s). The integration stops when time
-is larger than `tmax`, in which case the last returned
+the dependent variable(s), and if the method used involves `Val(true)` it also
+outputs the Taylor polynomial solutions obtained at each time step.
+The integration stops when time is larger than `tmax`, in which case the last returned
 value(s) correspond to that time, or when the number of saved steps is larger
 than `maxsteps`.
 
@@ -343,6 +345,8 @@ function f!(dx, x, p, t)
 end
 
 tv, xv = taylorinteg(f!, [3, 3], 0.0, 0.3, 25, 1.0e-20, maxsteps=100 )
+
+tv, xv, polynv = taylorinteg(f!, [3, 3], 0.0, 0.3, 25, 1.0e-20, maxsteps=100, Val(true) )
 ```
 
 """
@@ -391,6 +395,54 @@ function taylorinteg(f, x0::U, t0::T, tmax::T, order::Int, abstol::T,
     return view(tv,1:nsteps), view(xv,1:nsteps)
 end
 
+function taylorinteg(f, x0::U, t0::T, tmax::T, order::Int, abstol::T,
+        ::Val{true}, params = nothing;
+        maxsteps::Int=500, parse_eqs::Bool=true) where {T<:Real, U<:Number}
+
+    # Allocation
+    tv = Array{T}(undef, maxsteps+1)
+    xv = Array{U}(undef, maxsteps+1)
+    polynV = Array{Taylor1{U}}(undef, maxsteps+1)
+
+    # Initialize the Taylor1 expansions
+    t = Taylor1( T, order )
+    x = Taylor1( x0, order )
+
+    # Initial conditions
+    nsteps = 1
+    @inbounds t[0] = t0
+    @inbounds tv[1] = t0
+    @inbounds xv[1] = x0
+    @inbounds polynV[1] = deepcopy(x)
+    sign_tstep = copysign(1, tmax-t0)
+
+    # Determine if specialized jetcoeffs! method exists
+    parse_eqs = _determine_parsing!(parse_eqs, f, t, x, params)
+
+    # Integration
+    while sign_tstep*t0 < sign_tstep*tmax
+        δt = taylorstep!(f, t, x, abstol, params, parse_eqs) # δt is positive!
+        # Below, δt has the proper sign according to the direction of the integration
+        δt = sign_tstep * min(δt, sign_tstep*(tmax-t0))
+        x0 = evaluate(x, δt) # new initial condition
+        @inbounds polynV[nsteps+1] = deepcopy(x)
+        @inbounds x[0] = x0
+        t0 += δt
+        @inbounds t[0] = t0
+        nsteps += 1
+        @inbounds tv[nsteps] = t0
+        @inbounds xv[nsteps] = x0
+        if nsteps > maxsteps
+            @warn("""
+            Maximum number of integration steps reached; exiting.
+            """)
+            break
+        end
+    end
+
+    return view(tv,1:nsteps), view(xv,1:nsteps), view(polynV, 1:nsteps)
+end
+
 function taylorinteg(f!, q0::Array{U,1}, t0::T, tmax::T, order::Int, abstol::T,
         params = nothing; maxsteps::Int=500, parse_eqs::Bool=true) where {T<:Real, U<:Number}
 
@@ -447,6 +499,70 @@ function taylorinteg(f!, q0::Array{U,1}, t0::T, tmax::T, order::Int, abstol::T,
     return view(tv,1:nsteps), view(transpose(view(xv,:,1:nsteps)),1:nsteps,:)
 end
 
+# Method to also return (output) the Taylor polynomials of the integration
+function taylorinteg(f!, q0::Array{U,1}, t0::T, tmax::T, order::Int, abstol::T,
+        ::Val{true}, params = nothing;
+        maxsteps::Int=500, parse_eqs::Bool=true) where {T<:Real, U<:Number}
+
+    # Allocation
+    tv = Array{T}(undef, maxsteps+1)
+    dof = length(q0)
+    xv = Array{U}(undef, dof, maxsteps+1)
+    polynV = Array{Taylor1{U}}(undef, dof, maxsteps+1)
+
+    # Initialize the vector of Taylor1 expansions
+    t = Taylor1(T, order)
+    x = Array{Taylor1{U}}(undef, dof)
+    dx = Array{Taylor1{U}}(undef, dof)
+    xaux = Array{Taylor1{U}}(undef, dof)
+    for i in eachindex(q0)
+        @inbounds x[i] = Taylor1( q0[i], order )
+        @inbounds dx[i] = Taylor1( zero(q0[i]), order )
+    end
+
+    # Initial conditions
+    @inbounds t[0] = t0
+    x .= Taylor1.(q0, order)
+    x0 = deepcopy(q0)
+    @inbounds tv[1] = t0
+    @inbounds xv[:,1] .= q0
+    @inbounds polynV[:,1] .= deepcopy.(x)
+    sign_tstep = copysign(1, tmax-t0)
+
+    # Determine if specialized jetcoeffs! method exists
+    parse_eqs = _determine_parsing!(parse_eqs, f!, t, x, dx, params)
+
+    # Integration
+    nsteps = 1
+    while sign_tstep*t0 < sign_tstep*tmax
+        δt = taylorstep!(f!, t, x, dx, xaux, abstol, params, parse_eqs) # δt is positive!
+        # Below, δt has the proper sign according to the direction of the integration
+        δt = sign_tstep * min(δt, sign_tstep*(tmax-t0))
+        evaluate!(x, δt, x0) # new initial condition
+        # Store the Taylor polynomial solution
+        @inbounds polynV[:,nsteps+1] .= deepcopy.(x)
+        for i in eachindex(x0)
+            @inbounds x[i][0] = x0[i]
+            @inbounds dx[i] = Taylor1( zero(x0[i]), order )
+        end
+        t0 += δt
+        @inbounds t[0] = t0
+        nsteps += 1
+        @inbounds tv[nsteps] = t0
+        @inbounds xv[:,nsteps] .= x0
+        if nsteps > maxsteps
+            @warn("""
+            Maximum number of integration steps reached; exiting.
+            """)
+            break
+        end
+    end
+
+    return view(tv,1:nsteps), view(transpose(view(xv,:,1:nsteps)),1:nsteps,:),
+        view(transpose(view(polynV, :, 1:nsteps)), 1:nsteps, :)
+end
+
+
 # Integrate and return results evaluated at given time
 @doc doc"""
     taylorinteg(f, x0, trange, order, abstol, params[=nothing]; keyword... )
@@ -660,5 +776,31 @@ for R in (:Number, :Integer)
                 taylorinteg(f, q0_, vec(trange.*one(t0)), order, abstol,
                     params, maxsteps=maxsteps, parse_eqs=parse_eqs)
         end
+
+        function taylorinteg(f, xx0::S, tt0::T, ttmax::U, order::Int, aabstol::V,
+                ::Val{true}, params = nothing; maxsteps::Int=500, parse_eqs::Bool=true) where
+                    {S<:$R, T<:Real, U<:Real, V<:Real}
+
+            # In order to handle mixed input types, we promote types before integrating:
+            t0, tmax, abstol, bfloat = promote(tt0, ttmax, aabstol, one(Float64))
+            x0, tfloat1 = promote(xx0, t0)
+
+            return taylorinteg(f, x0, t0, tmax, order, abstol, Val(true), params,
+                maxsteps=maxsteps, parse_eqs=parse_eqs)
+        end
+
+        function taylorinteg(f, q0::Array{S,1}, tt0::T, ttmax::U, order::Int, aabstol::V,
+                ::Val{true}, params = nothing; maxsteps::Int=500, parse_eqs::Bool=true) where
+                    {S<:$R, T<:Real, U<:Real, V<:Real}
+
+            #promote to common type before integrating:
+            t0, tmax, abstol, afloat = promote(tt0, ttmax, aabstol, one(Float64))
+            elq0, tt0 = promote(q0[1], t0)
+            #convert the elements of q0 to the common, promoted type:
+            q0_ = convert(Array{typeof(elq0)}, q0)
+
+            return taylorinteg(f, q0_, t0, tmax, order, abstol, Val(true), params,
+                maxsteps=maxsteps, parse_eqs=parse_eqs)
+        end
     end
 end

test/many_ode.jl

@@ -59,6 +59,14 @@ using LinearAlgebra: Transpose, norm
         @test (isnan(xv2[end,1]) && isnan(xv2[end,2]))
         @test abs(xv2[3,2] - 4/3) ≤ eps(4/3)
         @test abs(xv2[2,1] - 4.8) ≤ eps(4.8)
+
+        # Output includes Taylor polynomial solution
+        tv, xv, polynV = taylorinteg(eqs_mov!, q0, 0, 0.5, _order, _abstol, Val(true), nothing, maxsteps=2)
+        @test size(polynV) == (3, 2)
+        @test xv[1,:] == q0
+        @test polynV[1,:] == Taylor1.(q0, _order)
+        @test xv[2,:] == evaluate.(polynV[2, :], tv[2]-tv[1])
+        @test polynV[2,2] == Taylor1(ones(_order+1))
     end
 
     @testset "Tests: dot{x}=x.^2, x(0) = [3, 3]" begin
@@ -98,8 +106,15 @@ using LinearAlgebra: Transpose, norm
         @test tv[end] < 1/3
         @test tv[end] == tmax
         @test xv[end,1] == xv[end,2]
-        @test abs(xv[end,1]-exactsol(tv[end], xv[1,1])) < 1e-11
-        @test abs(xv[end,2]-exactsol(tv[end], xv[1,2])) < 1e-11
+        @test abs(xv[end,1]-exactsol(tv[end], xv[1,1])) < 1.0e-11
+        @test abs(xv[end,2]-exactsol(tv[end], xv[1,2])) < 1.0e-11
+
+        # Output includes Taylor polynomial solution
+        tv, xv, polynV = taylorinteg(eqs_mov!, q0, 0, tmax, _order, _abstol, Val(true), nothing)
+        @test size(polynV) == size(xv)
+        @test polynV[1,:] == Taylor1.(q0, _order)
+        @test xv[end,1] == evaluate.(polynV[end, 1], tv[end]-tv[end-1])
+        @test polynV[:,1] == polynV[:,2]
     end
 
     @testset "Test non-autonomous ODE (2): dot{x}=cos(t)" begin

test/one_ode.jl

@@ -68,6 +68,14 @@ using LinearAlgebra: norm
         @test xv2[1] == x0
         @test isnan(xv2[end])
         @test abs(xv2[5] - 2.0) ≤ eps(2.0)
+
+        # Output includes Taylor polynomial solution
+        tv, xv, polynV = taylorinteg(eqs_mov, x0, 0, 0.5, _order, _abstol, Val(true), maxsteps=2)
+        @test length(polynV) == 3
+        @test xv[1] == x0
+        @test polynV[1] == Taylor1(x0, _order)
+        @test xv[2] == evaluate(polynV[2], tv[2]-tv[1])
+        @test polynV[2] == Taylor1(ones(_order+1))
     end
 
     @testset "Tests: dot{x}=x^2, x(0) = 3; nsteps <= maxsteps" begin

test/taylorize.jl

@@ -100,6 +100,14 @@ import Pkg
         @test iszero( norm(tv2-tv2e, Inf) )
         @test iszero( norm(xv2-xv2e, Inf) )
         @test_throws UndefVarError TaylorIntegration.__jetcoeffs!(Val(true), xdot2_err, tT, qT, similar(qT), similar(qT), [1.0])
+
+        # Output includes Taylor polynomial solution
+        tv3t, xv3t, polynV3t = taylorinteg(xdot3, x0, t0, tf, _order, _abstol, Val(true), 0.0, maxsteps=2)
+        @test length(polynV3t) == 3
+        @test xv3t[1] == x0
+        @test polynV3t[1] == Taylor1(x0, _order)
+        @test xv3t[2] == evaluate(polynV3t[2], tv3t[2]-tv3t[1])
+        @test polynV3t[2] == Taylor1([(-1.0)^i for i=0:_order])
     end