Introducing TaylorSolution

Project.toml

@@ -1,7 +1,7 @@
 name = "TaylorIntegration"
 uuid = "92b13dbe-c966-51a2-8445-caca9f8a7d42"
 repo = "https://github.com/PerezHz/TaylorIntegration.jl.git"
-version = "0.15.3"
+version = "0.16.0"
 
 [deps]
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
@@ -19,7 +19,7 @@ TaylorSeries = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 
 [extensions]
-TaylorIntegrationDiffEq = "OrdinaryDiffEq"
+TaylorIntegrationDiffEqExt = "OrdinaryDiffEq"
 
 [compat]
 Aqua = "0.8"

docs/src/api.md

@@ -14,6 +14,12 @@ lyap_taylorinteg
 @taylorize
 ```
 
+## Exported types
+
+```@docs
+TaylorSolution
+```
+
 ## Internal
 
 ```@autodocs

docs/src/kepler.md

@@ -74,7 +74,8 @@ We now perform the integration, using a 25 order expansion and
 absolute tolerance of ``10^{-20}``.
 ```@example kepler
 using TaylorIntegration, Plots
-t, q = taylorinteg(kepler_eqs!, q0, 0.0, 10000*2pi, 25, 1.0e-20, maxsteps=700_000);
+sol = taylorinteg(kepler_eqs!, q0, 0.0, 10000*2pi, 25, 1.0e-20, maxsteps=700_000);
+t, q = sol.t, sol.x;
 t[end], q[end,:]
 ```
 

docs/src/lorenz_lyapunov.md

@@ -114,12 +114,14 @@ one with the values of the Lyapunov spectrum.
 
 We first carry out the integration computing internally the Jacobian
 ```@example lorenz
-tv, xv, λv = lyap_taylorinteg(lorenz!, x0, t0, tmax, 28, 1e-20, params; maxsteps=2000000);
+sol = lyap_taylorinteg(lorenz!, x0, t0, tmax, 28, 1e-20, params; maxsteps=2000000);
+tv, xv, λv = sol.t, sol.x, sol.λ
 nothing # hide
 ```
 Now, the integration is obtained exploiting `lorenz_jac!`:
 ```@example lorenz
-tv_, xv_, λv_ = lyap_taylorinteg(lorenz!, x0, t0, tmax, 28, 1e-20, params, lorenz_jac!; maxsteps=2000000);
+sol_ = lyap_taylorinteg(lorenz!, x0, t0, tmax, 28, 1e-20, params, lorenz_jac!; maxsteps=2000000);
+tv_, xv_, λv_ = sol_.t, sol_.x, sol_.λ
 nothing # hide
 ```
 In terms of performance the second method is about ~50% faster than the first.

docs/src/pendulum.md

@@ -103,7 +103,7 @@ We perform the Taylor integration using the initial condition `x0TN`, during
 6 periods of the pendulum (i.e., ``6T``), exploiting multiple dispatch:
 ```@example pendulum
 tv = t0:integstep:tmax # the times at which the solution will be evaluated
-xv = taylorinteg(pendulum!, q0TN, tv, order, abstol)
+sol = taylorinteg(pendulum!, q0TN, tv, order, abstol)
 nothing # hide
 ```
 
@@ -128,19 +128,19 @@ nothing # hide
 ```
 
 We evaluate the jet at ``\partial U_{x(t)}`` (the boundary of ``U_{x(t)}``) at each
-value of the solution vector `xv`; we organize these values such that we can plot
+value of the solution array `sol.x`; we organize these values such that we can plot
 them later:
 ```@example pendulum
-xjet_plot = map(λ->λ.(ξv), xv[:,1])
-pjet_plot = map(λ->λ.(ξv), xv[:,2])
+xjet_plot = map(λ->λ.(ξv), sol.x[:,1])
+pjet_plot = map(λ->λ.(ξv), sol.x[:,2])
 nothing # hide
 ```
 Above, we have exploited the fact that `Array{TaylorN{Float64}}` variables are
 callable objects. Now, we evaluate the jet at the nominal solution, which
-corresponds to ``\xi=(0,0)``, at each value of the solution vector `xv`:
+corresponds to ``\xi=(0,0)``, at each value of the solution vector `sol.x`:
 ```@example pendulum
-x_nom = xv[:,1]()
-p_nom = xv[:,2]()
+x_nom = sol.x[:,1]()
+p_nom = sol.x[:,2]()
 nothing # hide
 ```
 

docs/src/root_finding.md

@@ -124,11 +124,11 @@ for i in 1:nconds
     x_ini .= x0 .+ 0.005 .* [0.0, sqrt(rand1)*cos(2pi*rand2), 0.0, sqrt(rand1)*sin(2pi*rand2)]
     px!(x_ini, E0)   # ensure initial energy is E0
 
-    tv_i, xv_i, tvS_i, xvS_i, gvS_i = taylorinteg(henonheiles!, g, x_ini, 0.0, 135.0,
+    sol_i = taylorinteg(henonheiles!, g, x_ini, 0.0, 135.0,
         25, 1e-25, maxsteps=30000);
-    tvSv[i] = vcat(0.0, tvS_i)
-    xvSv[i] = vcat(transpose(x_ini), xvS_i)
-    gvSv[i] = vcat(0.0, gvS_i)
+    tvSv[i] = vcat(0.0, sol_i.t)
+    xvSv[i] = vcat(transpose(x_ini), sol_i.x)
+    gvSv[i] = vcat(0.0, sol_i.gresids)
 end
 nothing # hide
 ```
