perezhz / taylorintegration.jl Goto Github PK

@JuliaRegistrator register

@JuliaRegistrator register

There is a subtle problem when using the macro in functions involving ^r (for r different from 0, 1, 2). I illustrate the problem here with the integration of the Duffing equation, using of the macro (using current master of TaylorSeries):

julia> using TaylorIntegration

julia> @taylorize function duffing!(du, u, params, t)
           du[1] = u[2]
           du[2] = u[1]-0.01*u[1]^3
           nothing
       end

julia> x = taylorinteg(duffing!, [0.0, 1.0], 0.0:0.25:20.0, 20, 1.e-20);

julia> x[2,:]
2-element Array{Float64,1}:
 NaN
 NaN

The problem is related with TaylorSeries.pow! and the fact that it requires one parameter (the order of the first non-zero entry of a polynomial), which by default is set to zero. That default is used by @taylorize, which produces the NaNs in this case, and it does not appear if one computes c = a^3.

julia> a = 0.0*Taylor1(20)   # Taylor polynomial identical to zero
 0.0 + 𝒪(t²¹)

julia> c = Taylor1(20)
 1.0 t + 𝒪(t²¹)

# We compute c=a^3 term by term, using pow!
julia> TaylorSeries.pow!(c, a, 3, 0)

julia> TaylorSeries.pow!(c, a, 3, 1)

julia> c
 NaN t + 𝒪(t²¹)

The problem is solved by JuliaDiff/TaylorSeries.jl#236; we will require a new tag of TaylorSeries.jl and related updates here.

using DifferentialEquations
using TaylorIntegration
@taylorize function benchmark7!(du, u, p, t)
    du[1] = u[3]^3 - u[2] + u[4]
    du[2] = u[3]
    du[3] = 2
    du[4] = u[4]
    return nothing
end
tspan = (0.0, 10.0)
q0 = [0.4, 0.5, 0.3, 0]
prob = ODEProblem(benchmark7!, q0, tspan)
sol = solve(prob, TaylorMethod(10), abstol=1e-20);

UndefVarError: false not defined

Stacktrace:
 [1] macro expansion at ./logging.jl:322 [inlined]
 [2] _determine_parsing!(::Bool, ::Function, ::Taylor1{Float64}, ::Array{Taylor1{Float64},1}, ::Array{Taylor1{Float64},1}, ::DiffEqBase.NullParameters) at /home/sguadalupe/.julia/packages/TaylorIntegration/KHxcX/src/explicitode.jl:276
 [3] taylorinteg(::Function, ::Array{Float64,1}, ::Float64, ::Float64, ::Int64, ::Float64, ::DiffEqBase.NullParameters; maxsteps::Int64, parse_eqs::Bool) at /home/sguadalupe/.julia/packages/TaylorIntegration/KHxcX/src/explicitode.jl:421
 [4] solve(::ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,DiffEqBase.NullParameters,ODEFunction{true,typeof(benchmark7!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}, ::TaylorMethod, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}; verbose::Bool, saveat::Array{Float64,1}, abstol::Float64, save_start::Bool, save_end::Bool, timeseries_errors::Bool, maxiters::Int64, callback::Nothing, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /home/sguadalupe/.julia/packages/TaylorIntegration/KHxcX/src/common.jl:107
 [5] top-level scope at In[15]:1

There is a workaround using 2.0+zero(u[1]) instead of only 2 for du[3]

This is necessary to spot easily performance regressions.

I am having trouble calculating a Poincare section for a positive momentum. I tried with the example you provided here. The specific portion of the code that gives me problems is:

function g(t, x, dx)
    px = constant_term(x[3])
    
    if px > zero(px)
        # ...return x
        return x[1]
    else
        #otherwise, discard the crossing
        return nothing
    end
end

Since, I got the error:

MethodError: no method matching surfacecrossing(::Taylor1{Float64}, ::Nothing, ::Int64)
Closest candidates are:
  surfacecrossing(::Taylor1{T<:Number}, !Matched::Taylor1{T<:Number}, ::Int64) where T<:Number at /Users/porras/.julia/dev/TaylorIntegration/src/rootfinding.jl:13

On current master, taylorinteg breaks down when using BigFloats, if the maxsteps keyword argument is greater than the actual number of steps taken during a given integration:

julia> using TaylorIntegration

julia> function harmosc!(t, x, dx) #the harmonic oscillator ODE
           dx[1] = x[2]
           dx[2] = -x[1]
           nothing
       end
harmosc! (generic function with 1 method)

julia> @time tv, xv = taylorinteg(harmosc!, [0, -1], 0, 20BigFloat(pi), 90, 1.0e-80, maxsteps=15 );
WARNING: Maximum number of integration steps reached; exiting.
  0.663779 seconds (452.26 k allocations: 25.480 MiB, 9.63% gc time)

julia> tv[end] == 20BigFloat(pi)
true

julia> @time tv, xv = taylorinteg(harmosc!, [0, -1], 0, 20BigFloat(pi), 90, 1.0e-80, maxsteps=16 );
ERROR: UndefRefError: access to undefined reference
Stacktrace:
 [1] transpose_f!(::Base.#transpose, ::Array{BigFloat,2}, ::Array{BigFloat,2}) at ./linalg/transpose.jl:54
 [2] transpose(::Array{BigFloat,2}) at ./linalg/transpose.jl:121
 [3] #taylorinteg#4(::Int64, ::Function, ::Function, ::Array{BigFloat,1}, ::BigFloat, ::BigFloat, ::Int64, ::BigFloat) at /Users/Jorge/TaylorIntegration.jl/src/explicitode.jl:453
 [4] (::TaylorIntegration.#kw##taylorinteg)(::Array{Any,1}, ::TaylorIntegration.#taylorinteg, ::Function, ::Array{BigFloat,1}, ::BigFloat, ::BigFloat, ::Int64, ::BigFloat) at ./<missing>:0
 [5] #taylorinteg#3(::Int64, ::Function, ::Function, ::Array{Int64,1}, ::Int64, ::BigFloat, ::Int64, ::Float64) at /Users/Jorge/TaylorIntegration.jl/src/explicitode.jl:408
 [6] (::TaylorIntegration.#kw##taylorinteg)(::Array{Any,1}, ::TaylorIntegration.#taylorinteg, ::Function, ::Array{Int64,1}, ::Int64, ::BigFloat, ::Int64, ::Float64) at ./<missing>:0

In the example above, if maxsteps=15 or less, then the answer is returned correctly, but if maxsteps is greater than 15, then it just breaks down. I think this is related with how transpose handles matrices with #undef components, when the elements of the matrix are BigFloats. When working with Float64, then the problem doesn't occur:

julia> @time tv, xv = taylorinteg(harmosc!, [0, -1], 0, 20pi, 25, 1.0e-20, maxsteps=200 );
  0.359674 seconds (200.31 k allocations: 10.576 MiB)

I'm having some trouble with jetcoeffs! and its relationship with TaylorSeries.derivative. Here is some code illustrating my confusion:

ord = 5
u_0 = [0.,1.]
T = typeof(u_0[1])
t =Taylor1(T,ord)
dof = length(u_0)
x = Array{Taylor1{T}}(undef,dof)
dx = Array{Taylor1{T}}(undef,dof)
xaux = Array{Taylor1{typeof(u_0[1])}}(undef,dof)
for i in eachindex(u_0)
    x[i] = Taylor1(u_0[i],ord)
    dx[i] = Taylor1(zero(u_0[i]),ord)
end

# Dynamics Equations 
f(t,u) = [1.,2.] .* [t].^5
f(t,u,du) = (du .= f(t,u)) 

# compute higher order coeffs into x
jetcoeffs!(f,t,x,dx,xaux)
x

# Test Derivatives
TaylorSeries.derivative.(f(t,x),5)

TaylorSeries.derivative.(x,5)

This is a simple situation where the dynamics only depend on t. This was chosen as an illustrative example because the 5th order derivative will just be [1.,2.] .* factorial(5).

I am confused because it seems that the jetcoeffs! computation is not working in this example. What I expected was that jetcoeffs! computes the higher order coefficients into x such that it's higher order derivatives correspond to those given by the dynamics function f.

However, jetcoeffs!(f,t,x,dx,xaux) is not computing the higher order coefficients into x.
As such, the derivative given by looking at the higher coefficients does not equal the derivative obtained by using TaylorSeries.derivative on the function to that order. If the jet coefficients were calculated correctly these would be equivalent, correct?

I suspect that I am confusing something here and appreciate some guidance. This confusion may be in how I'm understanding the independent variable.

E.g. how would this change if I wanted to evaluate for a time other than t=0.? Would t=1 correspond to t= Taylor1(T,ord) .+ 1.?

In the case of t=1. the jet coefficients are computed correctly, i.e.
TaylorSeries.derivative.(f(t,x),5) = TaylorSeries.derivative.(x,5) = [1.,2.] .* factorial(5)

However, in the case of t=2. (i.e. t =Taylor1(T,ord) .+ 2.) the jet coefficients are computed such that the final order derivative is double what it should be:
TaylorSeries.derivative.(f(t,x),5) = [1.,2.] .* factorial(5) still, as desired.
but TaylorSeries.derivative.(x,5) = 2 * TaylorSeries.derivative.(f(t,x),5)

For this reason I suspect I am incorrectly using the independent variable when computing jet coefficients.

@JuliaRegistrator register

On the upcoming 1.3 this package errors with

Complex dependent variable: Test Failed at /root/.julia/packages/TaylorIntegration/R6yrv/test/taylorize.jl:198
  Expression: norm(abs2(xv1p[end]) - abs2(exact_sol(tv1p[end], cc, cx0)), Inf) < 1.0e-15
   Evaluated: 2.1094237467877974e-15 < 1.0e-15
Stacktrace:
 [1] top-level scope at /root/.julia/packages/TaylorIntegration/R6yrv/test/taylorize.jl:198
 [2] top-level scope at /workspace/srcdir/julia/usr/share/julia/stdlib/v1.3/Test/src/Test.jl:1107
 [3] top-level scope at /root/.julia/packages/TaylorIntegration/R6yrv/test/taylorize.jl:153
Test Summary:              | Pass  Fail  Total
Complex dependent variable |   17     1     18

It seems that the test tolerance is put a bit too tightly?

@JuliaRegistrator register

I have a trouble with my installation. I tried to install TaylorIntegration package and got an error. Explicitly:

pkg> add TaylorIntegration Resolving package versions... ERROR: Unsatisfiable requirements detected for package DiffEqBase [2b5f629d]: DiffEqBase [2b5f629d] log: ├─possible versions are: [3.13.2-3.13.3, 4.0.0-4.0.1, 4.1.0, 4.2.0, 4.3.0-4.3.1, 4.4.0, 4.5.0, 4.6.0, 4.7.0, 4.8.0, 4.9.0, 4.10.0-4.10.1, 4.11.0-4.11.1, 4.12.0, 4.13.0, 4.14.0-4.14.1, 4.15.0, 4.16.0, 4.17.0, 4.18.0, 4.19.0, 4.20.0-4.20.3, 4.21.0, 4.21.2-4.21.3, 4.22.0-4.22.2, 4.23.0, 4.23.2-4.23.4, 4.24.0-4.24.3, 4.25.0-4.25.1, 4.26.0-4.26.3, 4.27.0-4.27.1, 4.28.0-4.28.1, 4.29.0-4.29.2, 4.30.0-4.30.2, 4.31.0-4.31.2, 4.32.0, 5.0.0-5.0.1, 5.1.0, 5.2.0-5.2.3, 5.3.0-5.3.2, 5.4.0-5.4.1, 5.5.0-5.5.2, 5.6.0-5.6.4, 5.7.0, 5.8.0-5.8.1, 5.9.0, 5.10.0-5.10.3, 5.11.0-5.11.1, 5.12.0, 5.13.0, 5.14.0-5.14.2, 5.15.0, 5.16.0-5.16.5, 5.17.0-5.17.1, 5.18.0, 5.19.0, 5.20.0-5.20.1, 6.0.0, 6.1.0, 6.2.0-6.2.4, 6.3.0-6.3.6, 6.4.0-6.4.2, 6.5.0-6.5.1, 6.6.0, 6.7.0, 6.8.0, 6.9.0-6.9.4, 6.10.0-6.10.2, 6.11.0, 6.12.0-6.12.5, 6.13.0-6.13.3, 6.14.0-6.14.2, 6.15.0-6.15.2, 6.16.0, 6.17.0-6.17.3, 6.18.0-6.18.1, 6.19.0, 6.20.0, 6.21.0-6.21.1, 6.22.0-6.22.2, 6.23.0, 6.24.0, 6.25.0-6.25.2, 6.26.0, 6.27.0, 6.28.0, 6.29.0-6.29.3, 6.30.0-6.30.4, 6.31.0-6.31.1, 6.32.0-6.32.2, 6.33.0-6.33.1, 6.34.0-6.34.3, 6.35.0-6.35.2, 6.36.0-6.36.4, 6.37.0, 6.38.0-6.38.4, 6.39.0-6.39.1, 6.40.0-6.40.9, 6.41.0-6.41.3, 6.42.0, 6.43.0-6.43.1, 6.44.0-6.44.3, 6.45.0-6.45.1, 6.46.0-6.46.1, 6.47.0-6.47.1, 6.48.0] or uninstalled ├─restricted by compatibility requirements with TaylorIntegration [92b13dbe] to versions: [4.21.3, 4.22.0-4.22.2, 4.23.0, 4.23.2-4.23.4, 4.24.0-4.24.3, 4.25.0-4.25.1, 4.26.0-4.26.3, 4.27.0-4.27.1, 4.28.0-4.28.1, 4.29.0-4.29.2, 4.30.0-4.30.2, 4.31.0-4.31.2, 4.32.0, 5.0.0-5.0.1, 5.1.0, 5.2.0-5.2.3, 5.3.0-5.3.2, 5.4.0-5.4.1, 5.5.0-5.5.2, 5.6.0-5.6.4, 5.7.0, 5.8.0-5.8.1, 5.9.0, 5.10.0-5.10.3, 5.11.0-5.11.1, 5.12.0, 5.13.0, 5.14.0-5.14.2, 5.15.0, 5.16.0-5.16.5, 5.17.0-5.17.1, 5.18.0, 5.19.0, 5.20.0-5.20.1, 6.0.0, 6.1.0, 6.2.0-6.2.4, 6.3.0-6.3.6, 6.4.0-6.4.2, 6.5.0-6.5.1, 6.6.0, 6.7.0, 6.8.0, 6.9.0-6.9.4, 6.10.0-6.10.2, 6.11.0, 6.12.0-6.12.5, 6.13.0-6.13.3, 6.14.0-6.14.2, 6.15.0-6.15.2, 6.16.0, 6.17.0-6.17.3, 6.18.0-6.18.1, 6.19.0, 6.20.0, 6.21.0-6.21.1, 6.22.0-6.22.2, 6.23.0, 6.24.0, 6.25.0-6.25.2, 6.26.0, 6.27.0, 6.28.0, 6.29.0-6.29.3, 6.30.0-6.30.4, 6.31.0-6.31.1, 6.32.0-6.32.2, 6.33.0-6.33.1, 6.34.0-6.34.3, 6.35.0-6.35.2, 6.36.0-6.36.4, 6.37.0, 6.38.0-6.38.4, 6.39.0-6.39.1, 6.40.0-6.40.9, 6.41.0-6.41.3, 6.42.0, 6.43.0-6.43.1, 6.44.0-6.44.3, 6.45.0-6.45.1, 6.46.0-6.46.1, 6.47.0-6.47.1, 6.48.0] │ └─TaylorIntegration [92b13dbe] log: │ ├─possible versions are: [0.4.0-0.4.1, 0.5.0-0.5.1, 0.6.0-0.6.1, 0.7.0-0.7.1, 0.8.0-0.8.4] or uninstalled │ └─restricted to versions * by an explicit requirement, leaving only versions [0.4.0-0.4.1, 0.5.0-0.5.1, 0.6.0-0.6.1, 0.7.0-0.7.1, 0.8.0-0.8.4] ├─restricted by julia compatibility requirements to versions: [3.13.2-3.13.3, 4.0.0-4.0.1, 4.1.0, 4.2.0, 4.3.0-4.3.1, 4.4.0, 4.5.0, 4.6.0, 4.7.0, 4.8.0, 4.9.0, 4.10.0-4.10.1, 4.11.0-4.11.1, 4.12.0, 4.13.0, 4.14.0-4.14.1, 4.15.0, 4.16.0, 4.17.0, 4.18.0, 4.19.0, 4.20.0-4.20.3, 4.21.0, 4.21.2-4.21.3, 4.22.0-4.22.2, 4.23.0, 4.23.2-4.23.4, 4.24.0-4.24.3, 4.25.0-4.25.1, 4.26.0-4.26.3, 4.27.0-4.27.1, 4.28.0-4.28.1, 4.29.0-4.29.2, 4.30.0-4.30.2, 4.31.0-4.31.2, 4.32.0, 5.0.0-5.0.1, 5.1.0, 5.2.0-5.2.3, 5.3.0-5.3.2, 5.4.0-5.4.1, 5.5.0-5.5.2, 5.6.0-5.6.4, 5.7.0, 5.8.0-5.8.1, 5.9.0, 5.10.0-5.10.3, 5.11.0-5.11.1, 5.12.0, 5.13.0, 5.14.0-5.14.2, 5.15.0, 5.16.0-5.16.5, 5.17.0-5.17.1, 5.18.0, 5.19.0, 5.20.0-5.20.1, 6.0.0, 6.1.0, 6.2.0-6.2.4, 6.3.0-6.3.6, 6.4.0-6.4.2, 6.5.0-6.5.1, 6.6.0, 6.7.0, 6.8.0, 6.9.0-6.9.4, 6.10.0-6.10.2, 6.11.0, 6.12.0-6.12.5, 6.13.0-6.13.3, 6.14.0-6.14.2, 6.15.0-6.15.2, 6.16.0, 6.17.0-6.17.3, 6.18.0-6.18.1, 6.19.0, 6.20.0] or uninstalled, leaving only versions: [4.21.3, 4.22.0-4.22.2, 4.23.0, 4.23.2-4.23.4, 4.24.0-4.24.3, 4.25.0-4.25.1, 4.26.0-4.26.3, 4.27.0-4.27.1, 4.28.0-4.28.1, 4.29.0-4.29.2, 4.30.0-4.30.2, 4.31.0-4.31.2, 4.32.0, 5.0.0-5.0.1, 5.1.0, 5.2.0-5.2.3, 5.3.0-5.3.2, 5.4.0-5.4.1, 5.5.0-5.5.2, 5.6.0-5.6.4, 5.7.0, 5.8.0-5.8.1, 5.9.0, 5.10.0-5.10.3, 5.11.0-5.11.1, 5.12.0, 5.13.0, 5.14.0-5.14.2, 5.15.0, 5.16.0-5.16.5, 5.17.0-5.17.1, 5.18.0, 5.19.0, 5.20.0-5.20.1, 6.0.0, 6.1.0, 6.2.0-6.2.4, 6.3.0-6.3.6, 6.4.0-6.4.2, 6.5.0-6.5.1, 6.6.0, 6.7.0, 6.8.0, 6.9.0-6.9.4, 6.10.0-6.10.2, 6.11.0, 6.12.0-6.12.5, 6.13.0-6.13.3, 6.14.0-6.14.2, 6.15.0-6.15.2, 6.16.0, 6.17.0-6.17.3, 6.18.0-6.18.1, 6.19.0, 6.20.0] └─restricted by compatibility requirements with RecipesBase [3cdcf5f2] to versions: [6.26.0, 6.27.0, 6.28.0, 6.29.0-6.29.3, 6.30.0-6.30.4, 6.31.0-6.31.1, 6.32.0-6.32.2, 6.33.0-6.33.1, 6.34.0-6.34.3, 6.35.0-6.35.2, 6.36.0-6.36.4, 6.37.0, 6.38.0-6.38.4, 6.39.0-6.39.1, 6.40.0-6.40.9, 6.41.0-6.41.3, 6.42.0, 6.43.0-6.43.1, 6.44.0-6.44.3, 6.45.0-6.45.1, 6.46.0-6.46.1, 6.47.0-6.47.1, 6.48.0] or uninstalled — no versions left └─RecipesBase [3cdcf5f2] log: ├─possible versions are: [0.4.0, 0.5.0, 0.6.0, 0.7.0, 0.8.0, 1.0.0-1.0.2, 1.1.0] or uninstalled └─restricted by compatibility requirements with Plots [91a5bcdd] to versions: [1.0.0-1.0.2, 1.1.0] └─Plots [91a5bcdd] log: ├─possible versions are: [0.12.1-0.12.4, 0.13.0-0.13.1, 0.14.0-0.14.2, 0.15.0-0.15.1, 0.16.0, 0.17.0-0.17.4, 0.18.0, 0.19.0-0.19.3, 0.20.0-0.20.6, 0.21.0, 0.22.0-0.22.5, 0.23.0-0.23.2, 0.24.0, 0.25.0-0.25.3, 0.26.0-0.26.3, 0.27.0-0.27.1, 0.28.0-0.28.4, 0.29.0-0.29.9, 1.0.0-1.0.14, 1.1.0-1.1.4, 1.2.0-1.2.6, 1.3.0-1.3.7, 1.4.0-1.4.4, 1.5.0-1.5.9, 1.6.0-1.6.11] or uninstalled └─restricted to versions 1.4.3 by an explicit requirement, leaving only versions 1.4.3

I already attempted the with pkg> update and pkg> add DiffEqBase. But I got the same error. Any ideas on how to solve this issue?

When using @taylorize inside a module, a warning appears about incremental compilation "maybe" being "fatally broken". I wasn't able to reproduce the problem when defining a module in the REPL, but figured out the following MWE:

I created the MyPkg julia package, which uses @taylorize to define an ODE. Then, MyPkg may be git clone'd on an empty folder, and from MyPkg root dir start a julia REPL and do:

# first, in an empty folder do:
# git clone https://github.com/PerezHz/MyPkg
# cd MyPkg
# julia ### open a Julia REPL
# then, in the REPL do:
import Pkg
Pkg.activate(".") # activate environment on MyPkg root dir
using MyPkg # warning appears

The warning messages that appear are:

[ Info: Precompiling MyPkg [87ba5f82-8271-4dfc-8d8f-dcd99beeaf9d]
WARNING: eval into closed module TaylorIntegration:
:-
  ** incremental compilation may be fatally broken for this module **

WARNING: eval into closed module TaylorIntegration:
:^
  ** incremental compilation may be fatally broken for this module **

WARNING: eval into closed module TaylorIntegration:
typeof(MyPkg.xdot2)()
  ** incremental compilation may be fatally broken for this module **

WARNING: eval into closed module TaylorIntegration:
nothing
  ** incremental compilation may be fatally broken for this module **

So far, I've reproduced the error on julia 1.3 and 1.4 locally and on a remote machine.

The incremental compilation warnings appear the first time MyPkg is loaded via using MyPkg, and do not appear on subsequent loadings, until MyPkg is changed internally (e.g., by uncommenting the MyPkg.greet() function), and using MyPkg is ran again so that the code is re-compiled. Also, as far as I can tell, things seem to work fine (taylorinteg works fine, etc.) even when the warnings appear, but I don't know if it's safe to live with these warnings hanging around.

A workaround is to set __precompile__(false), but then MyPkg is compiled each time it is loaded, and compilation times may become exceptionally slow (~ 10 min) in some cases where the defined ODE involves many operations.

This issue seems to be related to the use of eval in the definition of @taylorize, and found a post in Discourse which dealt with similar things, but wasn't able to figure out a way to get rid of these warnings so that it's possible to not use __precompile__(false). @lbenet, do you know what might be going on here?

Hi! I'm opening this issue as a to-do list of current WIP that we would like to get done for the upcoming releases of TaylorIntegration.jl. If you think something should be modified or added, feel free to comment. As always, feedback is very welcome!

• Add documentation
• Add JuliaDiffEq common interface bindings for (#19)
• Migrate tests from FactCheck to Base.Test; separate tests into several files (see [#33])
• Solve other pending open issues (#20, #27)

@lbenet @blas-ko

EDIT: Updated the PR where the tests are being migrated.
EDIT: Updated title and first comment of this issue, so that it reflects current development efforts.
EDIT: Updated issue title
EDIT: Reviewing, testing, and merging of #31 will be done after v0.4.0 tagging
EDIT: documentation was added in #57; update to-do list accordingly
EDIT: After v0.4.0 release, #20 was closed

This issue is related to JuliaLang/julia#42271, which changes the values computed for the time-step in the most common situation (using Float64) since that involves Float64^Float64 computations. Therefore, it slightly changes the results when using Julia 1.7.x or 1.8, in a way breaking reproducibility.

You probably also want to enable CompatHelper again (here).

CompatHelper.jl helps to keep the [compat] entries up-to-date. It automatically makes PRs to update the Project.toml file once installed (instructions here).

When I run the following code in julia 0.4.7-pre+3:

using TaylorIntegration
f(t, x) = x.^2
const x0 = [3.0,3.0]
const t0 = 0.0
const tmax = 0.3
const order = 28
const abstol = 1e-20
@time tT, xT = taylorinteg(f, x0, t0, tmax, order, abstol, maxsteps=50);

xT becomes a Array{Float64,2} with size 51x2 even when taylorinteg performs only 12 integration steps. When taylorinteg reaches the maximum number of allowed integration steps by the user (i.e., when the number of performed steps is equal to maxsteps), there is no issue. But I consider that when taylorinteg finishes before reaching maxsteps steps, xT should have size (nsteps+1) x dof, where nsteps is the number of integration steps, and dof=length(x0). Something similar happens when typeof(x0) has type T<:Number (i.e., when the initial condition is a scalar).

In the example above, xT should have size 13x2 (13 and not 12, because we are also returning the initial condition as the "0-th" step). But I understand that the current output format was chosen because of performance reasons (see #2). So a solution I propose to eliminate this ambiguity while preserving performance is:

Consider the following julia code:

a = 20000; b=100 #the size of a 2-dimensional array
xT = Array{Float64}(a, b) #allocating a 2-dimensional array
nsteps = 17691 #a random integer less than `a`

#filling the array with random values, up to the nsteps+1 row:
for i in 1:nsteps+1, j in 1:b
    xT[i,j] = rand()
end

The problem is then to obtain a Array{Float64,2}, with size (nsteps+1) x 100 from xT, which currently has size 20000x100. And to do this without performance penalties. Two possible solutions to this are either

getindex(slice(xT,1:nsteps+1,:),:,:)

or

getindex(sub(xT,1:nsteps+1,:),:,:)

If we @time them, we obtain the following:

julia> @time xT = getindex(sub(xT,1:nsteps+1,:),:,:)
  0.006763 seconds (24 allocations: 13.499 MB, 40.11% gc time)
@time xT = getindex(slice(xT,1:nsteps+1,:),:,:)
  0.007079 seconds (26 allocations: 13.499 MB, 59.29% gc time)

In this case, sub appears to perform slightly better than slice. I understand that sub essentially points to values, while slice creates a copy. Something similar could be done for tT, and also for xT, when the initial condition is a T<:Number (i.e., when the initial condition is a scalar). I will test this and report back here.

What are your thoughts on this, @lbenet ?

For the context see this discussion. In short, integrating a function such as

@taylorize function bball!(t, x, dx)
    dx[1] = x[2]
    dx[2] = -one(x[1]) # this one is constant
    return dx
end

with a taylor series in time whose order is larger than two (the exact solution has degree two), then TaylorIntegration.stepsize! returns a stepsize which is Inf.

This is maybe correct? Since one gets the exact solution for all t? So i don't know the expected action point here, but then some error happens down the line, in TaylorModels.validated_integ.

After #114 is merged, we will need to review the keyword warning list. Related to #114 (comment)

(This is a copy of the list that appeared originally in #3)
To do's before creating the v0.0 tag commit. Main changes are:

• Choose a suitable license, and add it into LICENSE.md #3
• Add .travis.yml #3
• Update jupyter notebooks found in examples/
• Add details to README.md; in particular links to concrete examples. To be discussed: include in README.md a fast, easy example for users who need to use the package right away. #8
• [ ] Documentation

Hi!

I tried to calculate the Lyapunov exponents with your code for a system with functions like tanh and found the following error. Maybe you have a suggestion.

Consider the next slightly modified Lorenz function example from the doc:

function lorenz!(dq, q, params, t)
    σ, ρ, β = params
    x, y, z = q
    dq[1] = σ*tanh(y-x)
    dq[2] = x*(ρ-z)-y
    dq[3] = x*y-β*z
    nothing
end

Then, the line

tv, xv, λv = lyap_taylorinteg(lorenz!, x0, t0, tmax, 28, 1e-20, params; maxsteps=2000000)

leads to the Stacktrace:

/ not defined for HomogeneousPolynomial{Taylor1{Float64}}

Stacktrace:
  [1] error(::String, ::String, ::Type)
    @ Base .\error.jl:42
  [2] no_op_err(name::String, T::Type)
    @ Base .\promotion.jl:395
  [3] /(x::HomogeneousPolynomial{Taylor1{Float64}}, y::HomogeneousPolynomial{Taylor1{Float64}})
    @ Base .\promotion.jl:399
  [4] /(x::HomogeneousPolynomial{Taylor1{Float64}}, y::Int64)
    @ Base .\promotion.jl:324
  [5] tanh!
    @ ~\.julia\packages\TaylorSeries\tGlgp\src\functions.jl:464 [inlined]
  [6] tanh(a::TaylorN{Taylor1{Float64}})
    @ TaylorSeries ~\.julia\packages\TaylorSeries\tGlgp\src\functions.jl:139
  [7] lorenz!(dq::Vector{TaylorN{Taylor1{Float64}}}, q::Vector{TaylorN{Taylor1{Float64}}}, params::Vector{Float64}, t::Taylor1{Float64})
    @ Main .\In[53]:5
  [8] stabilitymatrix!(eqsdiff!::typeof(lorenz!), t::Taylor1{Float64}, x::Vector{Taylor1{Float64}}, δx::Vector{TaylorN{Taylor1{Float64}}}, dδx::Vector{TaylorN{Taylor1{Float64}}}, jac::Matrix{Taylor1{Float64}}, _δv::Vector{TaylorN{Taylor1{Float64}}}, params::Vector{Float64}, jacobianfunc!::Nothing)
    @ TaylorIntegration ~\.julia\packages\TaylorIntegration\hUDy9\src\lyapunovspectrum.jl:26
  [9] lyap_taylorstep!(f!::Function, t::Taylor1{Float64}, x::Vector{Taylor1{Float64}}, dx::Vector{Taylor1{Float64}}, xaux::Vector{Taylor1{Float64}}, δx::Vector{TaylorN{Taylor1{Float64}}}, dδx::Vector{TaylorN{Taylor1{Float64}}}, jac::Matrix{Taylor1{Float64}}, abstol::Float64, _δv::Vector{TaylorN{Taylor1{Float64}}}, varsaux::Array{Taylor1{Float64}, 3}, params::Vector{Float64}, parse_eqs::Bool, jacobianfunc!::Nothing)
    @ TaylorIntegration ~\.julia\packages\TaylorIntegration\hUDy9\src\lyapunovspectrum.jl:184
 [10] lyap_taylorinteg(f!::Function, q0::Vector{Float64}, t0::Float64, tmax::Float64, order::Int64, abstol::Float64, params::Vector{Float64}, jacobianfunc!::Nothing; maxsteps::Int64, parse_eqs::Bool)
    @ TaylorIntegration ~\.julia\packages\TaylorIntegration\hUDy9\src\lyapunovspectrum.jl:267
 [11] top-level scope
    @ .\timing.jl:210 [inlined]
 [12] top-level scope
    @ .\In[56]:0
 [13] eval
    @ .\boot.jl:360 [inlined]
 [14] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
    @ Base .\loading.jl:1094

Thanks for the code, and congratulations on your job!

The current implementation is incompatible with DynamicalODEProblem specifications from DiffEq since

is specific to ODEFunctions and will fail for DynamicalODEFunctions (which have f1 and f2). I tried to address this in #108

I tested a numerical integration of the Kepler problem using TaylorIntegration.jl in the REPL with julia 0.4.7-pre+3, using the latest commit versions of branches masterand lb/changes. The initial conditions are the same as in the Kepler problem example included in TaylorSeries.jl. I ran the following code:

using the master branch I ran:

 julia> using TaylorIntegration, FastAnonymous
const μ = 1.0 #mass parameter
const q0 = [0.19999999999999996, 0.0, 0.0, 3.0] #x, y, vx, vy, initial conditions
const order = 28 #order of Taylor expansion
const t0 = 0.0 #initial time
const tmax = 10000(2π) #final time
const abstol = 1.0E-20 #absolute error tolerance

#the problem at hand
function kepler_problem(x)
    r32 = (x[1]^2+x[2]^2)^(3/2)
    return [x[3], x[4], -μ*x[1]/r32, -μ*x[2]/r32]
end

kepler_anon = @anon x_ -> kepler_problem(x_) #`FastAnonymous` version of `kepler problem`
stepsizeall_anon = @anon (q_, epsilon_) -> stepsizeall(q_, epsilon_) #`FastAnonymous` version of `stepsizeall`

and then...

julia> begin
           q = Array{Float64,1}[] #a vector of vectors which stores independent variable and state variables
           #fill `q` with Float64[]'s named t, x, y, u, v
           for var = (:t, :x, :y, :u, :v)
               @eval $var = Float64[]
               @eval push!($q, $var)
           end
           #Taylor integration of Kepler problem:
           @time integrate!(kepler_anon, stepsizeall_anon, q0, abstol, order, t0, tmax, q)
       end

...three times. The output for these three runs was:

40.873064 seconds (475.63 M allocations: 37.975 GB, 9.64% gc time)
41.547965 seconds (475.41 M allocations: 37.965 GB, 9.36% gc time)
42.512222 seconds (475.41 M allocations: 37.965 GB, 9.31% gc time)

And then, using the lb/changes branch I ran:

 julia> using TaylorIntegration, FastAnonymous
const μ = 1.0 #mass parameter
const q0 = [0.19999999999999996, 0.0, 0.0, 3.0] #x, y, vx, vy, initial conditions
const order = 28 #order of Taylor expansion
const t0 = 0.0 #initial time
const tmax = 10000(2π) #final time
const abstol = 1.0E-20 #absolute error tolerance
const iterations = 500000 #the max number of iterations

#the problem at hand
function kepler_problem(t, x)
    r32 = (x[1]^2+x[2]^2)^(3/2)
    return [x[3], x[4], -μ*x[1]/r32, -μ*x[2]/r32]
end

kepler_anon = @anon (t_, x_) -> kepler_problem(t_, x_) #`FastAnonymous` version of `kepler problem`

and then...

#Taylor integration of Kepler problem:
@time taylorinteg(kepler_anon, q0, t0, tmax, order, abstol, maxsteps=iterations);

...three times.

The output for these three runs was (again using julia 0.4.7-pre+3):

45.157840 seconds (462.13 M allocations: 36.790 GB, 19.75% gc time)
45.280109 seconds (461.89 M allocations: 36.779 GB, 19.93% gc time)
45.644434 seconds (461.89 M allocations: 36.779 GB, 20.49% gc time)

So the new interface is a hell of a lot more user-friendly; but the master version seems to be a bit faster than the lb/changes version (roughly 10%), at least for this example, even as the usage of memory improves! @lbenet what are your thoughts on this?

Edit: for lb/changes I also tried:

julia> @time t, q = taylorinteg(kepler_anon, q0, t0, tmax, order, abstol, maxsteps=iterations);

Again, I did three runs. The results were:

45.110408 seconds (462.13 M allocations: 36.790 GB, 19.24% gc time)
45.034779 seconds (461.89 M allocations: 36.779 GB, 19.72% gc time)
45.497351 seconds (461.89 M allocations: 36.779 GB, 20.16% gc time)

Perhaps the issue has to do with garbage collection time?

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

I have the following:

using TaylorSeries, TaylorIntegration
#eqs of motion
saddle(t,x) = [-x[1],x[2]]

#polynomial independent variables
δξ = set_variables("ξ", numvars=2, order=polyg_order)
x0 = [0.,0.0]
x0TN = x0 + δξ

const t0 = 0.0
const tmax = 4.
const Δt = 0.25
const order = 15
const tol = 1.0e-15

tv = t0:Δt:tmax
xv = taylorinteg(saddle,x0TN,tv,order,tol)

This code worked before, but now (I'm using julia0.5.3) it outputs the following error

MethodError: no method matching saddle(::TaylorSeries.Taylor1{Float64}, ::Array{TaylorSeries.Taylor1{TaylorSeries.TaylorN{Float64}},1}, ::Array{TaylorSeries.Taylor1{TaylorSeries.TaylorN{Float64}},1})
Closest candidates are:
  saddle(::Any, ::Any) at In[6]:2

 in jetcoeffs!(::#saddle, ::Float64, ::Array{TaylorSeries.Taylor1{TaylorSeries.TaylorN{Float64}},1}, ::Array{TaylorSeries.Taylor1{TaylorSeries.TaylorN{Float64}},1}, ::Array{TaylorSeries.Taylor1{TaylorSeries.TaylorN{Float64}},1}, ::Array{Float64,1}) at /Users/blaskolic/.julia/v0.5/TaylorIntegration/src/jettransport.jl:41
 in taylorstep!(::Function, ::Array{TaylorSeries.Taylor1{TaylorSeries.TaylorN{Float64}},1}, ::Array{TaylorSeries.Taylor1{TaylorSeries.TaylorN{Float64}},1}, ::Array{TaylorSeries.Taylor1{TaylorSeries.TaylorN{Float64}},1}, ::Float64, ::Float64, ::Array{TaylorSeries.TaylorN{Float64},1}, ::Int64, ::Float64, ::Array{Float64,1}) at /Users/blaskolic/.julia/v0.5/TaylorIntegration/src/jettransport.jl:97
 in #taylorinteg#14(::Int64, ::Function, ::Function, ::Array{TaylorSeries.TaylorN{Float64},1}, ::FloatRange{Float64}, ::Int64, ::Float64) at /Users/blaskolic/.julia/v0.5/TaylorIntegration/src/jettransport.jl:263
 in taylorinteg(::Function, ::Array{TaylorSeries.TaylorN{Float64},1}, ::FloatRange{Float64}, ::Int64, ::Float64) at /Users/blaskolic/.julia/v0.5/TaylorIntegration/src/jettransport.jl:234

I also tried to put t0 and tmax as separate arguments but it doesn't work either. Is there any substantial change on the internal function methods that doesn't make this work?

Thanks in advance!

Hey,
This package has really evolved and looks really cool! I was thinking it would be nice to add common interface bindings. To do this, it would just need to add a dependency on DiffEqBase and implement a solve function which translates how to solve an ODEProblem to the solver for this package.

Taking a quick look at the examples, it looks like it just would extract the information of the ODEProblem (the function, timespan, and initial condition), and use it to build the

taylorinteg(kepler_problem, q0, t0, 0.01, order, abs_tol, maxsteps)

command. The algorithm could be

alg = TaylorIntegration(order=SomeDefaultHere)

And then

solve(prob,alg,maxiters=...,abstol=...)

would be a solve command that translates the arguments and calls taylorinteg appropriately. Would you accept such a PR?

Now, the big issue that I see from the examples is that they're not using a inplace format. Your examples seem to use the ODE.jl du=f(t,u) version, but doesn't include support for f(t,u,du). Is that correct? If so, that hampers performance quite a bit, but we can translate any inplace function to the not-inplace form using the same translation as the ODE.jl common interface bindings.

@JuliaRegistrator register

I am trying to solve a nonlinear gasdynamics ODE system with high accuracy. The ODE system does have a 0/0 pole it approaches, and the goal is to be able to use TaylorIntegration.jl to approach that pole accurately and robustly. However, I began by checking the TaylorIntegration.jl solution against Vern7 from OrdinaryDiffEq.jl with abstol=1e-12, reltol=0.0. The Vern7 solver produces good results on the ODE system (as verified through BigFloat solutions), so I began by comparing a Vern7 and TaylorMethod solution away from the pole. However, even with similar or tighter tolerances on the TaylorMethod solution, the Euclidean distance between solution points grows to Order(1e-3), so ~9 orders of magnitude higher than the integration tolerance.

I also verified that I am able to get very comparable high-accuracy results between TaylorIntegration.jl and OrdinaryDiffEq.jl Vern7 solutions for a nonlinear test ODE system where no subtractive cancellation takes place in the denominators, so I don't believe this to be user error (but I could be wrong/missing something!).

Attached is a Pluto.jl notebook, but I will also paste code below. Note that I tried different values of the method order, and leaving on vs. omitting parse_eqs and they didn't seem to ultimately change the magnitude of solution differences. Again, the tspan right-side endpoint prescribed below is still some distance away from the singularity, and the plot shows significant differences building up quite early in the solution trajectory.

using OrdinaryDiffEq
using TaylorIntegration
using LinearAlgebra
using Plots

function TMM_odesys!(primes, dependent, params, independent)
    γ = params[1]
    θ = independent
    Mr = dependent[1]
    Mθ = dependent[2]
    gm1o2 = (γ - 1) / 2
    Mθsq = Mθ * Mθ
    Mθsqm1 = Mθsq - 1
    Mθcotθ = Mθ * cot(θ)
    # slightly excessive paretheses, but equations are correct
    primes[1] = ((Mθ * ((gm1o2 * ((Mr * Mr) + (Mr * Mθcotθ))) + Mθsq - 1)) / Mθsqm1)
    primes[2] = (((((gm1o2 * Mθsq) + 1) * (Mr + Mθcotθ)) - (Mr * Mθsqm1)) / Mθsqm1)
    return nothing
end

# Simple Euclidean distance between two points
function solptdiff(u0, u1)
  diff = u0 - u1
  sqrt(dot(diff, diff))
end

cf_test_prob = ODEProblem(TMM_odesys!, [2.5, 0.0], (deg2rad(0.1), 0.4), [1.4])

Vern7sol = solve(cf_test_prob, Vern7(), abstol = 1e-12, reltol = 0.0)

TIsol = solve(cf_test_prob, TaylorMethod(30), abstol = 1e-13, parse_eqs = false)

solptdiff(Vern7sol.u[end], TIsol.u[end]) # On my system, this returns 0.004776986921762344 , much higher than expected

# Next use taylorinteg() to match time points with the Vern7 solution so differences can be compared step-by-step
matching_t_sol = taylorinteg(TMM_odesys!, [2.5, 0.0], Vern7sol.t, 30, 1e-13, [1.4], parse_eqs = false)

# Convenience function to change the returned matrix to a vector of vectors so the difference function can be simply broadcast
function mat2vecvec(m)
  nr = size(m, 1)
  v = Vector{Vector{typeof(m[1,1])}}(undef, nr)
  for i = 1:nr
    v[i] = m[i, :]
  end
  return v
end

matching_t_sol_vv = mat2vecvec(matching_t_sol);

diffs = solptdiff.(Vern7sol.u, matching_t_sol_vv);

# Omit the first point since they are both the IC and the difference is a zero vector which doesn't work with the semilog plotting
plot(Vern7sol.t[2:end], diffs[2:end], yscale = :log10, legend_position = false)

taylorinteg_cf_test.jl.txt

Important changes are coming or in-process in JuliaDiffEq, as it has been recently announced). So I am opening this issue to discuss how to approach that, in view of the possibility of using TaylorIntegration from JuliaDiffEqs (#44) .

Possibilities are:

1. adapt the convention used for the definition of the equations of motion; currently we use the forms f(t,u) or f!(t,u,du);
2. pin the version of DiffEqs used (to v3.1.0);
3. adapt the call from JuliaDiffEqs, such that f(u,p,t) --> f(t,u) or f(du,u,p,t) --> f!(t, u, du);
4. something else.

From those options, I prefer the 3rd option though the 1st one could be a longer-term solution. (I simply like the convention independent and dependent/mutating variables in that order, but this is not the point.) In any case, I think we will need to use some sort of anonymous functions trick to get around the parameters, e.g., by using (du,u,t) -> f(du,u,p,t) in the function field.

@ChrisRackauckas Would that suffice?

In many applications, people is interested in the maximum Lyapunov exponent only. We should have a function which provides it.

   _       _ _(_)_     |  A fresh approach to technical computing
  (_)     | (_) (_)    |  Documentation: http://docs.julialang.org
   _ _   _| |_  __ _   |  Type "?help" for help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 0.5.1-pre+2 (2016-09-20 03:34 UTC)
 _/ |\__'_|_|_|\__'_|  |  Commit f0d40ec* (21 days old release-0.5)
|__/                   |  x86_64-apple-darwin14.5.0

julia> versioninfo()
Julia Version 0.5.1-pre+2
Commit f0d40ec* (2016-09-20 03:34 UTC)
Platform Info:
  System: Darwin (x86_64-apple-darwin14.5.0)
  CPU: Intel(R) Core(TM) i7-4750HQ CPU @ 2.00GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell)
  LAPACK: libopenblas64_
  LIBM: libopenlibm
  LLVM: libLLVM-3.7.1 (ORCJIT, haswell)

julia> using TaylorIntegration

julia> f(t, x) = x.^2
f (generic function with 1 method)

julia> @time tT, xT = taylorinteg(f, 3, 0, 0.3, 25, 1e-20, maxsteps=50);
ERROR: MethodError: no method matching taylorinteg(::#f, ::Int64, ::Int64, ::Float64, ::Int64, ::Float64; maxsteps=50)
Closest candidates are:
  taylorinteg{T<:Number}(::Any, ::T<:Number, ::T<:Number, ::T<:Number, ::Int64, ::T<:Number; maxsteps) at /Users/Jorge/.julia/v0.5/TaylorIntegration/src/integration_methods.jl:231
  taylorinteg{T<:Number}(::Any, ::Array{T<:Number,1}, ::T<:Number, ::T<:Number, ::Int64, ::T<:Number; maxsteps) at /Users/Jorge/.julia/v0.5/TaylorIntegration/src/integration_methods.jl:272
  taylorinteg{T<:Number}(::Any, ::T<:Number, ::Range{T<:Number}, ::Int64, ::T<:Number; maxsteps) at /Users/Jorge/.julia/v0.5/TaylorIntegration/src/integration_methods.jl:316
  ...

julia> 

Note that the initial condition x0 and initial time t0 are Int64's, whereas the final time is a Float64..

After merging #31, the documentation was not properly updated (deployed): everything related with @taylorize does not appear in the docs. I am opening this issue as a reminder to check if the documentation is properly built once #62 is merged. If not, we should probably open an issue in Documenter.jl.

I'm working in a 3-body problem and I got to the following situation where for some initial conditions (In particular [0,0]) an Assertion Error is thrown.

using TaylorSeries, TaylorIntegration

R_1 = -0.37037037037037035
R_2 = 9.62962962962963
Ω  = 0.16431676725154984
M_1 = 26.0
M_2 = 1.0

function test!(t,x,dx)
   
    r_13_cubed = ( (x[1]-R_1)^2 + x[2]^2 )^(3/2)
    r_23_cubed = ( (x[1]-R_2)^2 + x[2]^2 )^(3/2)
    
    dx[1] = -( M_1*(x[1]-R_1)/r_13_cubed + M_2*(x[1]-R_2)/r_23_cubed ) + Ω*(Ω*x[1])
    dx[2] = (-M_1/r_13_cubed -M_2/r_23_cubed + Ω^2)*x[2] 
    
    nothing
end

ord_taylor = 28
abstol = 1.0e-20

q0 = [0.0,0.0]
t0 = 0.0
tmax = 10
tv = linspace(t0,tmax,tmax*10)

taylorinteg(test!, q0, tv, ord_taylor, abstol);

AssertionError: t1 > t0

I tried simplifying the equations of motion but it seems that this error is sensitive to initial conditions and constants parameters also. Note that q0 = [0,0] doesn't blow up or anything...

Any ideas of what could be happening?

@JuliaRegistrator register

Opening this issue to remember that the gifs in the docs for Poincaré maps are not displayed.

Currently jetcoeffs!(eqsdiff!::Function, t, x, dx, xaux) expects eqsdiff! to follow with taylorinteg so f!(t, x, dx). However, it would be much cleaner if this followed the updated DifferentialEquations.jl convention, where in-place dynamics equations are given as f!(dx,x,p,t) where p are parameters for f!.

e.g. https://docs.juliadiffeq.org/stable/tutorials/ode_example.html#Example-2:-Solving-Systems-of-Equations-1

@JuliaRegistrator register

I'm trying to model a differential equation (in polar coordenates) that has a periodic orbit in r=a

using TaylorSeries, TaylorIntegration
periodic_orbit(t,x) = [x[1](x[1]-a),b]

with a and b constants.

Integrating this using

taylorinteg(periodic_orbit,x0,tv,15,1.0e-15)

with

x0 = [1.0 ξ₁ + 𝒪(‖x‖¹⁶) , 0.005 + 1.0 ξ₂ + 𝒪(‖x‖¹⁶)]
tv = t0:Δt:t

outputs the following error

UndefRefError: access to undefined reference

 in jetcoeffs!(::#per_orb_bad, ::Float64, ::Array{TaylorSeries.Taylor1{TaylorSeries.TaylorN{Float64}},1}) at /Users/<something>/.julia/v0.5/TaylorIntegration/src/jettransport.jl:54
 in taylorstep!(::Function, ::Float64, ::Array{TaylorSeries.TaylorN{Float64},1}, ::Int64, ::Float64) at /Users/<something>/.julia/v0.5/TaylorIntegration/src/jettransport.jl:147
 in #taylorinteg#10(::Int64, ::Function, ::Function, ::Array{TaylorSeries.TaylorN{Float64},1}, ::FloatRange{Float64}, ::Int64, ::Float64) at /Users/<something>/.julia/v0.5/TaylorIntegration/src/jettransport.jl:412
 in taylorinteg(::Function, ::Array{TaylorSeries.TaylorN{Float64},1}, ::FloatRange{Float64}, ::Int64, ::Float64) at /Users/<something>/.julia/v0.5/TaylorIntegration/src/jettransport.jl:392

if instead

periodic_orbit2(t,x) = [x[1](x[1]-a),b*x[2]]

is used, the integrator runs as it should.

I believe this has something to do with not making any operations for the second element of the array of the periodic_orbit function, but honestly I have no idea why his wouldn't work.

Thank you for this great package. I was wondering if it possible to extend the functionality to systems with control input.

Is it possible to use @taylorize with

x_n=f(x,u,t)

Currently, we have different ways to set the initial and final times for the integration, which defines which data is returned. DifferentialEquations distinguishes tspan (a tuple for the initial and final time) from saveat, which turns to be pretty handy.

In current master, when trying to use @taylorize with the following function definition

@taylorize function f!(dq, q, params, t)
    for i in 1:10
        if i == 1
            dq[i] = 2q[i]
        elseif i == 2
            dq[i] = -q[i]
        else
            dq[i] = q[i]
        end
    end
    nothing
end

An error occurs:

ERROR: LoadError: ArgumentError: collection must be non-empty

If, instead, the elseif statement is removed, things work fine:

@taylorize function f!(dq, q, params, t)
    for i in 1:10
        if i == 1
            dq[i] = 2q[i]
        else
            dq[i] = q[i]
        end
    end
    nothing
end

The problem seems to have been discussed already here.

Hello, very interesting package!

I just looked through the docs and found that you are comparing taylorintegration to Vern9. I was surprised to see that this higher order explicit method fails to visit the second mass and noticed that you are calling the solver with

solV = solve(prob, Vern9(), abstol=1e-20); #solve probwith theVern9 method

that is, the reltol value is not set, in which case it should default to reltol=1e-3 (see here), at least this is what happens for lower order methods. Do you get the same qualitative behaviour for

solV = solve(prob, Vern9(), reltol = 1e-20, abstol=1e-20);?

Greetings.

When using lyap_taylorinteg to work out the Lyapunov spectrum of the Lorenz system, an error was throw out.
(with Julia 1.6.1)

julia> using TaylorIntegration, TaylorSeries

julia> function lorenz!(t, x, dx)
dx[1] = σ*(x[2]-x[1])
dx[2] = x[1]*(ρ-x[3])-x[2]
dx[3] = x[1]x[2]-βx[3]
nothing
end
lorenz! (generic function with 1 method)

julia> const σ = 10.71
10.71

julia> const ρ = 119
119

julia> const β = 8.93
8.93

julia>

julia>

julia> const x0 = [19.0, 20.0, 50.0]
3-element Vector{Float64}:
19.0
20.0
50.0

julia> const t0 = 0.0
0.0

julia> const tmax = 100.0
100.0

julia> @time tv, xv, λv = lyap_taylorinteg(lorenz!, x0, t0, tmax, 28, 1e-20; maxsteps=2000000);
ERROR: AssertionError: length(q0) must be equal to number of variables set by TaylorN
Stacktrace:
[1] lyap_taylorinteg(f!::Function, q0::Vector{Float64}, t0::Float64, tmax::Float64, order::Int64, abstol::Float64, params::Nothing, jacobianfunc!::Nothing; maxsteps::Int64, parse_eqs::Bool)
@ TaylorIntegration C:\Users\gsjiang.julia\packages\TaylorIntegration\hUDy9\src\lyapunovspectrum.jl:241
[2] top-level scope
@ .\timing.jl:210 [inlined]
[3] top-level scope
@ .\REPL[9]:0

Why this error occurs?

Right now, it's not possible to integrate ODEs backwards. I know it's that way by design, but why is it that way?
Sometimes it's practical to integrate "to the past" and I think it would be a nice feature to have.
I found this by trying to do a Newton's method around an specific time.

In https://github.com/JuliaReach/Reachability.jl/issues/146 it was reported a problem when using @taylorize with one equation equal to zero. Opening this to keep track of the problem.
This may be related to #20.