@@ -202,14 +202,15 @@ nothing # hide
 
 We are now set to carry out the integration.
 ```@example poincare
-tvTN, xvTN, tvSTN, xvSTN, gvSTN = taylorinteg(henonheiles!, g, x0TN, 0.0, 135.0, 25, 1e-25, maxsteps=30000);
+solTN = taylorinteg(henonheiles!, g, x0TN, 0.0, 135.0, 25, 1e-25, maxsteps=30000);
 nothing # hide
 ```
 
 We define some auxiliary arrays, and then make an animation with the results for
 plotting.
 ```@example poincare
 #some auxiliaries:
+tvTN, xvTN, tvSTN, xvSTN, gvSTN = solTN.t, solTN.x, solTN.tevents, solTN.xevents, solTN.gresids
 xvSTNaa = Array{Array{TaylorN{Float64},1}}(undef, length(tvSTN)+1 );
 xvSTNaa[1] = x0TN
 for ind in 2:length(tvSTN)+1

docs/src/simple_example.md

@@ -95,24 +95,24 @@ below we shall *try* to compute it up to ``t_\textrm{end}=0.34``; as we shall
 see, Taylor's method takes care of this. For
 the integration presented below, we use a 25-th series expansion, with
 ``\epsilon_\textrm{tol} = 10^{-20}``, and compute up to 150
-integration steps.
+integration steps. The type of the solution returned by `taylorinteg` is [`TaylorSolution`](@ref).
 
 ```@example example1
-tT, xT = taylorinteg(diffeq, 3.0, 0.0, 0.34, 25, 1e-20, maxsteps=150) ;
+sol = taylorinteg(diffeq, 3.0, 0.0, 0.34, 25, 1e-20, maxsteps=150);
 ```
 
 We first note that the last point of the
 calculation does not exceed ``t_\textrm{max}``.
 ```@example example1
-tT[end]
+sol.t[end]
 ```
-Increasing the `maxsteps` parameter pushes `tT[end]` closer to ``t_\textrm{max}``
+Increasing the `maxsteps` parameter pushes `sol.t[end]` closer to ``t_\textrm{max}``
 but it actually does not reach this value.
 
 Figure 1 displays the computed solution as a function of
 time, in log scale.
 ```@example example1
-plot(tT, log10.(xT), shape=:circle)
+plot(sol.t, log10.(sol.x), shape=:circle)
 xlabel!("t")
 ylabel!("log10(x(t))")
 xlims!(0,0.34)
@@ -128,8 +128,8 @@ and the analytical solution in terms of time.
 
 ```@example example1
 exactsol(t, x0) = x0 / (1 - x0 * t)
-δxT = abs.(xT .- exactsol.(tT, 3.0)) ./ exactsol.(tT, 3.0);
-plot(tT[6:end], log10.(δxT[6:end]), shape=:circle)
+δxT = abs.(sol.x .- exactsol.(sol.t, 3.0)) ./ exactsol.(sol.t, 3.0);
+plot(sol.t[6:end], log10.(sol.x[6:end]), shape=:circle)
 xlabel!("t")
 ylabel!("log10(dx(t))")
 xlims!(0, 0.4)
@@ -141,8 +141,8 @@ impose (arbitrarily) a relative accuracy of ``10^{-13}``; the time until
 such accuracy is satisfied is given by:
 ```@example example1
 indx = findfirst(δxT .> 1.0e-13);
-esol = exactsol(tT[indx-1],3.0);
-tT[indx-1], esol, eps(esol)
+esol = exactsol(sol.t[indx-1],3.0);
+sol.t[indx-1], esol, eps(esol)
 ```
 Note that, the accuracy imposed in terms of the actual value
 of the exact solution means that the difference of the computed

docs/src/taylorize.md

@@ -56,7 +56,7 @@ tf = 10000.0
 q0 = [1.3, 0.0]
 
 # The actual integration
-t1, x1 = taylorinteg(pendulumNP!, q0, t0, tf, 25, 1e-20, maxsteps=50000); # warm-up run
+sol1 = taylorinteg(pendulumNP!, q0, t0, tf, 25, 1e-20, maxsteps=50000); # warm-up run
 e1 = @elapsed taylorinteg(pendulumNP!, q0, t0, tf, 25, 1e-20, maxsteps=50000);
 all1 = @allocated taylorinteg(pendulumNP!, q0, t0, tf, 25, 1e-20, maxsteps=50000);
 e1, all1
@@ -119,7 +119,7 @@ use the same instruction as above; the default value for the keyword argument `p
 is `true`, so we may omit it.
 
 ```@example taylorize
-t2, x2 = taylorinteg(pendulum!, q0, t0, tf, 25, 1e-20, maxsteps=50000); # warm-up run
+sol2 = taylorinteg(pendulum!, q0, t0, tf, 25, 1e-20, maxsteps=50000); # warm-up run
 e2 = @elapsed taylorinteg(pendulum!, q0, t0, tf, 25, 1e-20, maxsteps=50000);
 all2 = @allocated taylorinteg(pendulum!, q0, t0, tf, 25, 1e-20, maxsteps=50000);
 e2, all2
@@ -132,7 +132,7 @@ e1/e2, all1/all2
 
 We can check that both integrations yield the same results.
 ```@example taylorize
-t1 == t2 && x1 == x2
+sol1.t == sol2.t && sol1.x == sol2.x
 ```
 
 As stated above, in order to allow opting out of using the specialized method