diff --git a/Project.toml b/Project.toml index dd0d5da43..e4f276427 100644 --- a/Project.toml +++ b/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" diff --git a/docs/src/api.md b/docs/src/api.md index 8ecfe9145..0f47f88de 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -14,6 +14,12 @@ lyap_taylorinteg @taylorize ``` +## Exported types + +```@docs +TaylorSolution +``` + ## Internal ```@autodocs diff --git a/docs/src/kepler.md b/docs/src/kepler.md index 5675cda47..5a3307d21 100644 --- a/docs/src/kepler.md +++ b/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,:] ``` diff --git a/docs/src/lorenz_lyapunov.md b/docs/src/lorenz_lyapunov.md index 021f53373..f6af879e0 100644 --- a/docs/src/lorenz_lyapunov.md +++ b/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. diff --git a/docs/src/pendulum.md b/docs/src/pendulum.md index c3c8adb20..f445a4321 100644 --- a/docs/src/pendulum.md +++ b/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 ``` diff --git a/docs/src/root_finding.md b/docs/src/root_finding.md index c4bbe616c..14f4be621 100644 --- a/docs/src/root_finding.md +++ b/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,7 +202,7 @@ 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 ``` @@ -210,6 +210,7 @@ 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 diff --git a/docs/src/simple_example.md b/docs/src/simple_example.md index 6887bcf03..ef441f4ee 100644 --- a/docs/src/simple_example.md +++ b/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 diff --git a/docs/src/taylorize.md b/docs/src/taylorize.md index 76bdbd698..6d45ec8af 100644 --- a/docs/src/taylorize.md +++ b/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 diff --git a/ext/TaylorIntegrationDiffEq.jl b/ext/TaylorIntegrationDiffEqExt.jl similarity index 68% rename from ext/TaylorIntegrationDiffEq.jl rename to ext/TaylorIntegrationDiffEqExt.jl index f0afddeb3..9de736664 100644 --- a/ext/TaylorIntegrationDiffEq.jl +++ b/ext/TaylorIntegrationDiffEqExt.jl @@ -1,27 +1,25 @@ # This file is part of the TaylorIntegration.jl package; MIT licensed -module TaylorIntegrationDiffEq +module TaylorIntegrationDiffEqExt using TaylorIntegration if isdefined(Base, :get_extension) - using OrdinaryDiffEq - import OrdinaryDiffEq: OrdinaryDiffEqAdaptiveAlgorithm, - OrdinaryDiffEqConstantCache, OrdinaryDiffEqMutableCache, - alg_order, alg_cache, initialize!, perform_step!, @unpack, - @cache, stepsize_controller!, step_accept_controller!, _ode_addsteps! + using OrdinaryDiffEq: @unpack, @cache, OrdinaryDiffEqAdaptiveAlgorithm, + OrdinaryDiffEqConstantCache, OrdinaryDiffEqMutableCache, ODEFunction, + DynamicalODEFunction, check_keywords, warn_compat, ODEProblem, DynamicalODEProblem + import OrdinaryDiffEq else - using ..OrdinaryDiffEq - import ..OrdinaryDiffEq: OrdinaryDiffEqAdaptiveAlgorithm, - OrdinaryDiffEqConstantCache, OrdinaryDiffEqMutableCache, - alg_order, alg_cache, initialize!, perform_step!, @unpack, - @cache, stepsize_controller!, step_accept_controller!, _ode_addsteps! + using ..OrdinaryDiffEq: @unpack, @cache, OrdinaryDiffEqAdaptiveAlgorithm, + OrdinaryDiffEqConstantCache, OrdinaryDiffEqMutableCache, ODEFunction, + DynamicalODEFunction, check_keywords, warn_compat, ODEProblem, DynamicalODEProblem + import ..OrdinaryDiffEq end using StaticArrays: SVector, SizedArray -using RecursiveArrayTools: ArrayPartition +using RecursiveArrayTools: ArrayPartition, copyat_or_push! -import DiffEqBase: ODEProblem, solve, ODE_DEFAULT_NORM +import DiffEqBase # TODO: check which keywords work fine const warnkeywords = (:save_idxs, :d_discontinuities, :unstable_check, :save_everystep, @@ -29,7 +27,7 @@ const warnkeywords = (:save_idxs, :d_discontinuities, :unstable_check, :save_eve :dtmin, :force_dtmin, :internalnorm, :gamma, :beta1, :beta2, :qmax, :qmin, :qsteady_min, :qsteady_max, :qoldinit, :failfactor, :isoutofdomain, :unstable_check, - :calck, :progress, :timeseries_steps, :dense) + :calck, :progress, :timeseries_steps) global warnlist = Set(warnkeywords) @@ -41,16 +39,17 @@ struct TaylorMethodParams <: TaylorAlgorithm parse_eqs::Bool end -import TaylorIntegration: TaylorMethod +using TaylorIntegration: RetAlloc +import TaylorIntegration: TaylorMethod, update_jetcoeffs_cache! TaylorMethod(order; parse_eqs=true) = TaylorMethodParams(order, parse_eqs) # set `parse_eqs` to `true` by default -alg_order(alg::TaylorMethodParams) = alg.order +OrdinaryDiffEq.alg_order(alg::TaylorMethodParams) = alg.order TaylorMethod() = error("Maximum order must be specified for the Taylor method") # overload DiffEqBase.ODE_DEFAULT_NORM for Taylor1 arrays -ODE_DEFAULT_NORM(x::AbstractArray{Taylor1{T}, N},y) where {T<:Number, N} = norm(x, Inf) +DiffEqBase.ODE_DEFAULT_NORM(x::AbstractArray{<:AbstractSeries}, y) = norm(x, Inf) ### cache stuff struct TaylorMethodCache{uType, rateType, tTType, uTType} <: OrdinaryDiffEqMutableCache @@ -64,20 +63,16 @@ struct TaylorMethodCache{uType, rateType, tTType, uTType} <: OrdinaryDiffEqMutab duT::uTType uauxT::uTType parse_eqs::Ref{Bool} - rv::TaylorIntegration.RetAlloc + rv::RetAlloc end -# full_cache(c::TaylorMethodCache) = begin -# tuple(c.u, c.uprev, c.tmp, c.k, c.fsalfirst, c.tT, c.uT, c.duT, c.uauxT, c.parse_eqs, c.rv) -# end - struct TaylorMethodConstantCache{uTType} <: OrdinaryDiffEqConstantCache uT::uTType parse_eqs::Ref{Bool} - rv::TaylorIntegration.RetAlloc{uTType} + rv::RetAlloc{uTType} end -function alg_cache(alg::TaylorMethodParams, u, rate_prototype, uEltypeNoUnits, +function OrdinaryDiffEq.alg_cache(alg::TaylorMethodParams, u, rate_prototype, uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits, uprev, uprev2, f, t, dt, reltol, p, calck,::Val{true}) order = alg.order @@ -104,7 +99,7 @@ end # This method is used for DynamicalODEFunction's (`parse_eqs=false`): tmpT1 and arrT1 # must have the proper type to build `TaylorMethodCache` -function alg_cache(alg::TaylorMethodParams, u::ArrayPartition, rate_prototype, uEltypeNoUnits, +function OrdinaryDiffEq.alg_cache(alg::TaylorMethodParams, u::ArrayPartition, rate_prototype, uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits, uprev, uprev2, f, t, dt, reltol, p, calck, ::Val{true}) order = alg.order @@ -129,7 +124,7 @@ function alg_cache(alg::TaylorMethodParams, u::ArrayPartition, rate_prototype, u ) end -function alg_cache(alg::TaylorMethodParams, u, rate_prototype, uEltypeNoUnits, +function OrdinaryDiffEq.alg_cache(alg::TaylorMethodParams, u, rate_prototype, uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits, uprev, uprev2, f, t, dt, reltol, p, calck, ::Val{false}) order = alg.order @@ -140,7 +135,7 @@ function alg_cache(alg::TaylorMethodParams, u, rate_prototype, uEltypeNoUnits, return TaylorMethodConstantCache(Taylor1(u, alg.order), Ref(parse_eqs), rv) end -function initialize!(integrator, c::TaylorMethodConstantCache) +function OrdinaryDiffEq.initialize!(integrator, c::TaylorMethodConstantCache) @unpack u, t, f, p = integrator tT = Taylor1(typeof(t), integrator.alg.order) tT[0] = t @@ -157,7 +152,7 @@ function initialize!(integrator, c::TaylorMethodConstantCache) integrator.k[2] = integrator.fsallast end -function perform_step!(integrator,cache::TaylorMethodConstantCache) +function OrdinaryDiffEq.perform_step!(integrator,cache::TaylorMethodConstantCache) @unpack u, t, dt, f, p = integrator tT = Taylor1(typeof(t), integrator.alg.order) tT[0] = t+dt @@ -172,7 +167,7 @@ function perform_step!(integrator,cache::TaylorMethodConstantCache) integrator.u = u end -function initialize!(integrator, cache::TaylorMethodCache) +function OrdinaryDiffEq.initialize!(integrator, cache::TaylorMethodCache) @unpack u, t, f, p = integrator @unpack k, fsalfirst, tT, uT, duT, uauxT, parse_eqs, rv = cache TaylorIntegration.__jetcoeffs!(Val(parse_eqs.x), f, tT, uT, duT, uauxT, p, rv) @@ -187,23 +182,23 @@ function initialize!(integrator, cache::TaylorMethodCache) integrator.stats.nf += 1 end -function perform_step!(integrator, cache::TaylorMethodCache) +function OrdinaryDiffEq.perform_step!(integrator, cache::TaylorMethodCache) @unpack t, dt, u, f, p = integrator @unpack k, tT, uT, duT, uauxT, parse_eqs, rv = cache evaluate!(uT, dt, u) tT[0] = t+dt for i in eachindex(u) @inbounds uT[i][0] = u[i] - duT[i].coeffs .= zero(duT[i][0]) + @inbounds TaylorSeries.zero!(duT[i], 0) end TaylorIntegration.__jetcoeffs!(Val(parse_eqs.x), f, tT, uT, duT, uauxT, p, rv) k = constant_term.(duT) # For the interpolation, needs k at the updated point integrator.stats.nf += 1 end -stepsize_controller!(integrator,alg::TaylorMethodParams) = +OrdinaryDiffEq.stepsize_controller!(integrator,alg::TaylorMethodParams) = TaylorIntegration.stepsize(integrator.cache.uT, integrator.opts.abstol) -step_accept_controller!(integrator, alg::TaylorMethodParams, q) = q +OrdinaryDiffEq.step_accept_controller!(integrator, alg::TaylorMethodParams, q) = q function DiffEqBase.solve( prob::DiffEqBase.AbstractODEProblem{uType, tupType, isinplace}, @@ -237,17 +232,17 @@ function DiffEqBase.solve( _prob = DynamicalODEProblem(f1!, f2!, _u0.x[1], _u0.x[2], prob.tspan, prob.p; prob.kwargs...) else f! = (du, u, p, t) -> (du .= f(u, p, t); 0) - _alg = TaylorMethod(alg.order, parse_eqs = parse_eqs) + _alg = TaylorMethod(alg.order; parse_eqs) _prob = ODEProblem(f!, prob.u0, prob.tspan, prob.p; prob.kwargs...) end else - _alg = TaylorMethod(alg.order, parse_eqs = parse_eqs) + _alg = TaylorMethod(alg.order; parse_eqs) _prob = prob end # DiffEqBase.solve(prob, _alg, args...; kwargs...) integrator = DiffEqBase.__init(_prob, _alg, args...; kwargs...) - integrator.dt = TaylorIntegration.stepsize(integrator.cache.uT, integrator.opts.abstol) # override handle_dt! setting of initial dt + integrator.dt = integrator.tdir * TaylorIntegration.stepsize(integrator.cache.uT, integrator.opts.abstol) # override handle_dt! setting of initial dt DiffEqBase.solve!(integrator) integrator.sol end @@ -256,9 +251,8 @@ end function update_jetcoeffs_cache!(u,f,p,cache::TaylorMethodCache) @unpack tT, uT, duT, uauxT, parse_eqs, rv = cache @inbounds for i in eachindex(u) - uT[i][0] = u[i] - # duT[i].coeffs .= zero(duT[i][0]) - duT[i][0] = zero(uT[i][0]) + @inbounds uT[i][0] = u[i] + @inbounds TaylorSeries.zero!(duT[i], 0) end TaylorIntegration.__jetcoeffs!(Val(parse_eqs.x), f, tT, uT, duT, uauxT, p, rv) return nothing @@ -267,9 +261,9 @@ end # This function was modified from OrdinaryDiffEq.jl; MIT-licensed # _ode_addsteps! overload for ::TaylorMethodCache to handle continuous # and vector callbacks with TaylorIntegration.jl via the common interface -function _ode_addsteps!(k, t, uprev, u, dt, f, p, cache::TaylorMethodCache, +function OrdinaryDiffEq._ode_addsteps!(k, t, uprev, u, dt, f, p, cache::TaylorMethodCache, always_calc_begin = false, allow_calc_end = true,force_calc_end = false) - ### TODO: CHECK, AND IF NECESSARY, RE-SET TIME-STEP SIZE AFTER CALLBACK!!! + ### TODO: check, and if necessary, reset timestep after callback (!) if length(k)<2 || always_calc_begin if typeof(cache) <: OrdinaryDiffEqMutableCache rtmp = similar(u, eltype(eltype(k))) @@ -286,28 +280,66 @@ function _ode_addsteps!(k, t, uprev, u, dt, f, p, cache::TaylorMethodCache, nothing end +function DiffEqBase.interp_summary(::Type{cacheType}, dense::Bool) where {cacheType <: TaylorMethodCache} + dense ? "Taylor series polynomial evaluation" : "1st order linear" +end + +if VERSION < v"1.9" + # used when idxs gives back multiple values + function OrdinaryDiffEq._ode_interpolant!(out, Θ, dt, y₀, y₁, k, + cache::TaylorMethodCache, idxs, T::Type{Val{TI}}) where {TI} + Θm1 = Θ - 1 + @inbounds for i in eachindex(out) + out[i] = cache.uT[i](Θm1*dt) + end + out + end + # used when idxs gives back a single value + function OrdinaryDiffEq._ode_interpolant(Θ, dt, y₀, y₁, k, + cache::TaylorMethodCache, idxs, T::Type{Val{TI}}) where {TI} + Θm1 = Θ - 1 + return cache.uT[idxs](Θm1*dt) + end +else + # used when idxs gives back multiple values + function OrdinaryDiffEq._ode_interpolant!(out, Θ, dt, y₀, y₁, k, + cache::TaylorMethodCache, idxs, T::Type{Val{TI}}, differential_vars) where {TI} + Θm1 = Θ - 1 + @inbounds for i in eachindex(out) + out[i] = cache.uT[i](Θm1*dt) + end + out + end + # used when idxs gives back a single value + function OrdinaryDiffEq._ode_interpolant(Θ, dt, y₀, y₁, k, + cache::TaylorMethodCache, idxs, T::Type{Val{TI}}, differential_vars) where {TI} + Θm1 = Θ - 1 + return cache.uT[idxs](Θm1*dt) + end +end + @inline TaylorIntegration.__jetcoeffs!(::Val{false}, f::ODEFunction, t, x::Taylor1{U}, params, - rv::TaylorIntegration.RetAlloc{Taylor1{U}}) where {U} = TaylorIntegration.__jetcoeffs!(Val(false), f.f, t, x, params) + rv::RetAlloc{Taylor1{U}}) where {U} = TaylorIntegration.__jetcoeffs!(Val(false), f.f, t, x, params, rv) @inline TaylorIntegration.__jetcoeffs!(::Val{true}, f::ODEFunction, t, x::Taylor1{U}, params, - rv::TaylorIntegration.RetAlloc{Taylor1{U}}) where {U} = TaylorIntegration.__jetcoeffs!(Val(true), f.f, t, x, params, rv) + rv::RetAlloc{Taylor1{U}}) where {U} = TaylorIntegration.__jetcoeffs!(Val(true), f.f, t, x, params, rv) @inline TaylorIntegration.__jetcoeffs!(::Val{false}, f::ODEFunction, t, x::AbstractArray{Taylor1{U},N}, dx, xaux, params, - rv::TaylorIntegration.RetAlloc{Taylor1{U}}) where {U,N} = TaylorIntegration.__jetcoeffs!(Val(false), f.f, t, x, dx, xaux, params) + rv::RetAlloc{Taylor1{U}}) where {U,N} = TaylorIntegration.__jetcoeffs!(Val(false), f.f, t, x, dx, xaux, params, rv) @inline TaylorIntegration.__jetcoeffs!(::Val{true}, f::ODEFunction, t, x::AbstractArray{Taylor1{U},N}, dx, xaux, params, - rv::TaylorIntegration.RetAlloc{Taylor1{U}}) where {U,N} = TaylorIntegration.__jetcoeffs!(Val(true), f.f, t, x, dx, params, rv) + rv::RetAlloc{Taylor1{U}}) where {U,N} = TaylorIntegration.__jetcoeffs!(Val(true), f.f, t, x, dx, xaux, params, rv) # NOTE: DynamicalODEFunction assumes x is a vector # @inline TaylorIntegration.__jetcoeffs!(::Val{false}, f::DynamicalODEFunction, t, x::Taylor1{U}, params, -# rv::TaylorIntegration.RetAlloc{Taylor1{U}}) where {U} = TaylorIntegration.__jetcoeffs!(Val(false), f, t, x, params) +# rv::RetAlloc{Taylor1{U}}) where {U} = TaylorIntegration.__jetcoeffs!(Val(false), f, t, x, params) # @inline TaylorIntegration.__jetcoeffs!(::Val{true}, f::DynamicalODEFunction, t, x::Taylor1{U}, params, -# rv::TaylorIntegration.RetAlloc{Taylor1{U}}) where {U} = TaylorIntegration.__jetcoeffs!(Val(true), f, t, x, params, rv) +# rv::RetAlloc{Taylor1{U}}) where {U} = TaylorIntegration.__jetcoeffs!(Val(true), f, t, x, params, rv) @inline TaylorIntegration.__jetcoeffs!(::Val{false}, f::DynamicalODEFunction, t, x::ArrayPartition, dx, xaux, params, - rv::TaylorIntegration.RetAlloc) = TaylorIntegration.jetcoeffs!(f, t, vec(x), vec(dx), xaux, params) + rv::RetAlloc) = TaylorIntegration.jetcoeffs!(f, t, vec(x), vec(dx), xaux, params) # NOTE: `parse_eqs=true` not implemented for `DynamicalODEFunction` # @inline TaylorIntegration.__jetcoeffs!(::Val{true}, f::DynamicalODEFunction, t, x::ArrayPartition, dx, xaux, params, -# rv::TaylorIntegration.RetAlloc) = TaylorIntegration.__jetcoeffs!(Val(true), f, t, vec(x), vec(dx), params, rv) +# rv::RetAlloc) = TaylorIntegration.__jetcoeffs!(Val(true), f, t, vec(x), vec(dx), params, rv) TaylorIntegration._determine_parsing!(parse_eqs::Bool, f::ODEFunction, t, x, params) = TaylorIntegration._determine_parsing!(parse_eqs, f.f, t, x, params) TaylorIntegration._determine_parsing!(parse_eqs::Bool, f::ODEFunction, t, x, dx, params) = TaylorIntegration._determine_parsing!(parse_eqs, f.f, t, x, dx, params) -end \ No newline at end of file +end diff --git a/src/TaylorIntegration.jl b/src/TaylorIntegration.jl index 0686b00d7..64c693fe0 100644 --- a/src/TaylorIntegration.jl +++ b/src/TaylorIntegration.jl @@ -6,15 +6,19 @@ using Reexport @reexport using TaylorSeries using LinearAlgebra using Markdown -using InteractiveUtils: methodswith +using InteractiveUtils #: methodswith if !isdefined(Base, :get_extension) using Requires end -export taylorinteg, lyap_taylorinteg, @taylorize +export TaylorSolution, taylorinteg, lyap_taylorinteg, @taylorize include("parse_eqs.jl") -include("integrator.jl") +include("integrator/jetcoeffs.jl") +include("integrator/stepsize.jl") +include("integrator/taylorstep.jl") +include("integrator/taylorsolution.jl") +include("integrator/taylorinteg.jl") include("lyapunovspectrum.jl") include("rootfinding.jl") include("common.jl") @@ -23,19 +27,19 @@ function __init__() @static if !isdefined(Base, :get_extension) @require OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" begin - include("../ext/TaylorIntegrationDiffEq.jl") + include("../ext/TaylorIntegrationDiffEqExt.jl") end end end @inline function solcoeff!(a::Taylor1{T}, b::Taylor1{T}, k::Int) where {T<:TaylorSeries.NumberNotSeries} - a[k] = b[k-1]/k + @inbounds a[k] = b[k-1]/k return nothing end @inline function solcoeff!(res::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}}, k::Int) where {T<:TaylorSeries.NumberNotSeries} - for l in eachindex(a[k-1]) + @inbounds for l in eachindex(a[k-1]) res[k][l] = a[k-1][l]/k end return nothing @@ -43,7 +47,7 @@ end @inline function solcoeff!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, k::Int) where {T<:TaylorSeries.NumberNotSeries} - for l in eachindex(a[k-1]) + @inbounds for l in eachindex(a[k-1]) for m in eachindex(a[k-1][l]) res[k][l][m] = a[k-1][l][m]/k end diff --git a/src/common.jl b/src/common.jl index a982e06e0..c39108efe 100644 --- a/src/common.jl +++ b/src/common.jl @@ -3,3 +3,4 @@ export TaylorMethod function TaylorMethod end +function update_jetcoeffs_cache! end diff --git a/src/integrator.jl b/src/integrator.jl deleted file mode 100644 index 06fee0965..000000000 --- a/src/integrator.jl +++ /dev/null @@ -1,1053 +0,0 @@ -# This file is part of the TaylorIntegration.jl package; MIT licensed - -# jetcoeffs! -@doc doc""" - jetcoeffs!(eqsdiff::Function, t, x, params) - -Returns an updated `x` using the recursion relation of the -derivatives obtained from the differential equations -``\dot{x}=dx/dt=f(x, p, t)``. - -`eqsdiff` is the function defining the RHS of the ODE, -`x` contains the Taylor1 expansion of the dependent variable(s) and -`t` is the independent variable, and `params` are the parameters appearing on the -function defining the differential equation. See [`taylorinteg`](@ref) for examples -and convention for `eqsdiff`. Note that `x` is of type `Taylor1{U}` where -`U<:Number`; `t` is of type `Taylor1{T}` where `T<:Real`. - -Initially, `x` contains only the 0-th order Taylor coefficient of -the current system state (the initial conditions), and `jetcoeffs!` -computes recursively the high-order derivates back into `x`. - -""" -function jetcoeffs!(eqsdiff::Function, t::Taylor1{T}, x::Taylor1{U}, params) where - {T<:Real, U<:Number} - order = x.order - for ord in 0:order-1 - ordnext = ord+1 - - # # Set `taux`, `xaux`, auxiliary Taylor1 variables to order `ord` - @inbounds xaux = Taylor1( x.coeffs[1:ordnext] ) - - # Equations of motion - dx = eqsdiff(xaux, params, t) - - # Recursion relation - @inbounds solcoeff!(x, dx, ordnext) - end - nothing -end - -@doc doc""" - jetcoeffs!(eqsdiff!::Function, t, x, dx, xaux, params) - -Mutates `x` in-place using the recursion relation of the -derivatives obtained from the differential equations -``\dot{x}=dx/dt=f(x, p, t)``. - -`eqsdiff!` is the function defining the RHS of the ODE, -`x` contains the Taylor1 expansion of the dependent variables and -`t` is the independent variable, and `params` are the parameters appearing on the -function defining the differential equation. See [`taylorinteg`](@ref) for examples -and convention for `eqsdiff`. Note that `x` is of type `Vector{Taylor1{U}}` -where `U<:Number`; `t` is of type `Taylor1{T}` where `T<:Real`. In this case, -two auxiliary containers `dx` and `xaux` (both of the same type as `x`) are -needed to avoid allocations. - -Initially, `x` contains only the 0-th order Taylor coefficient of -the current system state (the initial conditions), and `jetcoeffs!` -computes recursively the high-order derivates back into `x`. - -""" -function jetcoeffs!(eqsdiff!::Function, t::Taylor1{T}, - x::AbstractArray{Taylor1{U}, N}, dx::AbstractArray{Taylor1{U}, N}, - xaux::AbstractArray{Taylor1{U}, N}, params) where {T<:Real, U<:Number, N} - order = x[1].order - for ord in 0:order-1 - ordnext = ord+1 - - # Set `xaux`, auxiliary vector of Taylor1 to order `ord` - for j in eachindex(x) - @inbounds xaux[j] = Taylor1( x[j].coeffs[1:ordnext] ) - end - - # Equations of motion - eqsdiff!(dx, xaux, params, t) - - # Recursion relations - for j in eachindex(x) - @inbounds solcoeff!(x[j], dx[j], ordnext) - end - end - nothing -end - - -# __jetcoeffs -""" - __jetcoeffs!(::Val{false}, f, t, x, params) - __jetcoeffs!(::Val{true}, f, t, x, params, rv) - __jetcoeffs!(::Val{false}, f, t, x, dx, xaux, params) - __jetcoeffs!(::Val{true}, f, t, x, dx, params, rv) - -Chooses a method of [`jetcoeffs!`](@ref) (hard-coded) or the generated by -[`@taylorize`](@ref)) depending on `Val{bool}` (`bool::Bool`). -""" -@inline __jetcoeffs!(::Val{false}, f, t, x::Taylor1{U}, params) where {U} = - jetcoeffs!(f, t, x, params) -@inline __jetcoeffs!(::Val{true}, f, t, x::Taylor1{U}, params, rv::RetAlloc{Taylor1{U}}) where {U} = - jetcoeffs!(Val(f), t, x, params, rv) - -@inline __jetcoeffs!(::Val{false}, f, t, x, dx, xaux, params) = - jetcoeffs!(f, t, x, dx, xaux, params) -@inline __jetcoeffs!(::Val{true}, f, t, x, dx, params, rv) = - jetcoeffs!(Val(f), t, x, dx, params, rv) - - -# _allocate_jetcoeffs! fallbacks; needed to define the corresponding parsed method -@inline _allocate_jetcoeffs!(::Taylor1{T}, x::Taylor1{S}, params) where {T,S} = RetAlloc{Taylor1{S}}() -@inline _allocate_jetcoeffs!(::Taylor1{T}, x::AbstractArray{Taylor1{S}, N}, - ::AbstractArray{Taylor1{S}, N}, params) where {T,S,N} = RetAlloc{Taylor1{S}}() - - -# stepsize -""" - stepsize(x, epsilon) -> h - -Returns a maximum time-step for a the Taylor expansion `x` -using a prescribed absolute tolerance `epsilon` and the last two -Taylor coefficients of (each component of) `x`. - -Note that `x` is of type `Taylor1{U}` or `Vector{Taylor1{U}}`, including -also the cases `Taylor1{TaylorN{U}}` and `Vector{Taylor1{TaylorN{U}}}`. - -Depending of `eltype(x)`, i.e., `U<:Number`, it may be necessary to overload -`stepsize`, specializing it on the type `U`, to avoid type instabilities. -""" -function stepsize(x::Taylor1{U}, epsilon::T) where {T<:Real, U<:Number} - R = promote_type(typeof(norm(constant_term(x), Inf)), T) - ord = x.order - h = convert(R, Inf) - z = zero(R) - for k in (ord-1, ord) - @inbounds aux = norm( x[k], Inf) - aux == z && continue - aux1 = _stepsize(aux, epsilon, k) - h = min(h, aux1) - end - return h::R -end - -function stepsize(q::AbstractArray{Taylor1{U}, N}, epsilon::T) where - {T<:Real, U<:Number, N} - R = promote_type(typeof(norm(constant_term(q[1]), Inf)), T) - h = convert(R, Inf) - for i in eachindex(q) - @inbounds hi = stepsize( q[i], epsilon ) - h = min( h, hi ) - end - - # If `isinf(h)==true`, we use the maximum (finite) - # step-size obtained from all coefficients as above. - # Note that the time step is independent from `epsilon`. - if isinf(h) - h = zero(R) - for i in eachindex(q) - @inbounds hi = _second_stepsize(q[i], epsilon) - h = max( h, hi ) - end - end - return h::R -end - -""" - _stepsize(aux1, epsilon, k) - -Helper function to avoid code repetition. -Returns ``(epsilon/aux1)^(1/k)``. -""" -@inline function _stepsize(aux1::U, epsilon::T, k::Int) where {T<:Real, U<:Number} - aux = epsilon / aux1 - kinv = 1 / k - return aux^kinv -end - -""" - _second_stepsize(x, epsilon) - -Corresponds to the "second stepsize control" in Jorba and Zou -(2005) paper. We use it if [`stepsize`](@ref) returns `Inf`. -""" -function _second_stepsize(x::Taylor1{U}, epsilon::T) where {T<:Real, U<:Number} - R = promote_type(typeof(norm(constant_term(x), Inf)), T) - x == zero(x) && return convert(R, Inf) - ord = x.order - z = zero(R) - u = one(R) - h = z - for k in 1:ord-2 - @inbounds aux = norm( x[k], Inf) - aux == z && continue - aux1 = _stepsize(aux, u, k) - h = max(h, aux1) - end - return h::R -end - -#taylorstep -@doc doc""" - taylorstep!(f, t, x, t0, order, abstol, params, tmpTaylor, arrTaylor, parse_eqs=true) -> δt - taylorstep!(f!, t, x, dx, xaux, t0, order, abstol, params, tmpTaylor, arrTaylor, parse_eqs=true) -> δt - -One-step Taylor integration for the one-dependent variable ODE ``\dot{x}=dx/dt=f(x, p, t)`` -with initial conditions ``x(t_0)=x_0``. -Returns the time-step `δt` of the actual integration carried out (δt is positive). - -Here, `f` is the function defining the RHS of the ODE (see -[`taylorinteg`](@ref)), `t` is the -independent variable, `x` contains the Taylor expansion of the dependent -variable, `order` is -the degree used for the `Taylor1` polynomials during the integration -`abstol` is the absolute tolerance used to determine the time step -of the integration, and `params` are the parameters entering the ODE -functions. -For several variables, `dx` and `xaux`, both of the same type as `x`, -are needed to save allocations. Finally, `parse_eqs` is a switch -to force *not* using (`parse_eqs=false`) the specialized method of `jetcoeffs!` -created with [`@taylorize`](@ref); the default is `true` (parse the equations). -Finally, `parse_eqs` is a switch -to force *not* using (`parse_eqs=false`) the specialized method of `jetcoeffs!` -created with [`@taylorize`](@ref); the default is `true` (parse the equations). - -""" -function taylorstep!(f, t::Taylor1{T}, x::Taylor1{U}, abstol::T, params) where {T<:Real, U<:Number} - - # Compute the Taylor coefficients - __jetcoeffs!(Val(false), f, t, x, params) - - # Compute the step-size of the integration using `abstol` - δt = stepsize(x, abstol) - if isinf(δt) - δt = _second_stepsize(x, abstol) - end - - return δt -end - -function taylorstep!(f, t::Taylor1{T}, x::Taylor1{U}, abstol::T, params, - rv::RetAlloc{Taylor1{U}}) where {T<:Real, U<:Number} - - # Compute the Taylor coefficients - __jetcoeffs!(Val(true), f, t, x, params, rv) - - # Compute the step-size of the integration using `abstol` - δt = stepsize(x, abstol) - if isinf(δt) - δt = _second_stepsize(x, abstol) - end - - return δt -end - -function taylorstep!(f!, t::Taylor1{T}, x::Vector{Taylor1{U}}, - dx::Vector{Taylor1{U}}, xaux::Vector{Taylor1{U}}, abstol::T, params) where - {T<:Real, U<:Number} - - # Compute the Taylor coefficients - __jetcoeffs!(Val(false), f!, t, x, dx, xaux, params) - - # Compute the step-size of the integration using `abstol` - δt = stepsize(x, abstol) - - return δt -end - -function taylorstep!(f!, t::Taylor1{T}, x::Vector{Taylor1{U}}, - dx::Vector{Taylor1{U}}, abstol::T, params, - rv::RetAlloc{Taylor1{U}}) where {T<:Real, U<:Number} - - # Compute the Taylor coefficients - __jetcoeffs!(Val(true), f!, t, x, dx, params, rv) - - # Compute the step-size of the integration using `abstol` - δt = stepsize(x, abstol) - - return δt -end - - - -""" - _determine_parsing!(parse_eqs::Bool, f, t, x, params) - _determine_parsing!(parse_eqs::Bool, f, t, x, dx, params) - -Check if the parsed method of `jetcoeffs!` exists and check it -runs without error. -""" -function _determine_parsing!(parse_eqs::Bool, f, t, x, params) - parse_eqs = parse_eqs && !isempty(methodswith(Val{f}, jetcoeffs!)) && - !isempty(methodswith(Val{f}, _allocate_jetcoeffs!)) - - # tmpTaylor, arrTaylor = _allocate_jetcoeffs!(t, x, params) - rv = _allocate_jetcoeffs!(t, x, params) - - if parse_eqs - try - rv = _allocate_jetcoeffs!(Val(f), t, x, params) - jetcoeffs!(Val(f), t, x, params, rv) - catch - @warn("""Unable to use the parsed method of `jetcoeffs!` for `$f`, - despite of having `parse_eqs=true`, due to some internal error. - Using `parse_eqs = false`.""") - parse_eqs = false - end - end - - return parse_eqs, rv -end -function _determine_parsing!(parse_eqs::Bool, f, t, x, dx, params) - - parse_eqs = parse_eqs && !isempty(methodswith(Val{f}, jetcoeffs!)) && - !isempty(methodswith(Val{f}, _allocate_jetcoeffs!)) - - rv = _allocate_jetcoeffs!(t, x, dx, params) - - if parse_eqs - try - rv = _allocate_jetcoeffs!(Val(f), t, x, dx, params) - jetcoeffs!(Val(f), t, x, dx, params, rv) - catch - @warn("""Unable to use the parsed method of `jetcoeffs!` for `$f`, - despite of having `parse_eqs=true`, due to some internal error. - Using `parse_eqs = false`.""") - parse_eqs = false - end - end - - return parse_eqs, rv -end - - -# taylorinteg -for V in (:(Val{true}), :(Val{false})) - @eval begin - function taylorinteg(f, x0::U, t0::T, tmax::T, order::Int, abstol::T, ::$V, params = nothing; - maxsteps::Int=500, parse_eqs::Bool=true) where {T<:Real, U<:Number} - - # Initialize the Taylor1 expansions - t = t0 + Taylor1( T, order ) - x = Taylor1( x0, order ) - - # Determine if specialized jetcoeffs! method exists - parse_eqs, rv = _determine_parsing!(parse_eqs, f, t, x, params) - - if parse_eqs - # Re-initialize the Taylor1 expansions - t = t0 + Taylor1( T, order ) - x = Taylor1( x0, order ) - return _taylorinteg!(f, t, x, x0, t0, tmax, abstol, rv, - $V(), params, maxsteps=maxsteps) - else - return _taylorinteg!(f, t, x, x0, t0, tmax, abstol, - $V(), params, maxsteps=maxsteps) - end - end - - function _taylorinteg!(f, t::Taylor1{T}, x::Taylor1{U}, - x0::U, t0::T, tmax::T, abstol::T, ::$V, params; - maxsteps::Int=500) where {T<:Real, U<:Number} - - # Allocation - tv = Array{T}(undef, maxsteps+1) - xv = Array{U}(undef, maxsteps+1) - if $V == Val{true} - psol = Array{Taylor1{U}}(undef, maxsteps) - end - - # Initial conditions - nsteps = 1 - @inbounds t[0] = t0 - @inbounds tv[1] = t0 - @inbounds xv[1] = x0 - sign_tstep = copysign(1, tmax-t0) - - # Integration - while sign_tstep*t0 < sign_tstep*tmax - δt = taylorstep!(f, t, x, abstol, params) # δ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 - if $V == Val{true} - @inbounds psol[nsteps] = deepcopy(x) - end - @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 - - if $V == Val{true} - return view(tv,1:nsteps), view(xv,1:nsteps), view(psol, 1:nsteps-1) - elseif $V == Val{false} - return view(tv,1:nsteps), view(xv,1:nsteps) - end - end - function _taylorinteg!(f, t::Taylor1{T}, x::Taylor1{U}, - x0::U, t0::T, tmax::T, abstol::T, rv::RetAlloc{Taylor1{U}}, ::$V, params; - maxsteps::Int=500) where {T<:Real, U<:Number} - - # Allocation - tv = Array{T}(undef, maxsteps+1) - xv = Array{U}(undef, maxsteps+1) - if $V == Val{true} - psol = Array{Taylor1{U}}(undef, maxsteps) - end - - # Initial conditions - nsteps = 1 - @inbounds t[0] = t0 - @inbounds tv[1] = t0 - @inbounds xv[1] = x0 - sign_tstep = copysign(1, tmax-t0) - - # Integration - while sign_tstep*t0 < sign_tstep*tmax - δt = taylorstep!(f, t, x, abstol, params, rv) # δ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 - if $V == Val{true} - @inbounds psol[nsteps] = deepcopy(x) - end - @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 - - if $V == Val{true} - return view(tv,1:nsteps), view(xv,1:nsteps), view(psol, 1:nsteps-1) - elseif $V == Val{false} - return view(tv,1:nsteps), view(xv,1:nsteps) - end - end - - - function taylorinteg(f!, q0::Array{U,1}, t0::T, tmax::T, order::Int, abstol::T, ::$V, params = nothing; - maxsteps::Int=500, parse_eqs::Bool=true) where {T<:Real, U<:Number} - - # Initialize the vector of Taylor1 expansions - dof = length(q0) - t = t0 + Taylor1( T, order ) - x = Array{Taylor1{U}}(undef, dof) - dx = Array{Taylor1{U}}(undef, dof) - @inbounds for i in eachindex(q0) - @inbounds x[i] = Taylor1( q0[i], order ) - @inbounds dx[i] = Taylor1( zero(q0[i]), order ) - end - - # Determine if specialized jetcoeffs! method exists - parse_eqs, rv = _determine_parsing!(parse_eqs, f!, t, x, dx, params) - - if parse_eqs - # Re-initialize the Taylor1 expansions - t = t0 + Taylor1( T, order ) - x .= Taylor1.( q0, order ) - dx .= Taylor1.( zero.(q0), order) - return _taylorinteg!(f!, t, x, dx, q0, t0, tmax, abstol, rv, - $V(), params; maxsteps) - else - return _taylorinteg!(f!, t, x, dx, q0, t0, tmax, abstol, - $V(), params; maxsteps) - end - end - - function _taylorinteg!(f!, t::Taylor1{T}, x::Array{Taylor1{U},1}, dx::Array{Taylor1{U},1}, - q0::Array{U,1}, t0::T, tmax::T, abstol::T, ::$V, params; - maxsteps::Int=500) where {T<:Real, U<:Number} - - # Initialize the vector of Taylor1 expansions - dof = length(q0) - - # Allocation - tv = Array{T}(undef, maxsteps+1) - xv = Array{U}(undef, dof, maxsteps+1) - if $V == Val{true} - psol = Array{Taylor1{U}}(undef, dof, maxsteps) - end - xaux = Array{Taylor1{U}}(undef, dof) - - # Initial conditions - @inbounds t[0] = t0 - # x .= Taylor1.(q0, order) - x0 = deepcopy(q0) - @inbounds tv[1] = t0 - @inbounds xv[:,1] .= q0 - sign_tstep = copysign(1, tmax-t0) - - # Integration - nsteps = 1 - while sign_tstep*t0 < sign_tstep*tmax - δt = taylorstep!(f!, t, x, dx, xaux, abstol, params) # δ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 - if $V == Val{true} - # Store the Taylor polynomial solution - @inbounds psol[:,nsteps] .= deepcopy.(x) - end - @inbounds for i in eachindex(x0) - x[i][0] = x0[i] - dx[i][0] = zero(x0[i]) - end - t0 += δt - @inbounds t[0] = t0 - nsteps += 1 - @inbounds tv[nsteps] = t0 - @inbounds xv[:,nsteps] .= deepcopy.(x0) - if nsteps > maxsteps - @warn(""" - Maximum number of integration steps reached; exiting. - """) - break - end - end - - if $V == Val{true} - return view(tv,1:nsteps), view(transpose(view(xv,:,1:nsteps)),1:nsteps,:), - view(transpose(view(psol, :, 1:nsteps-1)), 1:nsteps-1, :) - elseif $V == Val{false} - return view(tv,1:nsteps), view(transpose(view(xv,:,1:nsteps)),1:nsteps,:) - end - end - - function _taylorinteg!(f!, t::Taylor1{T}, x::Array{Taylor1{U},1}, dx::Array{Taylor1{U},1}, - q0::Array{U,1}, t0::T, tmax::T, abstol::T, rv::RetAlloc{Taylor1{U}}, ::$V, params; - maxsteps::Int=500) where {T<:Real, U<:Number} - - # Initialize the vector of Taylor1 expansions - dof = length(q0) - - # Allocation of output - tv = Array{T}(undef, maxsteps+1) - xv = Array{U}(undef, dof, maxsteps+1) - if $V == Val{true} - psol = Array{Taylor1{U}}(undef, dof, maxsteps) - end - - # Initial conditions - @inbounds t[0] = t0 - x0 = deepcopy(q0) - @inbounds tv[1] = t0 - @inbounds xv[:,1] .= q0 - sign_tstep = copysign(1, tmax-t0) - - # Integration - nsteps = 1 - while sign_tstep*t0 < sign_tstep*tmax - δt = taylorstep!(f!, t, x, dx, abstol, params, rv) # δ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 - if $V == Val{true} - # Store the Taylor polynomial solution - @inbounds psol[:,nsteps] .= deepcopy.(x) - end - @inbounds for i in eachindex(x0) - x[i][0] = x0[i] - dx[i][0] = zero(x0[i]) - end - t0 += δt - @inbounds t[0] = t0 - nsteps += 1 - @inbounds tv[nsteps] = t0 - @inbounds xv[:,nsteps] .= deepcopy.(x0) - if nsteps > maxsteps - @warn(""" - Maximum number of integration steps reached; exiting. - """) - break - end - end - - if $V == Val{true} - return view(tv,1:nsteps), view(transpose(view(xv,:,1:nsteps)),1:nsteps,:), - view(transpose(view(psol, :, 1:nsteps-1)), 1:nsteps-1, :) - elseif $V == Val{false} - return view(tv,1:nsteps), view(transpose(view(xv,:,1:nsteps)),1:nsteps,:) - end - end - - end -end -@doc doc""" - taylorinteg(f, x0, t0, tmax, order, abstol, params[=nothing]; kwargs... ) - taylorinteg(f, x0, t0, tmax, order, abstol, Val(false), 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)``, -where `p` are the parameters encoded in `params`. -The initial conditions are specified by `x0` at time `t0`; `x0` may be of type `T<:Number` -or `Vector{T}`, with `T` including `TaylorN{T}`; the latter case -is of interest for [jet transport applications](@ref jettransport). - -The equations of motion are specified by the function `f`; we follow the same -convention of `DifferentialEquations.jl` to define this function, i.e., -`f(x, p, t)` or `f!(dx, x, p, t)`; see the examples below. - -The functions return a vector with the values of time (independent variable), -and a vector with the computed values of -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 `tmax`, or when the number of saved steps is larger -than `maxsteps`. - -The integration method uses polynomial expansions on the independent variable -of order `order`; the parameter `abstol` serves to define the -time step using the last two Taylor coefficients of the expansions. -Make sure you use a *large enough* `order` to assure convergence. - -Currently, the recognized keyword arguments are: -- `maxsteps[=500]`: maximum number of integration steps. -- `parse_eqs[=true]`: use the specialized method of `jetcoeffs!` created - with [`@taylorize`](@ref). - -## Examples - -For one dependent variable the function `f` defines the RHS of the equation of -motion, returning the value of ``\dot{x}``. The arguments of -this function are `(x, p, t)`, where `x` are the dependent variables, `p` are -the paremeters and `t` is the independent variable. - -For several (two or more) dependent variables, the function `f!` defines -the RHS of the equations of motion, mutating (in-place) the (preallocated) vector -with components of ``\dot{x}``. The arguments of this function are `(dx, x, p, t)`, -where `dx` is the preallocated vector of ``\dot{x}``, `x` are the dependent -variables, `p` are the paremeters entering the ODEs and `t` is the independent -variable. The function may return this vector or simply `nothing`. - -```julia -using TaylorIntegration - -f(x, p, t) = x^2 - -tv, xv = taylorinteg(f, 3, 0.0, 0.3, 25, 1.0e-20, maxsteps=100 ) - -function f!(dx, x, p, t) - for i in eachindex(x) - dx[i] = x[i]^2 - end - return nothing -end - -tv, xv = taylorinteg(f!, [3, 3], 0.0, 0.3, 25, 1.0e-20, maxsteps=100 ) - -tv, xv, psol = taylorinteg(f!, [3, 3], 0.0, 0.3, 25, 1.0e-20, maxsteps=100, Val(true) ) -``` - -""" taylorinteg - - -# Integrate and return results evaluated at given time -@doc doc""" - taylorinteg(f, x0, trange, order, abstol, params[=nothing]; kwargs... ) - -General-purpose Taylor integrator for the explicit ODE -``\dot{x}=f(x,p,t)`` with initial condition specified by `x0::{T<:Number}` -or `x0::Vector{T}` at time `t0`. - -The method returns a vector with the computed values of -the dependent variable(s), evaluated *only* at the times specified by -the range `trange`. - -## Examples - -```julia -xv = taylorinteg(f, 3, 0.0:0.001:0.3, 25, 1.0e-20, maxsteps=100 ) - -xv = taylorinteg(f!, [3, 3], 0.0:0.001:0.3, 25, 1.0e-20, maxsteps=100 ); - -``` - -""" -function taylorinteg(f, x0::U, trange::AbstractVector{T}, - order::Int, abstol::T, params = nothing; - maxsteps::Int=500, parse_eqs::Bool=true) where {T<:Real, U<:Number} - - # Check if trange is increasingly or decreasingly sorted - @assert (issorted(trange) || - issorted(reverse(trange))) "`trange` or `reverse(trange)` must be sorted" - - # Initialize the Taylor1 expansions - t0 = trange[1] - t = t0 + Taylor1( T, order ) - x = Taylor1( x0, order ) - - # Determine if specialized jetcoeffs! method exists - parse_eqs, rv = _determine_parsing!(parse_eqs, f, t, x, params) - - if parse_eqs - # Re-initialize the Taylor1 expansions - t = t0 + Taylor1( T, order ) - x = Taylor1( x0, order ) - return _taylorinteg!(f, t, x, x0, trange, abstol, rv, - params, maxsteps=maxsteps) - else - return _taylorinteg!(f, t, x, x0, trange, abstol, params, maxsteps=maxsteps) - end -end - -function _taylorinteg!(f, t::Taylor1{T}, x::Taylor1{U}, x0::U, trange::AbstractVector{T}, - abstol::T, params; maxsteps::Int=500) where {T<:Real, U<:Number} - - # Allocation - nn = length(trange) - xv = Array{U}(undef, nn) - fill!(xv, T(NaN)) - - # Initial conditions - @inbounds t0, t1, tmax = trange[1], trange[2], trange[end] - sign_tstep = copysign(1, tmax-t0) - @inbounds t[0] = t0 - @inbounds xv[1] = x0 - - # Integration - iter = 2 - nsteps = 1 - while sign_tstep*t0 < sign_tstep*tmax - δt = taylorstep!(f, t, x, abstol, params)# δ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 - tnext = t0+δt - # Evaluate solution at times within convergence radius - while sign_tstep*t1 < sign_tstep*tnext - x1 = evaluate(x, t1-t0) - @inbounds xv[iter] = x1 - iter += 1 - @inbounds t1 = trange[iter] - end - if δt == tmax-t0 - @inbounds xv[iter] = x0 - break - end - @inbounds x[0] = x0 - t0 = tnext - @inbounds t[0] = t0 - nsteps += 1 - if nsteps > maxsteps - @warn(""" - Maximum number of integration steps reached; exiting. - """) - break - end - end - return xv -end -function _taylorinteg!(f, t::Taylor1{T}, x::Taylor1{U}, x0::U, trange::AbstractVector{T}, - abstol::T, rv::RetAlloc{Taylor1{U}}, params; maxsteps::Int=500) where {T<:Real, U<:Number} - - # Allocation - nn = length(trange) - xv = Array{U}(undef, nn) - fill!(xv, T(NaN)) - - # Initial conditions - @inbounds t0, t1, tmax = trange[1], trange[2], trange[end] - sign_tstep = copysign(1, tmax-t0) - @inbounds t[0] = t0 - @inbounds xv[1] = x0 - - # Integration - iter = 2 - nsteps = 1 - while sign_tstep*t0 < sign_tstep*tmax - δt = taylorstep!(f, t, x, abstol, params, rv)# δ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 - tnext = t0+δt - # Evaluate solution at times within convergence radius - while sign_tstep*t1 < sign_tstep*tnext - x1 = evaluate(x, t1-t0) - @inbounds xv[iter] = x1 - iter += 1 - @inbounds t1 = trange[iter] - end - if δt == tmax-t0 - @inbounds xv[iter] = x0 - break - end - @inbounds x[0] = x0 - t0 = tnext - @inbounds t[0] = t0 - nsteps += 1 - if nsteps > maxsteps - @warn(""" - Maximum number of integration steps reached; exiting. - """) - break - end - end - return xv -end - -function taylorinteg(f!, q0::Array{U,1}, trange::AbstractVector{T}, - order::Int, abstol::T, params = nothing; - maxsteps::Int=500, parse_eqs::Bool=true) where {T<:Real, U<:Number} - - # Check if trange is increasingly or decreasingly sorted - @assert (issorted(trange) || - issorted(reverse(trange))) "`trange` or `reverse(trange)` must be sorted" - - # Initialize the vector of Taylor1 expansions - dof = length(q0) - t0 = trange[1] - t = t0 + Taylor1( T, order ) - x = Array{Taylor1{U}}(undef, dof) - dx = Array{Taylor1{U}}(undef, dof) - x .= Taylor1.( q0, order ) - dx .= Taylor1.( zero.(q0), order ) - - # Determine if specialized jetcoeffs! method exists - parse_eqs, rv = _determine_parsing!(parse_eqs, f!, t, x, dx, params) - - if parse_eqs - # Re-initialize the Taylor1 expansions - t = t0 + Taylor1( T, order ) - x .= Taylor1.( q0, order ) - dx .= Taylor1.( zero.(q0), order ) - return _taylorinteg!(f!, t, x, dx, q0, trange, abstol, rv, - params, maxsteps=maxsteps) - else - return _taylorinteg!(f!, t, x, dx, q0, trange, abstol, params, maxsteps=maxsteps) - end -end - -function _taylorinteg!(f!, t::Taylor1{T}, x::Array{Taylor1{U},1}, dx::Array{Taylor1{U},1}, - q0::Array{U,1}, trange::AbstractVector{T}, abstol::T, params; - maxsteps::Int=500) where {T<:Real, U<:Number} - - # Allocation - nn = length(trange) - dof = length(q0) - x0 = similar(q0, eltype(q0), dof) - x1 = similar(x0) - fill!(x0, T(NaN)) - xv = Array{eltype(q0)}(undef, dof, nn) - for ind in 1:nn - @inbounds xv[:,ind] .= x0 - end - xaux = Array{Taylor1{U}}(undef, dof) - - # Initial conditions - @inbounds t[0] = trange[1] - @inbounds t0, t1, tmax = trange[1], trange[2], trange[end] - sign_tstep = copysign(1, tmax-t0) - # x .= Taylor1.(q0, order) - @inbounds x0 .= q0 - @inbounds xv[:,1] .= q0 - - # Integration - iter = 2 - nsteps = 1 - while sign_tstep*t0 < sign_tstep*tmax - δt = taylorstep!(f!, t, x, dx, xaux, abstol, params) # δ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 - tnext = t0+δt - # Evaluate solution at times within convergence radius - while sign_tstep*t1 < sign_tstep*tnext - evaluate!(x, t1-t0, x1) - @inbounds xv[:,iter] .= x1 - iter += 1 - @inbounds t1 = trange[iter] - end - if δt == tmax-t0 - @inbounds xv[:,iter] .= x0 - break - end - @inbounds for i in eachindex(x0) - x[i][0] = x0[i] - dx[i][0] = zero(x0[i]) - end - t0 = tnext - @inbounds t[0] = t0 - nsteps += 1 - if nsteps > maxsteps - @warn(""" - Maximum number of integration steps reached; exiting. - """) - break - end - end - - return transpose(xv) -end -function _taylorinteg!(f!, t::Taylor1{T}, x::Array{Taylor1{U},1}, dx::Array{Taylor1{U},1}, - q0::Array{U,1}, trange::AbstractVector{T}, abstol::T, rv::RetAlloc{Taylor1{U}}, params; - maxsteps::Int=500) where {T<:Real, U<:Number} - - # Allocation - nn = length(trange) - dof = length(q0) - x0 = similar(q0, eltype(q0), dof) - x1 = similar(x0) - fill!(x0, T(NaN)) - xv = Array{eltype(q0)}(undef, dof, nn) - for ind in 1:nn - @inbounds xv[:,ind] .= x0 - end - - # Initial conditions - @inbounds t[0] = trange[1] - @inbounds t0, t1, tmax = trange[1], trange[2], trange[end] - sign_tstep = copysign(1, tmax-t0) - # x .= Taylor1.(q0, order) - @inbounds x0 .= q0 - @inbounds xv[:,1] .= q0 - - # Integration - iter = 2 - nsteps = 1 - while sign_tstep*t0 < sign_tstep*tmax - δt = taylorstep!(f!, t, x, dx, abstol, params, rv) # δ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 - tnext = t0+δt - # Evaluate solution at times within convergence radius - while sign_tstep*t1 < sign_tstep*tnext - evaluate!(x, t1-t0, x1) - @inbounds xv[:,iter] .= x1 - iter += 1 - @inbounds t1 = trange[iter] - end - if δt == tmax-t0 - @inbounds xv[:,iter] .= x0 - break - end - @inbounds for i in eachindex(x0) - x[i][0] = x0[i] - dx[i][0] = zero(x0[i]) - end - t0 = tnext - @inbounds t[0] = t0 - nsteps += 1 - if nsteps > maxsteps - @warn(""" - Maximum number of integration steps reached; exiting. - """) - break - end - end - - return transpose(xv) -end - - -# Generic functions -for R in (:Number, :Integer) - @eval begin - taylorinteg(f, xx0::S, tt0::T, ttmax::U, order::Int, aabstol::V, params = nothing; - maxsteps::Int=500, parse_eqs::Bool=true) where {S<:$R, T<:Real, U<:Real, V<:Real} = - taylorinteg(f, xx0, tt0, ttmax, order, aabstol, Val(false), params, - maxsteps=maxsteps, parse_eqs=parse_eqs) - - taylorinteg(f, q0::Array{S,1}, tt0::T, ttmax::U, order::Int, aabstol::V, params = nothing; - maxsteps::Int=500, parse_eqs::Bool=true) where {S<:$R, T<:Real, U<:Real, V<:Real} = - taylorinteg(f, q0, tt0, ttmax, order, aabstol, Val(false), params, - maxsteps=maxsteps, parse_eqs=parse_eqs) - - function taylorinteg(f, xx0::S, trange::AbstractVector{T}, order::Int, aabstol::U, params = nothing; - maxsteps::Int=500, parse_eqs::Bool=true) where {S<:$R, T<:Real, U<:Real} - - t0, abstol, _ = promote(trange[1], aabstol, one(Float64)) - x0, _ = promote(xx0, t0) - - return taylorinteg(f, x0, trange .* one(t0), order, abstol, params, - maxsteps=maxsteps, parse_eqs=parse_eqs) - end - - function taylorinteg(f, q0::Array{S,1}, trange::AbstractVector{T}, order::Int, aabstol::U, params = nothing; - maxsteps::Int=500, parse_eqs::Bool=true) where {S<:$R, T<:Real, U<:Real} - - t0, abstol, _ = promote(trange[1], aabstol, one(Float64)) - elq0, _ = promote(q0[1], t0) - q0_ = convert(Array{typeof(elq0)}, q0) - - return taylorinteg(f, q0_, 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, _ = promote(tt0, ttmax, aabstol, one(Float64)) - x0, _ = 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, _ = promote(tt0, ttmax, aabstol, one(Float64)) - elq0, _ = 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 - - function taylorinteg(f, xx0::S, tt0::T, ttmax::U, order::Int, aabstol::V, - ::Val{false}, 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, _ = promote(tt0, ttmax, aabstol, one(Float64)) - x0, _ = promote(xx0, t0) - - return taylorinteg(f, x0, t0, tmax, order, abstol, Val(false), params, - maxsteps=maxsteps, parse_eqs=parse_eqs) - end - - function taylorinteg(f, q0::Array{S,1}, tt0::T, ttmax::U, order::Int, aabstol::V, - ::Val{false}, 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, _ = promote(tt0, ttmax, aabstol, one(Float64)) - elq0, _ = 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(false), params, - maxsteps=maxsteps, parse_eqs=parse_eqs) - end - end -end diff --git a/src/integrator/jetcoeffs.jl b/src/integrator/jetcoeffs.jl new file mode 100644 index 000000000..034f5a57b --- /dev/null +++ b/src/integrator/jetcoeffs.jl @@ -0,0 +1,162 @@ +# This file is part of the TaylorIntegration.jl package; MIT licensed + +# jetcoeffs! +@doc doc""" + jetcoeffs!(eqsdiff::Function, t, x, params) + +Returns an updated `x` using the recursion relation of the +derivatives obtained from the differential equations +``\dot{x}=dx/dt=f(x, p, t)``. + +`eqsdiff` is the function defining the RHS of the ODE, +`x` contains the Taylor1 expansion of the dependent variable(s) and +`t` is the independent variable, and `params` are the parameters appearing on the +function defining the differential equation. See [`taylorinteg`](@ref) for examples +and convention for `eqsdiff`. Note that `x` is of type `Taylor1{U}` where +`U<:Number`; `t` is of type `Taylor1{T}` where `T<:Real`. + +Initially, `x` contains only the 0-th order Taylor coefficient of +the current system state (the initial conditions), and `jetcoeffs!` +computes recursively the high-order derivates back into `x`. + +""" +function jetcoeffs!(eqsdiff::Function, t::Taylor1{T}, x::Taylor1{U}, params) where + {T<:Real, U<:Number} + order = x.order + for ord in 0:order-1 + ordnext = ord+1 + + # Set `xaux` auxiliary Taylor1 variables to order `ord` + @inbounds xaux = Taylor1( x.coeffs[1:ordnext] ) + + # Equations of motion + dx = eqsdiff(xaux, params, t) + + # Recursion relation + solcoeff!(x, dx, ordnext) + end + nothing +end + +@doc doc""" + jetcoeffs!(eqsdiff!::Function, t, x, dx, xaux, params) + +Mutates `x` in-place using the recursion relation of the +derivatives obtained from the differential equations +``\dot{x}=dx/dt=f(x, p, t)``. + +`eqsdiff!` is the function defining the RHS of the ODE, +`x` contains the Taylor1 expansion of the dependent variables and +`t` is the independent variable, and `params` are the parameters appearing on the +function defining the differential equation. See [`taylorinteg`](@ref) for examples +and convention for `eqsdiff`. Note that `x` is of type `Vector{Taylor1{U}}` +where `U<:Number`; `t` is of type `Taylor1{T}` where `T<:Real`. In this case, +two auxiliary containers `dx` and `xaux` (both of the same type as `x`) are +needed to avoid allocations. + +Initially, `x` contains only the 0-th order Taylor coefficient of +the current system state (the initial conditions), and `jetcoeffs!` +computes recursively the high-order derivates back into `x`. + +""" +function jetcoeffs!(eqsdiff!::Function, t::Taylor1{T}, + x::AbstractArray{Taylor1{U}, N}, dx::AbstractArray{Taylor1{U}, N}, + xaux::AbstractArray{Taylor1{U}, N}, params) where {T<:Real, U<:Number, N} + order = x[1].order + for ord in 0:order-1 + ordnext = ord+1 + + # Set `xaux`, auxiliary vector of Taylor1 to order `ord` + for j in eachindex(x) + @inbounds xaux[j] = Taylor1( x[j].coeffs[1:ordnext] ) + end + + # Equations of motion + eqsdiff!(dx, xaux, params, t) + + # Recursion relations + for j in eachindex(x) + @inbounds solcoeff!(x[j], dx[j], ordnext) + end + end + nothing +end + + +# __jetcoeffs +""" + __jetcoeffs!(::Val{false}, f, t, x, params, rv) + __jetcoeffs!(::Val{true}, f, t, x, params, rv) + __jetcoeffs!(::Val{false}, f, t, x, dx, xaux, params, rv) + __jetcoeffs!(::Val{true}, f, t, x, dx, xaux, params, rv) + +Chooses a method of [`jetcoeffs!`](@ref) (hard-coded) or the generated by +[`@taylorize`](@ref)) depending on `Val{bool}` (`bool::Bool`). +""" +@inline __jetcoeffs!(::Val{false}, f, t, x::Taylor1{U}, params, ::RetAlloc{Taylor1{U}}) where {U} = + jetcoeffs!(f, t, x, params) +@inline __jetcoeffs!(::Val{true}, f, t, x::Taylor1{U}, params, rv::RetAlloc{Taylor1{U}}) where {U} = + jetcoeffs!(Val(f), t, x, params, rv) + +@inline __jetcoeffs!(::Val{false}, f, t, x, dx, xaux, params, rv) = + jetcoeffs!(f, t, x, dx, xaux, params) +@inline __jetcoeffs!(::Val{true}, f, t, x, dx, xaux, params, rv) = + jetcoeffs!(Val(f), t, x, dx, params, rv) + + +# _allocate_jetcoeffs! fallbacks; needed to define the corresponding parsed method +@inline _allocate_jetcoeffs!(::Taylor1{T}, x::Taylor1{S}, params) where {T,S} = RetAlloc{Taylor1{S}}() +@inline _allocate_jetcoeffs!(::Taylor1{T}, x::AbstractArray{Taylor1{S}, N}, + ::AbstractArray{Taylor1{S}, N}, params) where {T,S,N} = RetAlloc{Taylor1{S}}() + + +# _determine_parsing! +""" + _determine_parsing!(parse_eqs::Bool, f, t, x, params) + _determine_parsing!(parse_eqs::Bool, f, t, x, dx, params) + +Check if the parsed method of `jetcoeffs!` exists and check it +runs without error. +""" +function _determine_parsing!(parse_eqs::Bool, f, t, x, params) + parse_eqs = parse_eqs && !isempty(methodswith(Val{f}, jetcoeffs!)) && + !isempty(methodswith(Val{f}, _allocate_jetcoeffs!)) + + # tmpTaylor, arrTaylor = _allocate_jetcoeffs!(t, x, params) + rv = _allocate_jetcoeffs!(t, x, params) + + if parse_eqs + try + rv = _allocate_jetcoeffs!(Val(f), t, x, params) + jetcoeffs!(Val(f), t, x, params, rv) + catch + @warn("""Unable to use the parsed method of `jetcoeffs!` for `$f`, + despite of having `parse_eqs=true`, due to some internal error. + Using `parse_eqs = false`.""") + parse_eqs = false + end + end + + return parse_eqs, rv +end +function _determine_parsing!(parse_eqs::Bool, f, t, x, dx, params) + + parse_eqs = parse_eqs && !isempty(methodswith(Val{f}, jetcoeffs!)) && + !isempty(methodswith(Val{f}, _allocate_jetcoeffs!)) + + rv = _allocate_jetcoeffs!(t, x, dx, params) + + if parse_eqs + try + rv = _allocate_jetcoeffs!(Val(f), t, x, dx, params) + jetcoeffs!(Val(f), t, x, dx, params, rv) + catch + @warn("""Unable to use the parsed method of `jetcoeffs!` for `$f`, + despite of having `parse_eqs=true`, due to some internal error. + Using `parse_eqs = false`.""") + parse_eqs = false + end + end + + return parse_eqs, rv +end diff --git a/src/integrator/stepsize.jl b/src/integrator/stepsize.jl new file mode 100644 index 000000000..e9a290ee3 --- /dev/null +++ b/src/integrator/stepsize.jl @@ -0,0 +1,85 @@ +# This file is part of the TaylorIntegration.jl package; MIT licensed + +# stepsize +""" + stepsize(x, epsilon) -> h + +Returns a maximum time-step for a the Taylor expansion `x` +using a prescribed absolute tolerance `epsilon` and the last two +Taylor coefficients of (each component of) `x`. + +Note that `x` is of type `Taylor1{U}` or `Vector{Taylor1{U}}`, including +also the cases `Taylor1{TaylorN{U}}` and `Vector{Taylor1{TaylorN{U}}}`. + +Depending of `eltype(x)`, i.e., `U<:Number`, it may be necessary to overload +`stepsize`, specializing it on the type `U`, to avoid type instabilities. +""" +function stepsize(x::Taylor1{U}, epsilon::T) where {T<:Real, U<:Number} + R = promote_type(typeof(norm(constant_term(x), Inf)), T) + ord = x.order + h = convert(R, Inf) + z = zero(R) + for k in (ord-1, ord) + @inbounds aux = norm( x[k], Inf) + aux == z && continue + aux1 = _stepsize(aux, epsilon, k) + h = min(h, aux1) + end + return h::R +end + +function stepsize(q::AbstractArray{Taylor1{U}, N}, epsilon::T) where + {T<:Real, U<:Number, N} + R = promote_type(typeof(norm(constant_term(q[1]), Inf)), T) + h = convert(R, Inf) + for i in eachindex(q) + @inbounds hi = stepsize( q[i], epsilon ) + h = min( h, hi ) + end + + # If `isinf(h)==true`, we use the maximum (finite) + # step-size obtained from all coefficients as above. + # Note that the time step is independent from `epsilon`. + if isinf(h) + h = zero(R) + for i in eachindex(q) + @inbounds hi = _second_stepsize(q[i], epsilon) + h = max( h, hi ) + end + end + return h::R +end + +""" + _stepsize(aux1, epsilon, k) + +Helper function to avoid code repetition. +Returns `(epsilon/aux1)^(1/k)`. +""" +@inline function _stepsize(aux1::U, epsilon::T, k::Int) where {T<:Real, U<:Number} + aux = epsilon / aux1 + kinv = 1 / k + return aux^kinv +end + +""" + _second_stepsize(x, epsilon) + +Corresponds to the "second stepsize control" in Jorba and Zou +(2005) paper. We use it if [`stepsize`](@ref) returns `Inf`. +""" +function _second_stepsize(x::Taylor1{U}, epsilon::T) where {T<:Real, U<:Number} + R = promote_type(typeof(norm(constant_term(x), Inf)), T) + x == zero(x) && return convert(R, Inf) + ord = x.order + z = zero(R) + u = one(R) + h = z + for k in 1:ord-2 + @inbounds aux = norm( x[k], Inf) + aux == z && continue + aux1 = _stepsize(aux, u, k) + h = max(h, aux1) + end + return h::R +end \ No newline at end of file diff --git a/src/integrator/taylorinteg.jl b/src/integrator/taylorinteg.jl new file mode 100644 index 000000000..c94905c9c --- /dev/null +++ b/src/integrator/taylorinteg.jl @@ -0,0 +1,490 @@ +# This file is part of the TaylorIntegration.jl package; MIT licensed + +# set_psol! +@doc doc""" + set_psol!(::Val{true}, psol::Array{Taylor1{U},1}, nsteps::Int, x::Taylor1{U}) where {U<:Number} + set_psol!(::Val{true}, psol::Array{Taylor1{U},2}, nsteps::Int, x::Vector{Taylor1{U}}) where {U<:Number} + set_psol!(::Val{false}, args...) + +Auxiliary function to save Taylor polynomials in a call to [`taylorinteg`](@ref). When the +first argument in the call signature is `Val(true)`, sets appropriate elements of argument +`psol`. Otherwise, when the first argument in the call signature is `Val(false)`, this +function simply returns `nothing`. Argument `psol` is the array +where the Taylor polynomials associated to the solution will be stored, corresponding to +field `:p` in [`TaylorSolution`](@ref). See also [`init_psol`](@ref). +""" +@inline function set_psol!(::Val{true}, psol::Array{Taylor1{U},1}, nsteps::Int, x::Taylor1{U}) where {U<:Number} + @inbounds for k in eachindex(x) + TaylorSeries.identity!(psol[nsteps], x, k) + end + return nothing +end +@inline function set_psol!(::Val{true}, psol::Array{Taylor1{U},2}, nsteps::Int, x::Vector{Taylor1{U}}) where {U<:Number} + @inbounds for i in eachindex(x) + for k in eachindex(x[i]) + TaylorSeries.identity!(psol[i, nsteps], x[i], k) + end + end + return nothing +end +@inline set_psol!(::Val{false}, args...) = nothing + +@doc doc""" + init_psol(::Val{true}, xv::Array{U,1}, x::Taylor1{U}) where {U<:Number} + init_psol(::Val{true}, xv::Array{U,2}, x::Array{Taylor1{U},1}) where {U<:Number} + init_psol(::Val{false}, ::Array{U,1}, ::Taylor1{U}) where {U<:Number} + init_psol(::Val{false}, ::Array{U,2}, ::Array{Taylor1{U},1}) where {U<:Number} + +Auxiliary function to initialize `psol` during a call to [`taylorinteg`](@ref). When the +first argument in the call signature is `Val(false)` this function simply returns `nothing`. +Otherwise, when the first argument in the call signature is `Val(true)`, then the +appropriate array is allocated and returned; this array is where the Taylor polynomials +associated to the solution will be stored, corresponding to field `:p` in +[`TaylorSolution`](@ref). +""" +@inline function init_psol(::Val{true}, xv::Array{U,1}, x::Taylor1{U}) where {U<:Number} + return [zero(x) for _ in axes(xv, 1)] +end +@inline function init_psol(::Val{true}, xv::Array{U,2}, x::Array{Taylor1{U},1}) where {U<:Number} + return [zero(x[1]) for _ in axes(xv, 1), _ in axes(xv, 2)] +end + +# init_psol +@inline init_psol(::Val{false}, ::Array{U,1}, ::Taylor1{U}) where {U<:Number} = nothing +@inline init_psol(::Val{false}, ::Array{U,2}, ::Array{Taylor1{U},1}) where {U<:Number} = nothing + +# taylorinteg +function taylorinteg(f, x0::U, t0::T, tmax::T, order::Int, abstol::T, params = nothing; + maxsteps::Int=500, parse_eqs::Bool=true, dense::Bool=true) where {T<:Real, U<:Number} + + # Initialize the Taylor1 expansions + t = t0 + Taylor1( T, order ) + x = Taylor1( x0, order ) + + # Determine if specialized jetcoeffs! method exists + parse_eqs, rv = _determine_parsing!(parse_eqs, f, t, x, params) + + # Re-initialize the Taylor1 expansions + t = t0 + Taylor1( T, order ) + x = Taylor1( x0, order ) + return _taylorinteg!(Val(dense), f, t, x, x0, t0, tmax, abstol, rv, params; parse_eqs, maxsteps) +end + +function _taylorinteg!(dense::Val{D}, f, t::Taylor1{T}, x::Taylor1{U}, + x0::U, t0::T, tmax::T, abstol::T, rv::RetAlloc{Taylor1{U}}, params; + parse_eqs::Bool=true, maxsteps::Int=500) where {T<:Real, U<:Number, D} + + # Allocation + tv = Array{T}(undef, maxsteps+1) + xv = Array{U}(undef, maxsteps+1) + psol = init_psol(dense, xv, x) + + # Initial conditions + nsteps = 1 + @inbounds t[0] = t0 + @inbounds tv[1] = t0 + @inbounds xv[1] = x0 + sign_tstep = copysign(1, tmax-t0) + + # Integration + while sign_tstep*t0 < sign_tstep*tmax + δt = taylorstep!(Val(parse_eqs), f, t, x, abstol, params, rv) # δ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 + set_psol!(dense, psol, nsteps, x) # Store the Taylor polynomial solution + @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 build_solution(tv, xv, psol, 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, dense::Bool=true) where {T<:Real, U<:Number} + + # Initialize the vector of Taylor1 expansions + dof = length(q0) + t = t0 + Taylor1( T, order ) + x = Array{Taylor1{U}}(undef, dof) + dx = Array{Taylor1{U}}(undef, dof) + @inbounds for i in eachindex(q0) + x[i] = Taylor1( q0[i], order ) + dx[i] = Taylor1( zero(q0[i]), order ) + end + + # Determine if specialized jetcoeffs! method exists + parse_eqs, rv = _determine_parsing!(parse_eqs, f!, t, x, dx, params) + + # Re-initialize the Taylor1 expansions + t = t0 + Taylor1( T, order ) + x .= Taylor1.( q0, order ) + dx .= Taylor1.( zero.(q0), order) + return _taylorinteg!(Val(dense), f!, t, x, dx, q0, t0, tmax, abstol, rv, + params; parse_eqs, maxsteps) +end + +function _taylorinteg!(dense::Val{D}, f!, t::Taylor1{T}, x::Array{Taylor1{U},1}, dx::Array{Taylor1{U},1}, + q0::Array{U,1}, t0::T, tmax::T, abstol::T, rv::RetAlloc{Taylor1{U}}, params; + parse_eqs::Bool=true, maxsteps::Int=500) where {T<:Real, U<:Number, D} + + # Initialize the vector of Taylor1 expansions + dof = length(q0) + + # Allocation of output + tv = Array{T}(undef, maxsteps+1) + xv = Array{U}(undef, dof, maxsteps+1) + psol = init_psol(dense, xv, x) + xaux = Array{Taylor1{U}}(undef, dof) + + # Initial conditions + @inbounds t[0] = t0 + x0 = deepcopy(q0) + @inbounds tv[1] = t0 + @inbounds xv[:,1] .= q0 + sign_tstep = copysign(1, tmax-t0) + + # Integration + nsteps = 1 + while sign_tstep*t0 < sign_tstep*tmax + δt = taylorstep!(Val(parse_eqs), f!, t, x, dx, xaux, abstol, params, rv) # δ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 + set_psol!(dense, psol, nsteps, x) # Store the Taylor polynomial solution + @inbounds for i in eachindex(x0) + x[i][0] = x0[i] + TaylorSeries.zero!(dx[i], 0) + end + t0 += δt + @inbounds t[0] = t0 + nsteps += 1 + @inbounds tv[nsteps] = t0 + @inbounds xv[:,nsteps] .= deepcopy.(x0) + if nsteps > maxsteps + @warn(""" + Maximum number of integration steps reached; exiting. + """) + break + end + end + + return build_solution(tv, xv, psol, nsteps) +end + +@doc doc""" + taylorinteg(f, x0, t0, tmax, order, abstol, params[=nothing]; kwargs... ) + +General-purpose Taylor integrator for the explicit ODE ``\dot{x}=f(x, p, t)``, +where `p` are the parameters encoded in `params`. +The initial conditions are specified by `x0` at time `t0`; `x0` may be of type `T<:Number` +or `Vector{T}`, with `T` including `TaylorN{T}`; the latter case +is of interest for [jet transport applications](@ref jettransport). + +The equations of motion are specified by the function `f`; we follow the same +convention of `DifferentialEquations.jl` to define this function, i.e., +`f(x, p, t)` or `f!(dx, x, p, t)`; see the examples below. + +The functions returns a `TaylorSolution`, whose fields are `t` and `x`; they represent, +respectively, a vector with the values of time (independent variable), +and a vector with the computed values of +the dependent variable(s). When the keyword argument `dense` is set to `true`, it also +outputs in the field `p` the Taylor polynomial expansion computed at each time step. +The integration stops when time is larger than `tmax`, in which case the last returned +value(s) correspond to `tmax`, or when the number of saved steps is larger +than `maxsteps`. + +The integration method uses polynomial expansions on the independent variable +of order `order`; the parameter `abstol` serves to define the +time step using the last two Taylor coefficients of the expansions. +Make sure you use a *large enough* `order` to assure convergence. + +Currently, the recognized keyword arguments are: +- `maxsteps[=500]`: maximum number of integration steps. +- `parse_eqs[=true]`: use the specialized method of `jetcoeffs!` created + with [`@taylorize`](@ref). +- `dense[=true]`: output the Taylor polynomial expansion at each time step. + +## Examples + +For one dependent variable the function `f` defines the RHS of the equation of +motion, returning the value of ``\dot{x}``. The arguments of +this function are `(x, p, t)`, where `x` are the dependent variables, `p` are +the paremeters and `t` is the independent variable. + +For several (two or more) dependent variables, the function `f!` defines +the RHS of the equations of motion, mutating (in-place) the (preallocated) vector +with components of ``\dot{x}``. The arguments of this function are `(dx, x, p, t)`, +where `dx` is the preallocated vector of ``\dot{x}``, `x` are the dependent +variables, `p` are the paremeters entering the ODEs and `t` is the independent +variable. The function may return this vector or simply `nothing`. + +```julia +using TaylorIntegration + +f(x, p, t) = x^2 + +sol = taylorinteg(f, 3, 0.0, 0.3, 25, 1.0e-20, maxsteps=100 ) + +function f!(dx, x, p, t) + for i in eachindex(x) + dx[i] = x[i]^2 + end + return nothing +end + +sol = taylorinteg(f!, [3, 3], 0.0, 0.3, 25, 1.0e-20, maxsteps=100 ) + +sol = taylorinteg(f!, [3, 3], 0.0, 0.3, 25, 1.0e-20, maxsteps=100, dense=true ) +``` + +""" taylorinteg + + +# Integrate and return results evaluated at given time +@doc doc""" + taylorinteg(f, x0, trange, order, abstol, params[=nothing]; kwargs... ) + +General-purpose Taylor integrator for the explicit ODE +``\dot{x}=f(x,p,t)`` with initial condition specified by `x0::{T<:Number}` +or `x0::Vector{T}` at time `t0`. + +The method returns a `TaylorSolution` whose field `x` represents the computed values of +the dependent variable(s), evaluated *only* at the times specified by +the range `trange`. + +## Examples + +```julia +sol = taylorinteg(f, 3, 0.0:0.001:0.3, 25, 1.0e-20, maxsteps=100 ) + +sol = taylorinteg(f!, [3, 3], 0.0:0.001:0.3, 25, 1.0e-20, maxsteps=100 ); + +``` + +""" +function taylorinteg(f, x0::U, trange::AbstractVector{T}, + order::Int, abstol::T, params = nothing; + maxsteps::Int=500, parse_eqs::Bool=true) where {T<:Real, U<:Number} + + # Check if trange is increasingly or decreasingly sorted + @assert (issorted(trange) || + issorted(reverse(trange))) "`trange` or `reverse(trange)` must be sorted" + + # Initialize the Taylor1 expansions + t0 = trange[1] + t = t0 + Taylor1( T, order ) + x = Taylor1( x0, order ) + + # Determine if specialized jetcoeffs! method exists + parse_eqs, rv = _determine_parsing!(parse_eqs, f, t, x, params) + + # Re-initialize the Taylor1 expansions + t = t0 + Taylor1( T, order ) + x = Taylor1( x0, order ) + return _taylorinteg!(f, t, x, x0, trange, abstol, rv, params; parse_eqs, maxsteps) +end + +function _taylorinteg!(f, t::Taylor1{T}, x::Taylor1{U}, x0::U, trange::AbstractVector{T}, + abstol::T, rv::RetAlloc{Taylor1{U}}, params; parse_eqs::Bool=true, maxsteps::Int=500) where {T<:Real, U<:Number} + + # Allocation + nn = length(trange) + xv = Array{U}(undef, nn) + fill!(xv, T(NaN)) + + # Initial conditions + @inbounds t0, t1, tmax = trange[1], trange[2], trange[end] + sign_tstep = copysign(1, tmax-t0) + @inbounds t[0] = t0 + @inbounds xv[1] = x0 + + # Integration + iter = 2 + nsteps = 1 + while sign_tstep*t0 < sign_tstep*tmax + δt = taylorstep!(Val(parse_eqs), f, t, x, abstol, params, rv)# δ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 + tnext = t0+δt + # Evaluate solution at times within convergence radius + while sign_tstep*t1 < sign_tstep*tnext + x1 = evaluate(x, t1-t0) + @inbounds xv[iter] = x1 + iter += 1 + @inbounds t1 = trange[iter] + end + if δt == tmax-t0 + @inbounds xv[iter] = x0 + break + end + @inbounds x[0] = x0 + t0 = tnext + @inbounds t[0] = t0 + nsteps += 1 + if nsteps > maxsteps + @warn(""" + Maximum number of integration steps reached; exiting. + """) + break + end + end + return build_solution(trange, xv) +end + +function taylorinteg(f!, q0::Array{U,1}, trange::AbstractVector{T}, + order::Int, abstol::T, params = nothing; + maxsteps::Int=500, parse_eqs::Bool=true) where {T<:Real, U<:Number} + + # Check if trange is increasingly or decreasingly sorted + @assert (issorted(trange) || + issorted(reverse(trange))) "`trange` or `reverse(trange)` must be sorted" + + # Initialize the vector of Taylor1 expansions + dof = length(q0) + t0 = trange[1] + t = t0 + Taylor1( T, order ) + x = Array{Taylor1{U}}(undef, dof) + dx = Array{Taylor1{U}}(undef, dof) + x .= Taylor1.( q0, order ) + dx .= Taylor1.( zero.(q0), order ) + + # Determine if specialized jetcoeffs! method exists + parse_eqs, rv = _determine_parsing!(parse_eqs, f!, t, x, dx, params) + + # Re-initialize the Taylor1 expansions + t = t0 + Taylor1( T, order ) + x .= Taylor1.( q0, order ) + dx .= Taylor1.( zero.(q0), order ) + return _taylorinteg!(f!, t, x, dx, q0, trange, abstol, rv, + params; parse_eqs, maxsteps) +end + +function _taylorinteg!(f!, t::Taylor1{T}, x::Array{Taylor1{U},1}, dx::Array{Taylor1{U},1}, + q0::Array{U,1}, trange::AbstractVector{T}, abstol::T, rv::RetAlloc{Taylor1{U}}, params; + parse_eqs::Bool=true, maxsteps::Int=500) where {T<:Real, U<:Number} + + # Allocation + nn = length(trange) + dof = length(q0) + x0 = similar(q0, eltype(q0), dof) + x1 = similar(x0) + fill!(x0, T(NaN)) + xv = Array{eltype(q0)}(undef, dof, nn) + for ind in 1:nn + @inbounds xv[:,ind] .= x0 + end + xaux = Array{Taylor1{U}}(undef, dof) + + # Initial conditions + @inbounds t[0] = trange[1] + @inbounds t0, t1, tmax = trange[1], trange[2], trange[end] + sign_tstep = copysign(1, tmax-t0) + # x .= Taylor1.(q0, order) + @inbounds x0 .= q0 + @inbounds xv[:,1] .= q0 + + # Integration + iter = 2 + nsteps = 1 + while sign_tstep*t0 < sign_tstep*tmax + δt = taylorstep!(Val(parse_eqs), f!, t, x, dx, xaux, abstol, params, rv) # δ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 + tnext = t0+δt + # Evaluate solution at times within convergence radius + while sign_tstep*t1 < sign_tstep*tnext + evaluate!(x, t1-t0, x1) + @inbounds xv[:,iter] .= x1 + iter += 1 + @inbounds t1 = trange[iter] + end + if δt == tmax-t0 + @inbounds xv[:,iter] .= x0 + break + end + @inbounds for i in eachindex(x0) + x[i][0] = x0[i] + dx[i][0] = zero(x0[i]) + end + t0 = tnext + @inbounds t[0] = t0 + nsteps += 1 + if nsteps > maxsteps + @warn(""" + Maximum number of integration steps reached; exiting. + """) + break + end + end + + return build_solution(trange, xv) +end + + +# Generic functions +for R in (:Number, :Integer) + @eval begin + + function taylorinteg(f, xx0::S, tt0::T, ttmax::U, order::Int, aabstol::V, + params = nothing; dense=false, 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, _ = promote(tt0, ttmax, aabstol, one(Float64)) + x0, _ = promote(xx0, t0) + + return taylorinteg(f, x0, t0, tmax, order, abstol, params, + dense=dense, maxsteps=maxsteps, parse_eqs=parse_eqs) + end + + function taylorinteg(f, q0::Array{S,1}, tt0::T, ttmax::U, order::Int, aabstol::V, + params = nothing; dense=false, 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, _ = promote(tt0, ttmax, aabstol, one(Float64)) + elq0, _ = 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, params, + dense=dense, maxsteps=maxsteps, parse_eqs=parse_eqs) + end + + function taylorinteg(f, xx0::S, trange::AbstractVector{T}, order::Int, aabstol::U, params = nothing; + maxsteps::Int=500, parse_eqs::Bool=true) where {S<:$R, T<:Real, U<:Real} + + t0, abstol, _ = promote(trange[1], aabstol, one(Float64)) + x0, _ = promote(xx0, t0) + + return taylorinteg(f, x0, trange .* one(t0), order, abstol, params, + maxsteps=maxsteps, parse_eqs=parse_eqs) + end + + function taylorinteg(f, q0::Array{S,1}, trange::AbstractVector{T}, order::Int, aabstol::U, params = nothing; + maxsteps::Int=500, parse_eqs::Bool=true) where {S<:$R, T<:Real, U<:Real} + + t0, abstol, _ = promote(trange[1], aabstol, one(Float64)) + elq0, _ = promote(q0[1], t0) + q0_ = convert(Array{typeof(elq0)}, q0) + + return taylorinteg(f, q0_, trange .* one(t0), order, abstol, params, + maxsteps=maxsteps, parse_eqs=parse_eqs) + end + + end +end diff --git a/src/integrator/taylorsolution.jl b/src/integrator/taylorsolution.jl new file mode 100644 index 000000000..5cd260b49 --- /dev/null +++ b/src/integrator/taylorsolution.jl @@ -0,0 +1,171 @@ +# This file is part of the TaylorIntegration.jl package; MIT licensed + +abstract type AbstractTaylorSolution{T<:Real, U<:Number} end + +## Constructors + +""" +TaylorSolution{T, U, N, VT<:AbstractVector{T}, AX<:AbstractArray{U,N}, + P<:Union{Nothing, AbstractArray{Taylor1{U}, N}}, + VTE<:Union{Nothing, AbstractVector{U}}, + AXE<:Union{Nothing, AbstractArray{U, N}}, + VΛ<:Union{Nothing, AbstractArray{U,N}}} <: AbstractTaylorSolution{T, U} + +This `struct` represents the return type for `taylorinteg`. Fields `t` and `x` represent, +respectively, a vector with the values of time (independent variable), and a vector with the +computed values of the dependent variable(s). When `taylorinteg` is called with `dense=true`, +then field `p` stores the Taylor polynomial expansion computed at each time step. Fields +`tevents`, `xevents` and `gresids` are related to root-finding methods of `taylorinteg`, while +`λ` is related to the output of [lyap_taylorinteg](@ref). +""" +struct TaylorSolution{T, U, N, VT<:AbstractVector{T}, AX<:AbstractArray{U,N}, + P<:Union{Nothing, AbstractArray{Taylor1{U}, N}}, + VTE<:Union{Nothing, AbstractVector{U}}, + AXE<:Union{Nothing, AbstractArray{U, N}}, + VΛ<:Union{Nothing, AbstractArray{U,N}}} <: AbstractTaylorSolution{T, U} + t::VT + x::AX + p::P + tevents::VTE + xevents::AXE + gresids::VTE + λ::VΛ + function TaylorSolution{T, U, N, VT, AX, P, VTE, AXE, VΛ}(t::VT, x::AX, p::P, tevents::VTE, xevents::AXE, gresids::VTE, λ::VΛ) where {T, U, + N, VT, AX, P, VTE, AXE, VΛ} + @assert length(t) == size(x, 1) + @assert issorted(t) || issorted(t, rev = true) + !isnothing(p) && begin + @assert size(x, 1) - 1 == size(p, 1) + @assert size(x)[2:end] == size(p)[2:end] + end + @assert isnothing(tevents) == isnothing(xevents) == isnothing(gresids) "`Nothing`-ness must be consistent across `tevents`, `xevents` and `gresids`." + !isnothing(tevents) && begin + @assert length(tevents) == size(xevents, 1) + @assert size(xevents, 2) == size(x, 2) + @assert length(tevents) == length(gresids) + end + !isnothing(λ) && @assert size(λ) == size(x) + return new{T, U, N, VT, AX, P, VTE, AXE, VΛ}(t, x, p, tevents, xevents, gresids, λ) + end +end +TaylorSolution(t::VT, x::AX, p::P, tevents::VTE, xevents::AXE, gresids::VTE, λ::VΛ) where { + T, U, N, VT<:AbstractVector{T}, AX<:AbstractArray{U,N}, + P<:Union{Nothing, AbstractArray{Taylor1{U},N}}, + VTE<:Union{Nothing, AbstractVector{U}}, AXE<:Union{Nothing, AbstractArray{U, N}}, + VΛ<:Union{Nothing, AbstractArray{U,N}}} = + TaylorSolution{T, U, N, VT, AX, P, VTE, AXE, VΛ}(t, x, p, tevents, xevents, gresids, λ) + +# 4-arg constructor (root-finding and Lyapunov fields are nothing) +TaylorSolution(t, x, p, ::Nothing) = TaylorSolution(t, x, p, nothing, nothing, nothing, nothing) +# 3-arg constructor (root-finding and Lyapunov fields are nothing; helps not to write too many nothings) +TaylorSolution(t, x, p) = TaylorSolution(t, x, p, nothing) +# 2-arg constructor (dense polynomial, root-finding and Lyapunov fields are nothing) +TaylorSolution(t, x) = TaylorSolution(t, x, nothing) + +# arraysol: auxiliary function for solution construction + +arraysol(::Nothing, ::Int) = nothing +arraysol(v::AbstractVector, n::Int) = view(v, 1:n) +arraysol(m::AbstractMatrix, n::Int) = view( transpose(view(m, :, 1:n)), 1:n, : ) + +# build_solution + +""" + build_solution(t, x, p, nsteps) + build_solution(t, x) + +Helper function to build a [`TaylorSolution`](@ref) from a call to +[`taylorinteg`](@ref). + +""" +build_solution(t::AbstractVector{T}, x::Vector{U}, ::Nothing, nsteps::Int) where {T, U} = + TaylorSolution(arraysol(t, nsteps), arraysol(x, nsteps)) +build_solution(t::AbstractVector{T}, x::Vector{U}, p::Vector{Taylor1{U}}, nsteps::Int) where {T, U} = + TaylorSolution(arraysol(t, nsteps), arraysol(x, nsteps), arraysol(p, nsteps-1)) +build_solution(t::AbstractVector{T}, x::Matrix{U}, ::Nothing, nsteps::Int) where {T, U} = + TaylorSolution(arraysol(t, nsteps), arraysol(x, nsteps)) +build_solution(t::AbstractVector{T}, x::Matrix{U}, p::Matrix{Taylor1{U}}, nsteps::Int) where {T, U} = + TaylorSolution(arraysol(t, nsteps), arraysol(x, nsteps), arraysol(p, nsteps-1)) + +build_solution(t::AbstractVector{T}, x::Vector{U}) where {T, U} = TaylorSolution(t, x) +build_solution(t::AbstractVector{T}, x::Matrix{U}) where {T, U} = TaylorSolution(t, transpose(x)) + +### Custom print + +function Base.show(io::IO, sol::TaylorSolution) + tspan = minmax(sol.t[1], sol.t[end]) + S = eltype(sol.x) + nvars = size(sol.x, 2) + plural = nvars > 1 ? "s" : "" + print(io, "tspan: ", tspan, ", x: ", nvars, " ", S, " variable"*plural) +end + +# Function-like (callability or "functor") methods + +@doc doc""" + timeindex(sol::TaylorSolution, t::TT) where TT + +Return the index of `sol.t` corresponding to `t` and the time elapsed from `sol.t0` +to `t`. +""" +function timeindex(sol::TaylorSolution, t::TT) where TT + t0 = sol.t[1] + _t = constant_term(constant_term(t)) # Current time + tmin, tmax = minmax(sol.t[end], t0) # Min and max time in sol + + @assert tmin ≤ _t ≤ tmax "Evaluation time outside range of interpolation" + + if _t == sol.t[end] # Compute solution at final time from last step expansion + ind = lastindex(sol.t) - 1 + elseif issorted(sol.t) # Forward integration + ind = searchsortedlast(sol.t, _t) + elseif issorted(sol.t, rev=true) # Backward integration + ind = searchsortedlast(sol.t, _t, rev=true) + end + # Time since the start of the ind-th timestep + δt = t - sol.t[ind] + # Return index and elapsed time since i-th timestep + return ind::Int, δt::TT +end + +const TaylorSolutionCallingArgs{T,U} = Union{T, U, Taylor1{U}, TaylorN{U}, Taylor1{TaylorN{U}}} where {T,U} + +@doc doc""" + (sol::TaylorSolution{T, U, 1})(t::T) where {T, U} + (sol::TaylorSolution{T, U, 1})(t::TT) where {T, U, TT<:TaylorSolutionCallingArgs{T,U}} + (sol::TaylorSolution{T, U, 2})(t::T) where {T, U} + (sol::TaylorSolution{T, U, 2})(t::TT) where {T, U, TT<:TaylorSolutionCallingArgs{T,U}} + +Evaluate `sol.x` at time `t`. + +See also [`timeindex`](@ref). +""" +(sol::TaylorSolution)(t) = sol(Val(isnothing(sol.p)), t) + +(sol::TaylorSolution{T, U, N})(::Val{true}, t::TT) where {T, U, N, TT<:TaylorSolutionCallingArgs{T,U}} = error("`TaylorSolution` objects computed from calls to `taylorinteg` with `dense=false` are not callable.") + +function (sol::TaylorSolution{T, U, 1})(::Val{false}, t::T) where {T, U} + # Get index of sol.x that interpolates at time t + ind::Int, δt::T = timeindex(sol, t) + # Evaluate sol.x[ind] at δt + return (sol.p[ind])(δt)::U +end +function (sol::TaylorSolution{T, U, 1})(::Val{false}, t::TT) where {T, U, TT<:TaylorSolutionCallingArgs{T,U}} + # Get index of sol.x that interpolates at time t + ind::Int, δt::TT = timeindex(sol, t) + # Evaluate sol.x[ind] at δt + return (sol.p[ind])(δt)::TT +end + +function (sol::TaylorSolution{T, U, 2})(::Val{false}, t::T) where {T, U} + # Get index of sol.x that interpolates at time t + ind::Int, δt::T = timeindex(sol, t) + # Evaluate sol.x[ind] at δt + return view(sol.p, ind, :)(δt)::Vector{U} +end +function (sol::TaylorSolution{T, U, 2})(::Val{false}, t::TT) where {T, U, TT<:TaylorSolutionCallingArgs{T,U}} + # Get index of sol.x that interpolates at time t + ind::Int, δt::TT = timeindex(sol, t) + # Evaluate sol.x[ind] at δt + return view(sol.p, ind, :)(δt)::Vector{TT} +end \ No newline at end of file diff --git a/src/integrator/taylorstep.jl b/src/integrator/taylorstep.jl new file mode 100644 index 000000000..7c0d60fe5 --- /dev/null +++ b/src/integrator/taylorstep.jl @@ -0,0 +1,59 @@ +# This file is part of the TaylorIntegration.jl package; MIT licensed + +#taylorstep +@doc doc""" + taylorstep!(::Val{false}, f, t, x, abstol, params, rv, parse_eqs=true) -> δt + taylorstep!(::Val{true}, f, t, x, abstol, params, rv, parse_eqs=true) -> δt + taylorstep!(::Val{false}, f!, t, x, dx, xaux, abstol, params, rv, parse_eqs=true) -> δt + taylorstep!(::Val{true}, f!, t, x, dx, xaux, abstol, params, rv, parse_eqs=true) -> δt + +One-step Taylor integration for the one-dependent variable ODE ``\dot{x}=dx/dt=f(x, p, t)`` +with initial conditions ``x(t_0)=x_0``. +Returns the time-step `δt` of the actual integration carried out (δt is positive). + +Here, `f` represents the function defining the RHS of the ODE (see +[`taylorinteg`](@ref)), `t` is the +independent variable, `x` contains the Taylor expansion of the dependent +variable, `order` is +the degree used for the `Taylor1` polynomials during the integration +`abstol` is the absolute tolerance used to determine the time step +of the integration, and `params` are the parameters entering the ODE +functions. +For several variables, `dx` and `xaux`, both of the same type as `x`, +are needed to save allocations. Finally, `parse_eqs` is a switch +to force *not* using (`parse_eqs=false`) the specialized method of `jetcoeffs!` +created with [`@taylorize`](@ref); the default is `true` (parse the equations). +Finally, `parse_eqs` is a switch +to force *not* using (`parse_eqs=false`) the specialized method of `jetcoeffs!` +created with [`@taylorize`](@ref); the default is `true` (parse the equations). +The first argument in the function call `Val{bool}` (`bool::Bool`) controls whether +a specialized [`jetcoeffs!](@ref) method is being used or not. + +""" +function taylorstep!(::Val{V}, f, t::Taylor1{T}, x::Taylor1{U}, abstol::T, params, + rv::RetAlloc{Taylor1{U}}) where {T<:Real, U<:Number, V} + + # Compute the Taylor coefficients + __jetcoeffs!(Val(V), f, t, x, params, rv) + + # Compute the step-size of the integration using `abstol` + δt = stepsize(x, abstol) + if isinf(δt) + δt = _second_stepsize(x, abstol) + end + + return δt +end + +function taylorstep!(::Val{V}, f!, t::Taylor1{T}, x::Vector{Taylor1{U}}, + dx::Vector{Taylor1{U}}, xaux::Vector{Taylor1{U}}, abstol::T, params, + rv::RetAlloc{Taylor1{U}}) where {T<:Real, U<:Number, V} + + # Compute the Taylor coefficients + __jetcoeffs!(Val(V), f!, t, x, dx, xaux, params, rv) + + # Compute the step-size of the integration using `abstol` + δt = stepsize(x, abstol) + + return δt +end \ No newline at end of file diff --git a/src/lyapunovspectrum.jl b/src/lyapunovspectrum.jl index bd50f4307..c3153df11 100644 --- a/src/lyapunovspectrum.jl +++ b/src/lyapunovspectrum.jl @@ -1,5 +1,40 @@ # This file is part of the TaylorIntegration.jl package; MIT licensed +# build_lyap_solution + +@doc doc""" + build_lyap_solution(t, x, λ, nsteps) + build_lyap_solution(t, x, λ) + +Helper function to build a [`TaylorSolution`](@ref) from a call to +[`lyap_taylorinteg`](@ref). + +""" +build_lyap_solution(t::AbstractVector{T}, + x::Matrix{U}, + λ::Matrix{U}, + nsteps::Int) where {T, U} = + TaylorSolution(arraysol(t, nsteps), + arraysol(x, nsteps), + nothing, + nothing, + nothing, + nothing, + arraysol(λ, nsteps)) + +build_lyap_solution(t::AbstractVector{T}, + x::Matrix{U}, + λ::Matrix{U}) where {T, U} = + TaylorSolution(t, + transpose(x), + nothing, + nothing, + nothing, + nothing, + transpose(λ)) + +# stabilitymatrix! + """ stabilitymatrix!(eqsdiff!, t, x, δx, dδx, jac, _δv, params[, jacobianfunc!=nothing]) @@ -15,21 +50,25 @@ Optionally, the user may provide a Jacobian function `jacobianfunc!` to compute function stabilitymatrix!(eqsdiff!, t::Taylor1{T}, x::Vector{Taylor1{U}}, δx::Vector{TaylorN{Taylor1{U}}}, dδx::Vector{TaylorN{Taylor1{U}}}, jac::Matrix{Taylor1{U}}, _δv::Vector{TaylorN{Taylor1{U}}}, params, - jacobianfunc! = nothing) where {T<:Real, U<:Number} - - if isa(jacobianfunc!, Nothing) - # Set δx equal to current value of x plus 1st-order variations - for ind in eachindex(δx) - @inbounds δx[ind] = x[ind] + _δv[ind] - end - # Equations of motion - eqsdiff!(dδx, δx, params, t) - TaylorSeries.jacobian!(jac, dδx) - else - jacobianfunc!(jac, x, params, t) + ::Nothing) where {T<:Real, U<:Number} + # Set δx equal to current value of x plus 1st-order variations + for ind in eachindex(δx) + @inbounds δx[ind] = x[ind] + _δv[ind] end + # Equations of motion + eqsdiff!(dδx, δx, params, t) + TaylorSeries.jacobian!(jac, dδx) + nothing +end +function stabilitymatrix!(eqsdiff!, t::Taylor1{T}, x::Vector{Taylor1{U}}, + δx::Vector{TaylorN{Taylor1{U}}}, dδx::Vector{TaylorN{Taylor1{U}}}, + jac::Matrix{Taylor1{U}}, _δv::Vector{TaylorN{Taylor1{U}}}, params, + jacobianfunc!::J) where {T<:Real, U<:Number, J} + jacobianfunc!(jac, x, params, t) nothing end +stabilitymatrix!(eqsdiff!, t, x, δx, dδx, jac, _δv, params) = + stabilitymatrix!(eqsdiff!, t, x, δx, dδx, jac, _δv, params, nothing) # Modified from `cgs` and `mgs`, obtained from: # http://nbviewer.jupyter.org/url/math.mit.edu/~stevenj/18.335/Gram-Schmidt.ipynb @@ -155,7 +194,7 @@ function lyap_jetcoeffs!(t::Taylor1{T}, x::AbstractVector{Taylor1{S}}, end """ - lyap_taylorstep!(f!, t, x, dx, xaux, δx, dδx, jac, t0, t1, order, abstol, _δv, varsaux, params[, jacobianfunc!]) + lyap_taylorstep!(::Val{V}, f!, t, x, dx, xaux, δx, dδx, jac, abstol, _δv, varsaux, params, rv[, jacobianfunc!]) Similar to [`taylorstep!`](@ref) for the calculation of the Lyapunov spectrum. `jac` is the Taylor expansion (wrt the independent variable) of the @@ -166,42 +205,19 @@ Jacobian function `jacobianfunc!` to compute `jac`. Otherwise, `jac` is computed via automatic differentiation using `TaylorSeries.jl`. """ -function lyap_taylorstep!(f!, t::Taylor1{T}, x::Vector{Taylor1{U}}, - dx::Vector{Taylor1{U}}, xaux::Vector{Taylor1{U}}, - δx::Array{TaylorN{Taylor1{U}},1}, dδx::Array{TaylorN{Taylor1{U}},1}, - jac::Array{Taylor1{U},2}, abstol::T, _δv::Vector{TaylorN{Taylor1{U}}}, - varsaux::Array{Taylor1{U},3}, params, jacobianfunc! =nothing) where {T<:Real, U<:Number} - - # Dimensions of phase-space: dof - nx = length(x) - dof = length(δx) - - # Compute the Taylor coefficients associated to trajectory - __jetcoeffs!(Val(false), f!, t, view(x, 1:dof), view(dx, 1:dof), view(xaux, 1:dof), params) - - # Compute stability matrix - stabilitymatrix!(f!, t, x, δx, dδx, jac, _δv, params, jacobianfunc!) - - # Compute the Taylor coefficients associated to variational equations - lyap_jetcoeffs!(t, view(x, dof+1:nx), view(dx, dof+1:nx), jac, varsaux) - - # Compute the step-size of the integration using `abstol` - δt = stepsize(view(x, 1:dof), abstol) - return δt -end -function lyap_taylorstep!(f!, t::Taylor1{T}, x::Vector{Taylor1{U}}, dx::Vector{Taylor1{U}}, - δx::Array{TaylorN{Taylor1{U}},1}, dδx::Array{TaylorN{Taylor1{U}},1}, +function lyap_taylorstep!(::Val{V}, f!, t::Taylor1{T}, x::Vector{Taylor1{U}}, dx::Vector{Taylor1{U}}, + xaux::Vector{Taylor1{U}}, δx::Array{TaylorN{Taylor1{U}},1}, dδx::Array{TaylorN{Taylor1{U}},1}, jac::Array{Taylor1{U},2}, abstol::T, _δv::Vector{TaylorN{Taylor1{U}}}, varsaux::Array{Taylor1{U},3}, params, rv::RetAlloc{Taylor1{U}}, - jacobianfunc! =nothing) where {T<:Real, U<:Number} + jacobianfunc! =nothing) where {T<:Real, U<:Number, V} # Dimensions of phase-space: dof nx = length(x) dof = length(δx) # Compute the Taylor coefficients associated to trajectory - __jetcoeffs!(Val(true), f!, t, view(x, 1:dof), view(dx, 1:dof), params, rv) + __jetcoeffs!(Val(V), f!, t, view(x, 1:dof), view(dx, 1:dof), view(xaux, 1:dof), params, rv) # Compute stability matrix stabilitymatrix!(f!, t, x, δx, dδx, jac, _δv, params, jacobianfunc!) @@ -260,32 +276,26 @@ function lyap_taylorinteg(f!, q0::Array{U,1}, t0::T, tmax::T, parse_eqs, rv = _determine_parsing!(parse_eqs, f!, t, view(x, 1:dof), view(dx, 1:dof), params) - if parse_eqs - # Re-initialize the Taylor1 expansions - t = t0 + Taylor1( T, order ) - x .= Taylor1.( x0, order ) - return _lyap_taylorinteg!(f!, t, x, dx, q0, t0, tmax, abstol, jt, _dv, rv, - params, jacobianfunc!, maxsteps=maxsteps) - else - return _lyap_taylorinteg!(f!, t, x, dx, q0, t0, tmax, abstol, jt, _dv, - params, jacobianfunc!, maxsteps=maxsteps) - end + # Re-initialize the Taylor1 expansions + t = t0 + Taylor1( T, order ) + x .= Taylor1.( x0, order ) + return _lyap_taylorinteg!(f!, t, x, dx, q0, t0, tmax, abstol, jt, _dv, rv, + params, jacobianfunc!; maxsteps, parse_eqs) end function _lyap_taylorinteg!(f!, t::Taylor1{T}, x::Array{Taylor1{U},1}, dx::Array{Taylor1{U},1}, q0::Array{U,1}, t0::T, tmax::T, abstol::T, jt::Matrix{U}, _δv::Array{TaylorN{Taylor1{U}},1}, - params, jacobianfunc!; maxsteps::Int=500) where {T<:Real, U<:Number} + rv::RetAlloc{Taylor1{U}}, params, jacobianfunc!; parse_eqs::Bool=true, + maxsteps::Int=500) where {T<:Real, U<:Number} # Allocation order = get_order(t) tv = Array{T}(undef, maxsteps+1) dof = length(q0) - nx0 = length(x) x0 = getcoeff.(x, 0) xv = Array{U}(undef, dof, maxsteps+1) λ = similar(xv) λtsum = similar(q0) - # q1 = similar(q0) # Initial conditions @inbounds t[0] = t0 @@ -300,6 +310,7 @@ function _lyap_taylorinteg!(f!, t::Taylor1{T}, x::Array{Taylor1{U},1}, dx::Array end #Allocate auxiliary arrays + nx0 = length(x0) xaux = Array{Taylor1{U}}(undef, nx0) δx = Array{TaylorN{Taylor1{U}}}(undef, dof) dδx = Array{TaylorN{Taylor1{U}}}(undef, dof) @@ -315,81 +326,7 @@ function _lyap_taylorinteg!(f!, t::Taylor1{T}, x::Array{Taylor1{U},1}, dx::Array # Integration nsteps = 1 while sign_tstep*t0 < sign_tstep*tmax - δt = lyap_taylorstep!(f!, t, x, dx, xaux, δx, dδx, jac, abstol, _δv, - varsaux, params, jacobianfunc!) # δ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) # Update x0 - for ind in eachindex(jt) - @inbounds jt[ind] = x0[dof+ind] - end - modifiedGS!( jt, QH, RH, aⱼ, qᵢ, vⱼ ) - t0 += δt - @inbounds t[0] = t0 - tspan = t0-t00 - nsteps += 1 - @inbounds tv[nsteps] = t0 - @inbounds for ind in eachindex(q0) - xv[ind,nsteps] = x0[ind] - λtsum[ind] += log(RH[ind,ind]) - λ[ind,nsteps] = λtsum[ind]/tspan - end - for ind in eachindex(QH) - @inbounds x0[dof+ind] = QH[ind] - end - x .= Taylor1.( x0, order ) - if nsteps > maxsteps - @warn(""" - Maximum number of integration steps reached; exiting. - """) - break - end - end - - return view(tv,1:nsteps), view(transpose(xv),1:nsteps,:), view(transpose(λ),1:nsteps,:) -end - -function _lyap_taylorinteg!(f!, t::Taylor1{T}, x::Array{Taylor1{U},1}, dx::Array{Taylor1{U},1}, - q0::Array{U,1}, t0::T, tmax::T, abstol::T, jt::Matrix{U}, _δv::Array{TaylorN{Taylor1{U}},1}, - rv::RetAlloc{Taylor1{U}}, params, jacobianfunc!; maxsteps::Int=500) where {T<:Real, U<:Number} - - # Allocation - order = get_order(t) - tv = Array{T}(undef, maxsteps+1) - dof = length(q0) - x0 = getcoeff.(x, 0) - xv = Array{U}(undef, dof, maxsteps+1) - λ = similar(xv) - λtsum = similar(q0) - - # Initial conditions - @inbounds t[0] = t0 - sign_tstep = copysign(1, tmax-t0) - t00 = t0 - tspan = zero(T) - @inbounds tv[1] = t0 - @inbounds for ind in eachindex(q0) - xv[ind,1] = q0[ind] - λ[ind,1] = zero(U) - λtsum[ind] = zero(U) - end - - #Allocate auxiliary arrays - δx = Array{TaylorN{Taylor1{U}}}(undef, dof) - dδx = Array{TaylorN{Taylor1{U}}}(undef, dof) - jac = Array{Taylor1{U}}(undef, dof, dof) - varsaux = Array{Taylor1{U}}(undef, dof, dof, dof) - fill!(jac, zero(x[1])) - QH = Array{U}(undef, dof, dof) - RH = Array{U}(undef, dof, dof) - aⱼ = Array{U}(undef, dof ) - qᵢ = similar(aⱼ) - vⱼ = similar(aⱼ) - - # Integration - nsteps = 1 - while sign_tstep*t0 < sign_tstep*tmax - δt = lyap_taylorstep!(f!, t, x, dx, δx, dδx, jac, abstol, _δv, + δt = lyap_taylorstep!(Val(parse_eqs), f!, t, x, dx, xaux, δx, dδx, jac, abstol, _δv, varsaux, params, rv, jacobianfunc!) # δ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)) @@ -420,7 +357,7 @@ function _lyap_taylorinteg!(f!, t::Taylor1{T}, x::Array{Taylor1{U},1}, dx::Array end end - return view(tv,1:nsteps), view(transpose(xv),1:nsteps,:), view(transpose(λ),1:nsteps,:) + return build_lyap_solution(tv, xv, λ, nsteps) end @@ -445,7 +382,7 @@ function lyap_taylorinteg(f!, q0::Array{U,1}, trange::AbstractVector{T}, _dv = Array{TaylorN{Taylor1{U}}}(undef, dof) # If user does not provide Jacobian, check number of TaylorN variables and initialize _dv - if isa(jacobianfunc!, Nothing) + if isnothing(jacobianfunc!) @assert get_numvars() == dof "`length(q0)` must be equal to number of variables set by `TaylorN`" for ind in eachindex(q0) _dv[ind] = one(x[1])*TaylorN(Taylor1{U}, ind, order=1) @@ -456,129 +393,19 @@ function lyap_taylorinteg(f!, q0::Array{U,1}, trange::AbstractVector{T}, parse_eqs, rv = _determine_parsing!(parse_eqs, f!, t, view(x, 1:dof), view(dx, 1:dof), params) - if parse_eqs - # Re-initialize the Taylor1 expansions - t = trange[1] + Taylor1( T, order ) - x .= Taylor1.( x0, order ) - tmpTaylor, arrTaylor = rv.v0, rv.v1 - return _lyap_taylorinteg!(f!, t, x, dx, q0, trange, abstol, jt, _dv, rv, - params, jacobianfunc!, maxsteps=maxsteps) - else - return _lyap_taylorinteg!(f!, t, x, dx, q0, trange, abstol, jt, _dv, - params, jacobianfunc!, maxsteps=maxsteps) - end - -end - -function _lyap_taylorinteg!(f!, t::Taylor1{T}, x::Array{Taylor1{U},1}, dx::Array{Taylor1{U},1}, - q0::Array{U,1}, trange::AbstractVector{T}, abstol::T, jt::Matrix{U}, _δv::Array{TaylorN{Taylor1{U}},1}, - params, jacobianfunc!; maxsteps::Int=500) where {T<:Real, U<:Number} - - # Allocation - order = get_order(t) - nn = length(trange) - dof = length(q0) - nx0 = length(x) - x0 = getcoeff.(x, 0) - xv = Array{U}(undef, dof, nn) - fill!(xv, U(NaN)) - λ = Array{U}(undef, dof, nn) - fill!(λ, U(NaN)) - λtsum = similar(q0) - q1 = similar(q0) - - # Initial conditions - @inbounds t0, t1, tmax = trange[1], trange[2], trange[end] - @inbounds t[0] = t0 - sign_tstep = copysign(1, tmax-t0) - t00 = t0 - tspan = zero(T) - @inbounds for ind in eachindex(q0) - xv[ind,1] = q0[ind] - λ[ind,1] = zero(U) - λtsum[ind] = zero(U) - end - - #Allocate auxiliary arrays - xaux = Array{Taylor1{U}}(undef, nx0) - δx = Array{TaylorN{Taylor1{U}}}(undef, dof) - dδx = Array{TaylorN{Taylor1{U}}}(undef, dof) - jac = Array{Taylor1{U}}(undef, dof, dof) - varsaux = Array{Taylor1{U}}(undef, dof, dof, dof) - fill!(jac, zero(x[1])) - QH = Array{U}(undef, dof, dof) - RH = Array{U}(undef, dof, dof) - aⱼ = Array{U}(undef, dof ) - qᵢ = similar(aⱼ) - vⱼ = similar(aⱼ) - - # Integration - iter = 2 - nsteps = 1 - while sign_tstep*t0 < sign_tstep*tmax - δt = lyap_taylorstep!(f!, t, x, dx, xaux, δx, dδx, jac, abstol, _δv, - varsaux, params, jacobianfunc!) # δ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) # Update x0 - tnext = t0+δt - # # Evaluate solution at times within convergence radius - while t1 < tnext - evaluate!(x[1:dof], t1-t0, q1) - @inbounds xv[:,iter] .= q1 - for ind in eachindex(jt) - @inbounds jt[ind] = evaluate(x[dof+ind], δt) - end - modifiedGS!( jt, QH, RH, aⱼ, qᵢ, vⱼ ) - tspan = t1-t00 - @inbounds for ind in eachindex(q0) - λ[ind,iter] = (λtsum[ind]+log(RH[ind,ind]))/tspan - end - iter += 1 - @inbounds t1 = trange[iter] - end - if δt == tmax-t0 - @inbounds xv[:,iter] .= x0[1:dof] - for ind in eachindex(jt) - @inbounds jt[ind] = x0[dof+ind] - end - modifiedGS!( jt, QH, RH, aⱼ, qᵢ, vⱼ ) - tspan = tmax-t00 - @inbounds for ind in eachindex(q0) - λ[ind,iter] = (λtsum[ind]+log(RH[ind,ind]))/tspan - end - break - end - - for ind in eachindex(jt) - @inbounds jt[ind] = x0[dof+ind] - end - modifiedGS!( jt, QH, RH, aⱼ, qᵢ, vⱼ ) - - t0 = tnext - @inbounds t[0] = t0 - nsteps += 1 - @inbounds for ind in eachindex(q0) - λtsum[ind] += log(RH[ind,ind]) - end - for ind in eachindex(QH) - @inbounds x0[dof+ind] = QH[ind] - end - x .= Taylor1.( x0, order ) - if nsteps > maxsteps - @warn(""" - Maximum number of integration steps reached; exiting. - """) - break - end - end + # Re-initialize the Taylor1 expansions + t = trange[1] + Taylor1( T, order ) + x .= Taylor1.( x0, order ) + tmpTaylor, arrTaylor = rv.v0, rv.v1 + return _lyap_taylorinteg!(f!, t, x, dx, q0, trange, abstol, jt, _dv, rv, + params, jacobianfunc!; parse_eqs, maxsteps) - return transpose(xv), transpose(λ) end function _lyap_taylorinteg!(f!, t::Taylor1{T}, x::Array{Taylor1{U},1}, dx::Array{Taylor1{U},1}, q0::Array{U,1}, trange::AbstractVector{T}, abstol::T, jt::Matrix{U}, _δv::Array{TaylorN{Taylor1{U}},1}, - rv::RetAlloc{Taylor1{U}}, params, jacobianfunc!; maxsteps::Int=500) where {T<:Real, U<:Number} + rv::RetAlloc{Taylor1{U}}, params, jacobianfunc!; parse_eqs::Bool=true, + maxsteps::Int=500) where {T<:Real, U<:Number} # Allocation order = get_order(t) @@ -605,6 +432,8 @@ function _lyap_taylorinteg!(f!, t::Taylor1{T}, x::Array{Taylor1{U},1}, dx::Array end #Allocate auxiliary arrays + nx0 = length(x0) + xaux = Array{Taylor1{U}}(undef, nx0) δx = Array{TaylorN{Taylor1{U}}}(undef, dof) dδx = Array{TaylorN{Taylor1{U}}}(undef, dof) jac = Array{Taylor1{U}}(undef, dof, dof) @@ -620,7 +449,7 @@ function _lyap_taylorinteg!(f!, t::Taylor1{T}, x::Array{Taylor1{U},1}, dx::Array iter = 2 nsteps = 1 while sign_tstep*t0 < sign_tstep*tmax - δt = lyap_taylorstep!(f!, t, x, dx, δx, dδx, jac, abstol, _δv, + δt = lyap_taylorstep!(Val(parse_eqs), f!, t, x, dx, xaux, δx, dδx, jac, abstol, _δv, varsaux, params, rv, jacobianfunc!) # δ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)) @@ -677,5 +506,5 @@ function _lyap_taylorinteg!(f!, t::Taylor1{T}, x::Array{Taylor1{U},1}, dx::Array end end - return transpose(xv), transpose(λ) + return build_lyap_solution(trange, xv, λ) end diff --git a/src/parse_eqs.jl b/src/parse_eqs.jl index c4f85979a..5a83acaf4 100644 --- a/src/parse_eqs.jl +++ b/src/parse_eqs.jl @@ -46,13 +46,15 @@ end """ - RetAlloc{Taylor1{T}} + RetAlloc{T <: Number} Struct related to the returned variables that are pre-allocated when `@taylorize` is used. - - `v0` : Vector{Taylor1{T}} - - `v1` : Vector{Vector{Taylor1{T}}} - + - `v0` : Array{T,1} + - `v1` : Vector{Array{T,1}} + - `v2` : Vector{Array{T,2}} + - `v3` : Vector{Array{T,3}} + - `v4` : Vector{Array{T,4}} """ struct RetAlloc{T <: Number} v0 :: Array{T,1} @@ -810,7 +812,7 @@ function _replacecalls!(bkkeep::BookKeeping, fnold::Expr, newvar::Symbol) def_fnexpr = Expr(:block, :(TaylorSeries.zero!(_res)), - :(_res.coeffs[1] = $(def_fnexpr.args[2])) ) + :(_res[0] = $(def_fnexpr.args[2])) ) def_fnexpr = subs(def_fnexpr, Dict(:_res => newvar, :_arg1 => :(constant_term($(newarg1))), :_k => :ord)) # def_fnexpr = Expr(:block, diff --git a/src/rootfinding.jl b/src/rootfinding.jl index fbe904fc5..e560d6dfd 100644 --- a/src/rootfinding.jl +++ b/src/rootfinding.jl @@ -1,3 +1,46 @@ +### `build_solution` method for root-finding + +@doc doc""" + build_solution(t, x, p, tevents, xevents, gresids, nsteps, nevents) + build_solution(t, x, tevents, xevents, gresids, nsteps, nevents) + +Helper function to build a [`TaylorSolution`](@ref) from a call to a +root-finding method of [`taylorinteg`](@ref). + +""" +build_solution(t::AbstractVector{T}, + x::Matrix{U}, + p::P, + tevents::AbstractVector{U}, + xevents::Matrix{U}, + gresids::AbstractVector{U}, + nsteps::Int, + nevents::Int) where {T, U, P<:Union{Nothing, Matrix{Taylor1{U}}}} = + TaylorSolution(arraysol(t, nsteps), + arraysol(x, nsteps), + arraysol(p, nsteps-1), + arraysol(tevents, nevents-1), + arraysol(xevents, nevents-1), + arraysol(gresids, nevents-1), + nothing) + +#### `build_solution` method for root-finding with time-ranges + +build_solution(t::AbstractVector{T}, + x::Matrix{U}, + tevents::AbstractVector{U}, + xevents::Matrix{U}, + gresids::AbstractVector{U}, + nevents::Int) where {T, U} = + TaylorSolution(t, + transpose(x), + nothing, + arraysol(tevents, nevents-1), + arraysol(xevents, nevents-1), + arraysol(gresids, nevents-1), + nothing) + + """ surfacecrossing(g_old, g_now, eventorder::Int) @@ -101,8 +144,8 @@ end """ taylorinteg(f, g, x0, t0, tmax, order, abstol, params[=nothing]; kwargs... ) - taylorinteg(f, g, x0, t0, tmax, order, abstol, Val(false), params[=nothing]; kwargs... ) - taylorinteg(f, g, x0, t0, tmax, order, abstol, Val(true), params[=nothing]; kwargs... ) + taylorinteg(f, g, x0, t0, tmax, order, abstol, params[=nothing]; kwargs... ) + taylorinteg(f, g, x0, t0, tmax, order, abstol, params[=nothing]; kwargs... ) taylorinteg(f, g, x0, trange, order, abstol, params[=nothing]; kwargs... ) Root-finding method of `taylorinteg`. @@ -115,8 +158,7 @@ tuple (cond1, cond2). Then, `taylorinteg` attempts to find that root (or event, or crossing) by performing a Newton-Raphson process. When called with the `eventorder=n` keyword argument, `taylorinteg` searches for the roots of the `n`-th derivative of `cond2`, which is computed via automatic -differentiation. When the method used involves `Val(true)`, it also -outputs the Taylor polynomial solutions obtained at each time step. +differentiation. The current keyword argument are: - `maxsteps[=500]`: maximum number of integration steps. @@ -126,6 +168,7 @@ The current keyword argument are: - `newtoniter[=10]`: maximum Newton-Raphson iterations per detected root. - `nrabstol[=eps(T)]`: allowed tolerance for the Newton-Raphson process; T is the common type of `t0`, `tmax` (or `eltype(trange)`) and `abstol`. +- `dense[=true]`: output the Taylor polynomial expansion at each time step. ## Examples: @@ -144,240 +187,33 @@ g(dx, x, params, t) = (true, x[2]) x0 = [1.3, 0.0] # find the roots of `g` along the solution -tv, xv, tvS, xvS, gvS = taylorinteg(pendulum!, g, x0, 0.0, 22.0, 28, 1.0E-20) +sol = taylorinteg(pendulum!, g, x0, 0.0, 22.0, 28, 1.0E-20) # find the roots of the 2nd derivative of `g` along the solution -tv, xv, tvS, xvS, gvS = taylorinteg(pendulum!, g, x0, 0.0, 22.0, 28, 1.0E-20; eventorder=2) +sol = taylorinteg(pendulum!, g, x0, 0.0, 22.0, 28, 1.0E-20; eventorder=2) -# find the roots of `g` along the solution, with dense solution output `psol` -tv, xv, psol, tvS, xvS, gvS = taylorinteg(pendulum!, g, x0, 0.0, 22.0, 28, 1.0E-20, Val(true)) +# find the roots of `g` along the solution, with dense solution output +sol = taylorinteg(pendulum!, g, x0, 0.0, 22.0, 28, 1.0E-20, dense=true) # times at which the solution will be returned tv = 0.0:1.0:22.0 # find the roots of `g` along the solution; return the solution *only* at each value of `tv` -xv, tvS, xvS, gvS = taylorinteg(pendulum!, g, x0, tv, 28, 1.0E-20) +sol = taylorinteg(pendulum!, g, x0, tv, 28, 1.0E-20) # find the roots of the 2nd derivative of `g` along the solution; return the solution *only* at each value of `tv` -xv, tvS, xvS, gvS = taylorinteg(pendulum!, g, x0, tv, 28, 1.0E-20; eventorder=2) +sol = taylorinteg(pendulum!, g, x0, tv, 28, 1.0E-20; eventorder=2) ``` """ -taylorinteg(f!, g, q0::Array{U,1}, t0::T, tmax::T, - order::Int, abstol::T, params = nothing; maxsteps::Int=500, parse_eqs::Bool=true, - eventorder::Int=0, newtoniter::Int=10, nrabstol::T=eps(T)) where {T <: Real,U <: Number} = -taylorinteg(f!, g, q0, t0, tmax, order, abstol, Val(false), params; maxsteps, parse_eqs, - eventorder, newtoniter, nrabstol) - -for V in (:(Val{true}), :(Val{false})) - @eval begin - function taylorinteg(f!, g, q0::Array{U,1}, t0::T, tmax::T, - order::Int, abstol::T, ::$V, params = nothing; maxsteps::Int=500, parse_eqs::Bool=true, - eventorder::Int=0, newtoniter::Int=10, nrabstol::T=eps(T)) where {T <: Real,U <: Number} - - @assert order ≥ eventorder "`eventorder` must be less than or equal to `order`" - - # Initialize the vector of Taylor1 expansions - dof = length(q0) - t = t0 + Taylor1( T, order ) - x = Array{Taylor1{U}}(undef, dof) - dx = Array{Taylor1{U}}(undef, dof) - @inbounds for i in eachindex(q0) - x[i] = Taylor1( q0[i], order ) - dx[i] = Taylor1( zero(q0[i]), order ) - end - - # Determine if specialized jetcoeffs! method exists - parse_eqs, rv = _determine_parsing!(parse_eqs, f!, t, x, dx, params) - - if parse_eqs - # Re-initialize the Taylor1 expansions - t = t0 + Taylor1( T, order ) - x .= Taylor1.( q0, order ) - return _taylorinteg!(f!, g, t, x, dx, q0, t0, tmax, abstol, rv, $V(), params; - maxsteps, eventorder, newtoniter, nrabstol) - else - return _taylorinteg!(f!, g, t, x, dx, q0, t0, tmax, abstol, $V(), params; - maxsteps, eventorder, newtoniter, nrabstol) - end - end - - function _taylorinteg!(f!, g, t::Taylor1{T}, x::Array{Taylor1{U},1}, dx::Array{Taylor1{U},1}, - q0::Array{U,1}, t0::T, tmax::T, abstol::T, ::$V, params; - maxsteps::Int=500, eventorder::Int=0, newtoniter::Int=10, nrabstol::T=eps(T)) where {T <: Real,U <: Number} - - # Allocation - tv = Array{T}(undef, maxsteps+1) - dof = length(q0) - xv = Array{U}(undef, dof, maxsteps+1) - if $V == Val{true} - psol = Array{Taylor1{U}}(undef, dof, maxsteps) - end - xaux = Array{Taylor1{U}}(undef, dof) - - # Initial conditions - order = get_order(t) - @inbounds t[0] = t0 - x0 = deepcopy(q0) - x .= Taylor1.(q0, order) - dx .= zero.(x) - @inbounds tv[1] = t0 - @inbounds xv[:,1] .= q0 - sign_tstep = copysign(1, tmax-t0) - - # Some auxiliary arrays for root-finding/event detection/Poincaré surface of section evaluation - g_tupl = g(dx, x, params, t) - g_tupl_old = g(dx, x, params, t) - δt = zero(x[1]) - δt_old = zero(x[1]) - - x_dx = vcat(x, dx) - g_dg = vcat(g_tupl[2], g_tupl_old[2]) - x_dx_val = Array{U}(undef, length(x_dx) ) - g_dg_val = vcat(evaluate(g_tupl[2]), evaluate(g_tupl_old[2])) - - tvS = Array{U}(undef, maxsteps+1) - xvS = similar(xv) - gvS = similar(tvS) - - - # Integration - nsteps = 1 - nevents = 1 #number of detected events - while sign_tstep*t0 < sign_tstep*tmax - δt_old = δt - δt = taylorstep!(f!, t, x, dx, xaux, abstol, params) # δ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 - if $V == Val{true} - # Store the Taylor polynomial solution - @inbounds psol[:,nsteps] .= deepcopy.(x) - end - g_tupl = g(dx, x, params, t) - nevents = findroot!(t, x, dx, g_tupl_old, g_tupl, eventorder, - tvS, xvS, gvS, t0, δt_old, x_dx, x_dx_val, g_dg, g_dg_val, - nrabstol, newtoniter, nevents) - g_tupl_old = deepcopy(g_tupl) - for i in eachindex(x0) - @inbounds x[i][0] = x0[i] - 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 - - if $V == Val{true} - return view(tv,1:nsteps), view(transpose(view(xv,:,1:nsteps)),1:nsteps,:), - view(transpose(view(psol, :, 1:nsteps-1)), 1:nsteps-1, :), view(tvS,1:nevents-1), view(transpose(view(xvS,:,1:nevents-1)),1:nevents-1,:), view(gvS,1:nevents-1) - elseif $V == Val{false} - return return view(tv,1:nsteps), view(transpose(view(xv,:,1:nsteps)),1:nsteps,:), view(tvS,1:nevents-1), view(transpose(view(xvS,:,1:nevents-1)),1:nevents-1,:), view(gvS,1:nevents-1) - end - end - function _taylorinteg!(f!, g, t::Taylor1{T}, x::Array{Taylor1{U},1}, dx::Array{Taylor1{U},1}, - q0::Array{U,1}, t0::T, tmax::T, abstol::T, rv::RetAlloc{Taylor1{U}}, ::$V, params; - maxsteps::Int=500, eventorder::Int=0, newtoniter::Int=10, nrabstol::T=eps(T)) where {T <: Real,U <: Number} - - # Allocation - tv = Array{T}(undef, maxsteps+1) - dof = length(q0) - xv = Array{U}(undef, dof, maxsteps+1) - if $V == Val{true} - psol = Array{Taylor1{U}}(undef, dof, maxsteps) - end - - # Initial conditions - order = get_order(t) - @inbounds t[0] = t0 - x0 = deepcopy(q0) - x .= Taylor1.(q0, order) - dx .= zero.(x) - @inbounds tv[1] = t0 - @inbounds xv[:,1] .= q0 - sign_tstep = copysign(1, tmax-t0) - - # Some auxiliary arrays for root-finding/event detection/Poincaré surface of section evaluation - g_tupl = g(dx, x, params, t) - g_tupl_old = g(dx, x, params, t) - δt = zero(x[1]) - δt_old = zero(x[1]) - - x_dx = vcat(x, dx) - g_dg = vcat(g_tupl[2], g_tupl_old[2]) - x_dx_val = Array{U}(undef, length(x_dx) ) - g_dg_val = vcat(evaluate(g_tupl[2]), evaluate(g_tupl_old[2])) - - tvS = Array{U}(undef, maxsteps+1) - xvS = similar(xv) - gvS = similar(tvS) - - - # Integration - nsteps = 1 - nevents = 1 #number of detected events - while sign_tstep*t0 < sign_tstep*tmax - δt_old = δt - # δt = taylorstep!(f!, t, x, dx, xaux, abstol, params) # δt is positive! - δt = taylorstep!(f!, t, x, dx, abstol, params, rv) # δ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 - if $V == Val{true} - # Store the Taylor polynomial solution - @inbounds psol[:,nsteps] .= deepcopy.(x) - end - g_tupl = g(dx, x, params, t) - nevents = findroot!(t, x, dx, g_tupl_old, g_tupl, eventorder, - tvS, xvS, gvS, t0, δt_old, x_dx, x_dx_val, g_dg, g_dg_val, - nrabstol, newtoniter, nevents) - g_tupl_old = deepcopy(g_tupl) - for i in eachindex(x0) - @inbounds x[i][0] = x0[i] - 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 - - if $V == Val{true} - return view(tv,1:nsteps), view(transpose(view(xv,:,1:nsteps)),1:nsteps,:), - view(transpose(view(psol, :, 1:nsteps-1)), 1:nsteps-1, :), view(tvS,1:nevents-1), view(transpose(view(xvS,:,1:nevents-1)),1:nevents-1,:), view(gvS,1:nevents-1) - elseif $V == Val{false} - return view(tv,1:nsteps), view(transpose(view(xv,:,1:nsteps)),1:nsteps,:), view(tvS,1:nevents-1), view(transpose(view(xvS,:,1:nevents-1)),1:nevents-1,:), view(gvS,1:nevents-1) - end - end - end -end - -function taylorinteg(f!, g, q0::Array{U,1}, trange::AbstractVector{T}, +function taylorinteg(f!, g, q0::Array{U,1}, t0::T, tmax::T, order::Int, abstol::T, params = nothing; maxsteps::Int=500, parse_eqs::Bool=true, - eventorder::Int=0, newtoniter::Int=10, nrabstol::T=eps(T)) where {T <: Real,U <: Number} + dense::Bool=true, eventorder::Int=0, newtoniter::Int=10, nrabstol::T=eps(T)) where {T <: Real,U <: Number} @assert order ≥ eventorder "`eventorder` must be less than or equal to `order`" - # Check if trange is increasingly or decreasingly sorted - @assert (issorted(trange) || - issorted(reverse(trange))) "`trange` or `reverse(trange)` must be sorted" - # Initialize the vector of Taylor1 expansions dof = length(q0) - @inbounds t0 = trange[1] t = t0 + Taylor1( T, order ) x = Array{Taylor1{U}}(undef, dof) dx = Array{Taylor1{U}}(undef, dof) @@ -389,45 +225,40 @@ function taylorinteg(f!, g, q0::Array{U,1}, trange::AbstractVector{T}, # Determine if specialized jetcoeffs! method exists parse_eqs, rv = _determine_parsing!(parse_eqs, f!, t, x, dx, params) - if parse_eqs - # Re-initialize the Taylor1 expansions - t = t0 + Taylor1( T, order ) - x .= Taylor1.(q0, order) - return _taylorinteg!(f!, g, t, x, dx, q0, trange, abstol, rv, params; - maxsteps, eventorder, newtoniter, nrabstol) - else - return _taylorinteg!(f!, g, t, x, dx, q0, trange, abstol,params; - maxsteps, eventorder, newtoniter, nrabstol) - end + # Re-initialize the Taylor1 expansions + t = t0 + Taylor1( T, order ) + x .= Taylor1.( q0, order ) + return _taylorinteg!(Val(dense), f!, g, t, x, dx, q0, t0, tmax, abstol, rv, params; + parse_eqs, maxsteps, eventorder, newtoniter, nrabstol) end -function _taylorinteg!(f!, g, t::Taylor1{T}, x::Array{Taylor1{U},1}, dx::Array{Taylor1{U},1}, - q0::Array{U,1}, trange::AbstractVector{T}, abstol::T, params; - maxsteps::Int=500, eventorder::Int=0, newtoniter::Int=10, nrabstol::T=eps(T)) where {T <: Real,U <: Number} +function _taylorinteg!(dense::Val{D}, f!, g, t::Taylor1{T}, x::Array{Taylor1{U},1}, dx::Array{Taylor1{U},1}, + q0::Array{U,1}, t0::T, tmax::T, abstol::T, rv::RetAlloc{Taylor1{U}}, params; + parse_eqs::Bool=true, maxsteps::Int=500, eventorder::Int=0, + newtoniter::Int=10, nrabstol::T=eps(T)) where {T <: Real,U <: Number, D} # Allocation - nn = length(trange) + tv = Array{T}(undef, maxsteps+1) dof = length(q0) - x0 = similar(q0, eltype(q0), dof) - fill!(x0, T(NaN)) - xv = Array{eltype(q0)}(undef, dof, nn) - for ind in 1:nn - @inbounds xv[:,ind] .= x0 - end + xv = Array{U}(undef, dof, maxsteps+1) + psol = init_psol(dense, xv, x) xaux = Array{Taylor1{U}}(undef, dof) # Initial conditions - @inbounds t0, t1, tmax = trange[1], trange[2], trange[end] - sign_tstep = copysign(1, tmax-t0) + order = get_order(t) + @inbounds t[0] = t0 x0 = deepcopy(q0) - x1 = similar(x0) - @inbounds xv[:,1] .= q0 + x .= Taylor1.(q0, order) + dx .= zero.(x) + @inbounds tv[1] = t0 + @inbounds xv[:,1] .= deepcopy(q0) + sign_tstep = copysign(1, tmax-t0) # Some auxiliary arrays for root-finding/event detection/Poincaré surface of section evaluation g_tupl = g(dx, x, params, t) g_tupl_old = g(dx, x, params, t) - δt = zero(U) - δt_old = zero(U) + δt = zero(x[1]) + δt_old = zero(x[1]) x_dx = vcat(x, dx) g_dg = vcat(g_tupl[2], g_tupl_old[2]) @@ -439,39 +270,29 @@ function _taylorinteg!(f!, g, t::Taylor1{T}, x::Array{Taylor1{U},1}, dx::Array{T gvS = similar(tvS) # Integration - iter = 2 nsteps = 1 nevents = 1 #number of detected events while sign_tstep*t0 < sign_tstep*tmax δt_old = δt - # δt = taylorstep!(f!, t, x, dx, xaux, abstol, params, tmpTaylor, arrTaylor, parse_eqs) # δt is positive! - δt = taylorstep!(f!, t, x, dx, xaux, abstol, params) # δt is positive! + δt = taylorstep!(Val(parse_eqs), f!, t, x, dx, xaux, abstol, params, rv) # δ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 - tnext = t0+δt - # Evaluate solution at times within convergence radius - while sign_tstep*t1 < sign_tstep*tnext - evaluate!(x, t1-t0, x1) - @inbounds xv[:,iter] .= x1 - iter += 1 - @inbounds t1 = trange[iter] - end - if δt == tmax-t0 - @inbounds xv[:,iter] .= x0 - break - end + set_psol!(dense, psol, nsteps, x) # Store the Taylor polynomial solution g_tupl = g(dx, x, params, t) nevents = findroot!(t, x, dx, g_tupl_old, g_tupl, eventorder, tvS, xvS, gvS, t0, δt_old, x_dx, x_dx_val, g_dg, g_dg_val, nrabstol, newtoniter, nevents) g_tupl_old = deepcopy(g_tupl) - for i in eachindex(x0) - @inbounds x[i][0] = x0[i] + @inbounds for i in eachindex(x0) + x[i][0] = x0[i] + TaylorSeries.zero!(dx[i], 0) end - t0 = tnext + t0 += δt @inbounds t[0] = t0 nsteps += 1 + @inbounds tv[nsteps] = t0 + @inbounds xv[:,nsteps] .= deepcopy(x0) if nsteps > maxsteps @warn(""" Maximum number of integration steps reached; exiting. @@ -480,12 +301,43 @@ function _taylorinteg!(f!, g, t::Taylor1{T}, x::Array{Taylor1{U},1}, dx::Array{T end end - return transpose(xv), view(tvS,1:nevents-1), view(transpose(view(xvS,:,1:nevents-1)),1:nevents-1,:), view(gvS,1:nevents-1) + return build_solution(tv, xv, psol, tvS, xvS, gvS, nsteps, nevents) +end + +function taylorinteg(f!, g, q0::Array{U,1}, trange::AbstractVector{T}, + order::Int, abstol::T, params = nothing; maxsteps::Int=500, parse_eqs::Bool=true, + eventorder::Int=0, newtoniter::Int=10, nrabstol::T=eps(T)) where {T <: Real,U <: Number} + + @assert order ≥ eventorder "`eventorder` must be less than or equal to `order`" + + # Check if trange is increasingly or decreasingly sorted + @assert (issorted(trange) || + issorted(reverse(trange))) "`trange` or `reverse(trange)` must be sorted" + + # Initialize the vector of Taylor1 expansions + dof = length(q0) + @inbounds t0 = trange[1] + t = t0 + Taylor1( T, order ) + x = Array{Taylor1{U}}(undef, dof) + dx = Array{Taylor1{U}}(undef, dof) + @inbounds for i in eachindex(q0) + x[i] = Taylor1( q0[i], order ) + dx[i] = Taylor1( zero(q0[i]), order ) + end + + # Determine if specialized jetcoeffs! method exists + parse_eqs, rv = _determine_parsing!(parse_eqs, f!, t, x, dx, params) + + # Re-initialize the Taylor1 expansions + t = t0 + Taylor1( T, order ) + x .= Taylor1.(q0, order) + return _taylorinteg!(f!, g, t, x, dx, q0, trange, abstol, rv, params; + parse_eqs, maxsteps, eventorder, newtoniter, nrabstol) end function _taylorinteg!(f!, g, t::Taylor1{T}, x::Array{Taylor1{U},1}, dx::Array{Taylor1{U},1}, q0::Array{U,1}, trange::AbstractVector{T}, abstol::T, rv::RetAlloc{Taylor1{U}}, params; - maxsteps::Int=500, eventorder::Int=0, newtoniter::Int=10, nrabstol::T=eps(T)) where {T <: Real,U <: Number} + parse_eqs::Bool=true, maxsteps::Int=500, eventorder::Int=0, newtoniter::Int=10, nrabstol::T=eps(T)) where {T <: Real,U <: Number} # Allocation nn = length(trange) @@ -494,15 +346,16 @@ function _taylorinteg!(f!, g, t::Taylor1{T}, x::Array{Taylor1{U},1}, dx::Array{T fill!(x0, T(NaN)) xv = Array{eltype(q0)}(undef, dof, nn) for ind in 1:nn - @inbounds xv[:,ind] .= x0 + @inbounds xv[:,ind] .= deepcopy(x0) end + xaux = Array{Taylor1{U}}(undef, dof) # Initial conditions @inbounds t0, t1, tmax = trange[1], trange[2], trange[end] sign_tstep = copysign(1, tmax-t0) x0 = deepcopy(q0) x1 = similar(x0) - @inbounds xv[:,1] .= q0 + @inbounds xv[:,1] .= deepcopy(q0) # Some auxiliary arrays for root-finding/event detection/Poincaré surface of section evaluation g_tupl = g(dx, x, params, t) @@ -525,7 +378,7 @@ function _taylorinteg!(f!, g, t::Taylor1{T}, x::Array{Taylor1{U},1}, dx::Array{T nevents = 1 #number of detected events while sign_tstep*t0 < sign_tstep*tmax δt_old = δt - δt = taylorstep!(f!, t, x, dx, abstol, params, rv) # δt is positive! + δt = taylorstep!(Val(parse_eqs), f!, t, x, dx, xaux, abstol, params, rv) # δ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 @@ -533,12 +386,12 @@ function _taylorinteg!(f!, g, t::Taylor1{T}, x::Array{Taylor1{U},1}, dx::Array{T # Evaluate solution at times within convergence radius while sign_tstep*t1 < sign_tstep*tnext evaluate!(x, t1-t0, x1) - @inbounds xv[:,iter] .= x1 + @inbounds xv[:,iter] .= deepcopy(x1) iter += 1 @inbounds t1 = trange[iter] end if δt == tmax-t0 - @inbounds xv[:,iter] .= x0 + @inbounds xv[:,iter] .= deepcopy(x0) break end g_tupl = g(dx, x, params, t) @@ -546,8 +399,8 @@ function _taylorinteg!(f!, g, t::Taylor1{T}, x::Array{Taylor1{U},1}, dx::Array{T tvS, xvS, gvS, t0, δt_old, x_dx, x_dx_val, g_dg, g_dg_val, nrabstol, newtoniter, nevents) g_tupl_old = deepcopy(g_tupl) - for i in eachindex(x0) - @inbounds x[i][0] = x0[i] + @inbounds for i in eachindex(x0) + x[i][0] = x0[i] end t0 = tnext @inbounds t[0] = t0 @@ -560,5 +413,5 @@ function _taylorinteg!(f!, g, t::Taylor1{T}, x::Array{Taylor1{U},1}, dx::Array{T end end - return transpose(xv), view(tvS,1:nevents-1), view(transpose(view(xvS,:,1:nevents-1)),1:nevents-1,:), view(gvS,1:nevents-1) + return build_solution(trange, xv, tvS, xvS, gvS, nevents) end diff --git a/test/aqua.jl b/test/aqua.jl index cde247e64..dde52925a 100644 --- a/test/aqua.jl +++ b/test/aqua.jl @@ -28,6 +28,6 @@ end Aqua.test_all( TaylorIntegration; ambiguities=false, - stale_deps=(ignore=[:DiffEqBase, :RecursiveArrayTools, :Requires, :StaticArrays],), + stale_deps=(ignore=[:DiffEqBase, :RecursiveArrayTools, :Requires, :StaticArrays],) ) end diff --git a/test/bigfloats.jl b/test/bigfloats.jl index 15bd30a53..9ede11a9b 100644 --- a/test/bigfloats.jl +++ b/test/bigfloats.jl @@ -22,8 +22,10 @@ import Logging: Warn # T is the pendulum's librational period == 4Elliptic.K(sin(q0[1]/2)^2) # we will evaluate the elliptic integral K using TaylorIntegration.jl: G(x, p, t) = 1/sqrt(1-((sin(big"1.3"/2))^2)*(sin(t)^2)) # K elliptic integral kernel - tvk, xvk = (@test_logs min_level=Logging.Warn taylorinteg( + solk = (@test_logs min_level=Logging.Warn taylorinteg( G, 0.0, 0.0, BigFloat(π)/2, _order, _abstol)) + tvk = solk.t + xvk = solk.x @test eltype(tvk) == BigFloat @test eltype(xvk) == BigFloat T = 4xvk[end] # T = 4Elliptic.K(sin(q0[1]/2)^2) @@ -32,8 +34,10 @@ import Logging: Warn t0 = 0.0 #the initial time - tv, xv = (@test_logs (Warn, max_iters_reached()) taylorinteg( + sol = (@test_logs (Warn, max_iters_reached()) taylorinteg( pendulum!, q0, t0, T, _order, _abstol; maxsteps=1)) + tv = sol.t + xv = sol.x @test eltype(tv) == BigFloat @test eltype(xv) == BigFloat @test length(tv) == 2 @@ -41,8 +45,10 @@ import Logging: Warn @test length(xv[:,2]) == 2 #note that T is a BigFloat - tv, xv = (@test_logs min_level=Logging.Warn taylorinteg( + sol = (@test_logs min_level=Logging.Warn taylorinteg( pendulum!, q0, t0, T, _order, _abstol)) + tv = sol.t + xv = sol.x @test length(tv) < 501 @test length(xv[:,1]) < 501 @test length(xv[:,2]) < 501 diff --git a/test/common.jl b/test/common.jl index 8d3ddcc94..b38d1541e 100644 --- a/test/common.jl +++ b/test/common.jl @@ -9,6 +9,11 @@ import Logging: Warn local TI = TaylorIntegration max_iters_reached() = "Maximum number of integration steps reached; exiting.\n" + larger_maxiters_needed() = "Interrupted. Larger maxiters is needed. If you "* + "are using an integrator for non-stiff ODEs or an automatic switching "* + "algorithm (the default), you may want to consider using a method for "* + "stiff equations. See the solver pages for more details (e.g. "* + "https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Stiff-Problems)." f(u,p,t) = u g(u,p,t) = cos(t) @@ -95,10 +100,96 @@ import Logging: Warn u0 = [1.0; 0.0] prob = ODEProblem(harmosc!, u0, tspan) sol = (@test_logs min_level=Logging.Warn solve(prob, TaylorMethod(order), abstol=abstol)) - tv1, xv1 = (@test_logs min_level=Logging.Warn taylorinteg( + solti = (@test_logs min_level=Logging.Warn taylorinteg( harmosc!, u0, tspan[1], tspan[2], order, abstol)) - @test sol.t == tv1 - @test xv1[end,:] == sol.u[end] + @test solti.t == sol.t + @test solti.x[end,:] == sol.u[end] + @test sol.t == solti.t + @test solti.x[end,:] == sol.u[end] + @test DiffEqBase.interp_summary(sol.interp) == "Taylor series polynomial evaluation" + ### backwards integration + tspanb = (0.0, -10pi) + probb = ODEProblem(harmosc!, u0, tspanb) + solb = (@test_logs min_level=Logging.Warn solve(probb, TaylorMethod(order), abstol=abstol)) + soltib = (@test_logs min_level=Logging.Warn taylorinteg( + harmosc!, u0, tspanb[1], tspanb[2], order, abstol)) + @test soltib.t == solb.t + @test soltib.x[end,:] == solb.u[end] + @test solb.t == soltib.t + @test soltib.x[end,:] == solb.u[end] + @test DiffEqBase.interp_summary(solb.interp) == "Taylor series polynomial evaluation" + ### with dense output + sol2 = (@test_logs min_level=Logging.Warn taylorinteg( + harmosc!, u0, tspan[1], tspan[2], order, abstol, dense=true)) + tv2, xv2, psol2 = sol2.t, sol2.x, sol2.p + Θ = rand() + dt_test = (tv2[end]-tv2[end-1]) + t_test = tv2[end-1] + Θ*dt_test + @test norm(sol(t_test) .- psol2[end,:]( Θ*dt_test ), Inf) < 1e-14 + Θ = rand() + Θm1 = Θ - 1 + dt_test = (tv2[end-1]-tv2[end-2]) + t_test = tv2[end-2] + Θ*dt_test + @test norm( sol(t_test, idxs=1:1:2) .- psol2[end,:]( Θm1*dt_test ) ) < 1e-14 + # Kepler problem + @taylorize function kepler1!(dq, q, p, t) + local μ = -1.0 + r_p2 = ((q[1]^2)+(q[2]^2)) + r_p3d2 = r_p2^1.5 + dq[1] = q[3] + dq[2] = q[4] + newtonianCoeff = μ / r_p3d2 + dq[3] = q[1] * newtonianCoeff + dq[4] = q[2] * newtonianCoeff + return nothing + end + p = set_variables("ξ", numvars=4, order=2) + q0 = [0.2, 0.0, 0.0, 3.0] + q0TN = q0 + p + tspan = (0.0, 10pi) + (@test_logs (Warn, max_iters_reached()) taylorinteg(kepler1!, q0, tspan[1], tspan[2], order, abstol, + dense=true, maxsteps=1)) + solp = taylorinteg(kepler1!, q0, tspan[1], tspan[2], order, abstol, + dense=true, maxsteps=2000) + tvp, xvp, psolp = solp.t, solp.x, solp.p + prob = ODEProblem(kepler1!, q0, tspan) + @test isinplace(prob) == true + sol1 = solve(prob, TaylorMethod(order), abstol=abstol, maxiters=2000, parse_eqs=false) + (@test_logs (Warn, larger_maxiters_needed()) solve(prob, TaylorMethod(order), abstol=abstol, maxiters=1)) + sol2 = solve(prob, TaylorMethod(order), abstol=abstol, maxiters=2000) + @test sol1.alg.parse_eqs == false + @test sol2.alg.parse_eqs == true + @test tvp == sol2.t + @test sol1.u == sol2.u + @test transpose(Array(sol2)) == xvp + @test norm( solp(tspan[1]) - sol1(tspan[1]), Inf ) < 1e-16 + @test norm( solp(0.1) - sol1(0.1), Inf ) < 1e-13 + @test norm( solp(tspan[2]/2) - sol1(tspan[2]/2), Inf ) < 1e-13 + @test norm( solp(solp.t[end-2] + 0.01) - sol1(sol1.t[end-2] + 0.01), Inf ) < 1e-13 + @test solp(tspan[2]) == sol1(tspan[2]) + + (@test_logs (Warn, max_iters_reached()) taylorinteg(kepler1!, q0TN, tspan[1], tspan[2], order, abstol, + dense = true, maxsteps=1)) + solTN = taylorinteg(kepler1!, q0TN, tspan[1], tspan[2], order, abstol, + dense = true, maxsteps=2000) + tTN, xTN, psolTN = solTN.t, solTN.x, solTN.p + + probTN = ODEProblem(kepler1!, q0TN, tspan) + (@test_logs (Warn,larger_maxiters_needed()) solve(probTN, TaylorMethod(order), abstol=abstol, parse_eqs=true, dense=true, maxiters=1)) + sol2TN = solve(probTN, TaylorMethod(order), abstol=abstol, parse_eqs=true, dense=true, maxiters=2000) + @test sol2TN.alg.parse_eqs == true + @test DiffEqBase.interp_summary(sol2TN.interp) == "Taylor series polynomial evaluation" + @test size(xTN) == size(Array(sol2TN)') + @test xTN == Array(sol2TN)' + @test tTN == sol2TN.t + @test norm( solTN(tspan[1]) - sol2TN(tspan[1]), Inf ) < 1e-14 + @test psolTN[end,:]( tTN[end] - tTN[end-1] ) == sol2TN(sol2TN.t[end]) + @test solTN(solTN.t[end]) == sol2TN(sol2TN.t[end]) + @test psolTN[2,:]( 0.005 - tTN[2] ) == sol2TN(0.005) + @test norm( solTN(0.1) - sol2TN(0.1), Inf ) < 1e-13 + @test norm( solTN(tspan[2]/2)() - sol2TN(tspan[2]/2)(), Inf ) < 1e-15 + @test norm( solTN(tspan[2]/2) - sol2TN(tspan[2]/2), Inf ) < 1e-10 + @test solTN(tspan[2]) == sol2TN(tspan[2]) end @testset "Test discrete callback in common interface" begin @@ -238,11 +329,11 @@ import Logging: Warn @test length(sol1.t) == length(sol2.t) @test sol1.t == sol2.t @test sol1.u == sol2.u - tv, xv = (@test_logs min_level=Logging.Warn taylorinteg( + sol1ti = (@test_logs min_level=Logging.Warn taylorinteg( integ_vec, x0, tspan[1], tspan[2], order, abstol, [1.0])) - @test sol1.t == tv - @test sol1[1,:] == xv[:,1] - @test sol1[2,:] == xv[:,2] + @test sol1.t == sol1ti.t + @test sol1[1,:] == sol1ti.x[:,1] + @test sol1[2,:] == sol1ti.x[:,2] @test norm(sol1.u[end][1] - sin(sol1.t[end]), Inf) < 1.0e-14 @test norm(sol1.u[end][2] - cos(sol1.t[end]), Inf) < 1.0e-14 diff --git a/test/complex.jl b/test/complex.jl index 60bef45a3..453f6f19e 100644 --- a/test/complex.jl +++ b/test/complex.jl @@ -26,33 +26,37 @@ import Logging: Warn ts = 0.0:pi:2pi zsol1 = (@test_logs (Warn, max_iters_reached()) taylorinteg( eqs1, z0, ts, _order, _abstol, maxsteps=1)) - @test size(zsol1) == (length(ts),) + @test size(zsol1.x) == (length(ts),) zsol1 = (@test_logs (Warn, max_iters_reached()) taylorinteg( eqs1, z0, tr, _order, _abstol, maxsteps=1, nothing)) - @test length(zsol1) == length(tr) + @test length(zsol1.x) == length(tr) ta = vec(tr) zsol1 = (@test_logs (Warn, max_iters_reached()) taylorinteg( eqs1, z0, ta, _order, _abstol, maxsteps=1)) - @test length(zsol1) == length(ta) - zsol1 = (@test_logs min_level=Logging.Warn taylorinteg( + @test length(zsol1.x) == length(ta) + sol1 = (@test_logs min_level=Logging.Warn taylorinteg( eqs1, z0, tr, _order, _abstol)) + zsol1 = sol1.x @test zsol1[1] == z0 @test isapprox( zsol1[2] , z0*exp(-tr[2]) ) @test isapprox( zsol1[6], z0*exp(-tr[6]) ) @test isapprox( zsol1[end], z0*exp(-tr[end]) ) - zsol11 = (@test_logs min_level=Logging.Warn taylorinteg( + sol11 = (@test_logs min_level=Logging.Warn taylorinteg( eqs1, z0, ta, _order, _abstol, nothing)) + zsol11 = sol11.x @test zsol11 == zsol1 @test zsol11[1] == z0 @test isapprox( zsol11[2] , z0*exp(-tr[2]) ) @test isapprox( zsol11[6], z0*exp(-tr[6]) ) @test isapprox( zsol11[end], z0*exp(-tr[end]) ) - tt, zsol1t = (@test_logs (Warn, max_iters_reached()) taylorinteg( + sol1t = (@test_logs (Warn, max_iters_reached()) taylorinteg( eqs1, z0, 0.0, 2pi, _order, _abstol, maxsteps=1)) + tt = sol1t.t; zsol1t = sol1t.x @test length(tt) == 2 @test length(zsol1t) == 2 - tt, zsol1t = (@test_logs min_level=Logging.Warn taylorinteg( + sol1t = (@test_logs min_level=Logging.Warn taylorinteg( eqs1, z0, 0.0, 2pi, _order, _abstol)) + tt = sol1t.t; zsol1t = sol1t.x @test zsol1t[1] == z0 @test isapprox( zsol1t[2] , z0*exp(-tt[2]) ) @test isapprox( zsol1t[end], z0*exp(-tt[end]) ) @@ -65,31 +69,35 @@ import Logging: Warn ts = 0.0:pi:2pi zsol2 = (@test_logs (Warn, max_iters_reached()) taylorinteg( eqs2, z0, ts, _order, _abstol, maxsteps=1, nothing)) - @test size(zsol2) == (length(ts),) + @test size(zsol2.x) == (length(ts),) zsol2 = (@test_logs (Warn, max_iters_reached()) taylorinteg( eqs2, z0, tr, _order, _abstol, maxsteps=1)) - @test length(zsol2) == length(tr) + @test length(zsol2.x) == length(tr) ta = vec(tr) zsol2 = (@test_logs (Warn, max_iters_reached()) taylorinteg( eqs2, z0, ta, _order, _abstol, maxsteps=1)) - @test length(zsol2) == length(ta) - zsol2 = (@test_logs min_level=Logging.Warn taylorinteg( + @test length(zsol2.x) == length(ta) + sol2 = (@test_logs min_level=Logging.Warn taylorinteg( eqs2, z0, tr, _order, _abstol)) + zsol2 = sol2.x @test zsol2[1] == z0 @test isapprox( zsol2[3], z0*exp( complex(0.0, tr[3])) ) @test isapprox( zsol2[5], z0*exp( complex(0.0, tr[5])) ) - zsol22 = (@test_logs min_level=Logging.Warn taylorinteg( + sol22 = (@test_logs min_level=Logging.Warn taylorinteg( eqs2, z0, ta, _order, _abstol, nothing)) + zsol22 = sol22.x @test zsol22 == zsol2 @test zsol22[1] == z0 @test isapprox( zsol22[3], z0*exp( complex(0.0, tr[3])) ) @test isapprox( zsol22[5], z0*exp( complex(0.0, tr[5])) ) - tt, zsol2 = (@test_logs (Warn, max_iters_reached()) taylorinteg( + sol2 = (@test_logs (Warn, max_iters_reached()) taylorinteg( eqs2, z0, 0.0, 2pi, _order, _abstol, maxsteps=1)) + tt = sol2.t; zsol2 = sol2.x @test length(tt) == 2 @test length(zsol2) == 2 - tt, zsol2 = (@test_logs min_level=Logging.Warn taylorinteg( + sol2 = (@test_logs min_level=Logging.Warn taylorinteg( eqs2, z0, 0.0, 2pi, _order, _abstol)) + tt = sol2.t; zsol2 = sol2.x @test zsol2[1] == z0 @test isapprox( zsol2[2], z0*exp( complex(0.0, tt[2])) ) @test isapprox( zsol2[end], z0*exp( complex(0.0, tt[end])) ) @@ -102,33 +110,35 @@ import Logging: Warn ts = 0.0:pi:2pi zsol3 = (@test_logs (Warn, max_iters_reached()) taylorinteg( eqs3!, zz0, ts, _order, _abstol, maxsteps=1)) - @test size(zsol3) == (length(ts), length(zz0)) + @test size(zsol3.x) == (length(ts), length(zz0)) zsol3 = (@test_logs (Warn, max_iters_reached()) taylorinteg( eqs3!, zz0, tr, _order, _abstol, nothing, maxsteps=1)) - @test size(zsol3) == (length(tr), length(zz0)) + @test size(zsol3.x) == (length(tr), length(zz0)) ta = vec(tr) zsol3 = (@test_logs (Warn, max_iters_reached()) taylorinteg( eqs3!, zz0, ta, _order, _abstol, maxsteps=1)) - @test size(zsol3) == (length(ta), length(zz0)) + @test size(zsol3.x) == (length(ta), length(zz0)) zsol3 = (@test_logs min_level=Logging.Warn taylorinteg( eqs3!, [z0, z0], tr, _order, _abstol)) - @test isapprox( zsol3[4,1], z0*exp(-tr[4]) ) - @test isapprox( zsol3[7,1], z0*exp(-tr[7]) ) - @test isapprox( zsol3[4,2], z0*exp( complex(0.0, tr[4])) ) - @test isapprox( zsol3[7,2], z0*exp( complex(0.0, tr[7])) ) + @test isapprox( zsol3.x[4,1], z0*exp(-tr[4]) ) + @test isapprox( zsol3.x[7,1], z0*exp(-tr[7]) ) + @test isapprox( zsol3.x[4,2], z0*exp( complex(0.0, tr[4])) ) + @test isapprox( zsol3.x[7,2], z0*exp( complex(0.0, tr[7])) ) zsol33 = (@test_logs min_level=Logging.Warn taylorinteg( eqs3!, [z0, z0], ta, _order, _abstol, nothing)) - @test zsol33 == zsol3 - @test isapprox( zsol33[4,1], z0*exp(-tr[4]) ) - @test isapprox( zsol33[7,1], z0*exp(-tr[7]) ) - @test isapprox( zsol33[4,2], z0*exp( complex(0.0, tr[4])) ) - @test isapprox( zsol33[7,2], z0*exp( complex(0.0, tr[7])) ) - tt, zsol3 = (@test_logs (Warn, max_iters_reached()) taylorinteg( + @test zsol33.x == zsol3.x + @test isapprox( zsol33.x[4,1], z0*exp(-tr[4]) ) + @test isapprox( zsol33.x[7,1], z0*exp(-tr[7]) ) + @test isapprox( zsol33.x[4,2], z0*exp( complex(0.0, tr[4])) ) + @test isapprox( zsol33.x[7,2], z0*exp( complex(0.0, tr[7])) ) + sol3 = (@test_logs (Warn, max_iters_reached()) taylorinteg( eqs3!, zz0, 0.0, 2pi, _order, _abstol, maxsteps=1)) + tt = sol3.t; zsol3 = sol3.x @test length(tt) == 2 @test size(zsol3) == (length(tt), length(zz0)) - tt, zsol3 = (@test_logs min_level=Logging.Warn taylorinteg( + sol3 = (@test_logs min_level=Logging.Warn taylorinteg( eqs3!, zz0, 0.0, 2pi, _order, _abstol)) + tt = sol3.t; zsol3 = sol3.x @test zsol3[1,1] == z0 @test zsol3[1,2] == z0 @test isapprox( zsol3[2,1], z0*exp( -tt[2]) ) diff --git a/test/jettransport.jl b/test/jettransport.jl index 670562f21..10d55664b 100644 --- a/test/jettransport.jl +++ b/test/jettransport.jl @@ -22,14 +22,17 @@ import Logging: Warn x0TN = x0 + p[1] #jet transport initial condition t0=0.0 tmax=0.3 - tvTN, xvTN = (@test_logs (Warn, max_iters_reached()) taylorinteg( + solTN = (@test_logs (Warn, max_iters_reached()) taylorinteg( f, x0TN, t0, tmax, _order, _abstol, maxsteps=1)) + tvTN = solTN.t; xvTN = solTN.x @test size(tvTN) == (2,) @test size(xvTN) == (2,) - tvTN, xvTN = (@test_logs min_level=Logging.Warn taylorinteg( + solTN = (@test_logs min_level=Logging.Warn taylorinteg( f, x0TN, t0, tmax, _order, _abstol)) - tv, xv = (@test_logs min_level=Logging.Warn taylorinteg( + tvTN = solTN.t; xvTN = solTN.x + sol = (@test_logs min_level=Logging.Warn taylorinteg( f, x0, t0, tmax, _order, _abstol)) + tv = sol.t; xv = sol.x exactsol(t, x0, t0) = x0/(1.0-x0*(t-t0)) #the analytical solution δsol = exactsol(tvTN[end], x0TN, t0)-xvTN[end] δcoeffs = map(y->y[1], map(x->x.coeffs, δsol.coeffs)) @@ -43,9 +46,10 @@ import Logging: Warn disp = 0.001*rand() #a small, random displacement x0_disp = x0+disp dv = map(x->[disp], tvTN) #a vector of identical displacements - tv_disp, xv_disp = (@test_logs min_level=Logging.Warn taylorinteg( + sol_disp = (@test_logs min_level=Logging.Warn taylorinteg( f, x0_disp, t0, tmax, _order, _abstol)) xvTN_disp = evaluate.(xvTN, dv) + tv_disp = sol_disp.t; xv_disp = sol_disp.x @test norm(exactsol.(tvTN, x0_disp, t0)-xvTN_disp, Inf) < 1E-12 #analytical vs jet transport @test norm(x0_disp-evaluate(xvTN[1], [disp]), Inf) < 1E-12 @test norm(xv_disp[1]-evaluate(xvTN[1], [disp]), Inf) < 1E-12 @@ -55,14 +59,17 @@ import Logging: Warn y0 = 1.0 #"nominal" initial condition u0 = 0.0 #initial time y0TN = y0 + p[1] #jet transport initial condition - uvTN, yvTN = (@test_logs (Warn, max_iters_reached()) taylorinteg( + solTN = (@test_logs (Warn, max_iters_reached()) taylorinteg( g, y0TN, u0, 10/0.3, _order, _abstol, maxsteps=1)) + uvTN = solTN.t; yvTN = solTN.x @test size(uvTN) == (2,) @test size(yvTN) == (2,) - uvTN, yvTN = (@test_logs min_level=Logging.Warn taylorinteg( + solTN = (@test_logs min_level=Logging.Warn taylorinteg( g, y0TN, u0, 10/0.3, _order, _abstol)) - uv, yv = (@test_logs min_level=Logging.Warn taylorinteg( + uvTN = solTN.t; yvTN = solTN.x + sol = (@test_logs min_level=Logging.Warn taylorinteg( g, y0, u0, 10/0.3, _order, _abstol)) + uv = sol.t; yv = sol.x exactsol_g(u, y0, u0) = y0*exp(0.3(u-u0)) δsol_g = exactsol_g(uvTN[end], y0TN, u0)-yvTN[end] δcoeffs_g = map(y->y[1], map(x->x.coeffs, δsol_g.coeffs)) @@ -76,8 +83,9 @@ import Logging: Warn disp = 0.001*rand() #a small, random displacement y0_disp = y0+disp dv = map(x->[disp], uvTN) #a vector of identical displacements - uv_disp, yv_disp = (@test_logs min_level=Logging.Warn taylorinteg( + sol_disp = (@test_logs min_level=Logging.Warn taylorinteg( g, y0_disp, u0, 10/0.3, _order, _abstol)) + uv_disp = sol_disp.t; yv_disp = sol_disp.x yvTN_disp = evaluate.(yvTN, dv) @test norm(exactsol_g.(uvTN, y0_disp, u0)-yvTN_disp, Inf) < 1E-9 #analytical vs jet transport @test norm(y0_disp-yvTN[1]([disp]), Inf) < 1E-9 @@ -92,14 +100,17 @@ import Logging: Warn x0T1 = x0 + p #jet transport initial condition t0 = 0.0 tmax = 0.3 - tvT1, xvT1 = (@test_logs (Warn, max_iters_reached()) taylorinteg( + solT1 = (@test_logs (Warn, max_iters_reached()) taylorinteg( f, x0T1, t0, tmax, _order, _abstol, maxsteps=1)) + tvT1 = solT1.t; xvT1 = solT1.x @test size(tvT1) == (2,) @test size(xvT1) == (2,) - tvT1, xvT1 = (@test_logs min_level=Logging.Warn taylorinteg( + solT1 = (@test_logs min_level=Logging.Warn taylorinteg( f, x0T1, t0, tmax, _order, _abstol)) - tv, xv = (@test_logs min_level=Logging.Warn taylorinteg( + tvT1 = solT1.t; xvT1 = solT1.x + sol = (@test_logs min_level=Logging.Warn taylorinteg( f, x0, t0, tmax, _order, _abstol)) + tv = sol.t; xv = sol.x exactsol(t, x0, t0) = x0/(1.0-x0*(t-t0)) #the analytical solution δsol = exactsol(tvT1[end], x0T1, t0)-xvT1[end] #analytical vs jet transport diff at end of integration δcoeffs = δsol.coeffs @@ -112,8 +123,9 @@ import Logging: Warn for i in 1:5 disp = 0.001*rand() #a small, random displacement x0_disp = x0+disp - tv_disp, xv_disp = (@test_logs min_level=Logging.Warn taylorinteg( + sol_disp = (@test_logs min_level=Logging.Warn taylorinteg( f, x0_disp, t0, tmax, _order, _abstol)) + tv_disp = sol_disp.t; xv_disp = sol_disp.x xvT1_disp = evaluate.(xvT1, disp) @test norm(exactsol.(tvT1, x0_disp, t0)-xvT1_disp, Inf) < 1E-12 #analytical vs jet transport @test norm(x0_disp-evaluate(xvT1[1], disp), Inf) < 1E-12 @@ -124,15 +136,17 @@ import Logging: Warn y0 = 1.0 #"nominal" initial condition u0 = 0.0 #initial time y0T1 = y0 + p #jet transport initial condition - uvT1, yvT1 = (@test_logs (Warn, max_iters_reached()) taylorinteg( + solT1 = (@test_logs (Warn, max_iters_reached()) taylorinteg( g, y0T1, t0, 10/0.3, _order, _abstol, maxsteps=1)) #warmup lap # test maxsteps break - @test size(uvT1) == (2,) - @test size(yvT1) == (2,) - uvT1, yvT1 = (@test_logs min_level=Logging.Warn taylorinteg( + @test size(solT1.t) == (2,) + @test size(solT1.x) == (2,) + solT1 = (@test_logs min_level=Logging.Warn taylorinteg( g, y0T1, t0, 10/0.3, _order, _abstol)) #Taylor1 jet transport integration - uv, yv = (@test_logs min_level=Logging.Warn taylorinteg( + uvT1 = solT1.t; yvT1 = solT1.x + sol = (@test_logs min_level=Logging.Warn taylorinteg( g, y0, t0, 10/0.3, _order, _abstol)) #reference integration + uv = sol.t; yv = sol.x exactsol_g(u, y0, u0) = y0*exp(0.3(u-u0)) #the analytical solution δsol_g = exactsol_g(uvT1[end], y0T1, t0)-yvT1[end] #analytical vs jet transport diff at end of integration δcoeffs_g = δsol_g.coeffs @@ -145,8 +159,9 @@ import Logging: Warn for i in 1:5 disp = 0.001*rand() #a small, random displacement y0_disp = y0+disp - uv_disp, yv_disp = (@test_logs min_level=Logging.Warn taylorinteg( + sol_disp = (@test_logs min_level=Logging.Warn taylorinteg( g, y0_disp, u0, 10/0.3, _order, _abstol)) + uv_disp = sol_disp.t; yv_disp = sol_disp.x yvT1_disp = evaluate.(yvT1, disp) @test norm(exactsol_g.(uvT1, y0_disp, u0)-yvT1_disp, Inf) < 1E-9 #analytical vs jet transport @test norm(y0_disp-evaluate(yvT1[1], disp), Inf) < 1E-9 @@ -162,25 +177,25 @@ import Logging: Warn tv = 0.0:0.05:0.33 xvT1 = (@test_logs (Warn, max_iters_reached()) taylorinteg( f, x0T1, tv, _order, _abstol, maxsteps=1)) - @test size(xvT1) == (7,) + @test size(xvT1.x) == (7,) ta = vec(tv) xvT1 = (@test_logs (Warn, max_iters_reached()) taylorinteg( f, x0T1, ta, _order, _abstol, maxsteps=1)) - @test size(xvT1) == (7,) + @test size(xvT1.x) == (7,) xvT1 = (@test_logs min_level=Logging.Warn taylorinteg( f, x0T1, tv, _order, _abstol)) xvT11 = (@test_logs min_level=Logging.Warn taylorinteg( f, x0T1, ta, _order, _abstol)) - @test xvT11 == xvT1 + @test xvT11.x == xvT1.x xv = (@test_logs min_level=Logging.Warn taylorinteg( f, x0, tv, _order, _abstol)) exactsol(t, x0, t0) = x0/(1.0-x0*(t-t0)) #the analytical solution - δsol = exactsol(tv[end], x0T1, tv[1])-xvT1[end] + δsol = exactsol(tv[end], x0T1, tv[1])-xvT1.x[end] δcoeffs = δsol.coeffs - @test length(tv) == length(xvT1) + @test length(tv) == length(xvT1.x) @test isapprox(δcoeffs, zeros(6), atol=1e-10, rtol=0) xv_analytical = exactsol.(tv, x0, tv[1]) - xvT1_0 = xvT1.() + xvT1_0 = xvT1.x.() @test isapprox(xv_analytical, xvT1_0, atol=1e-10, rtol=0) for i in 1:5 disp = 0.001*rand() #a small, random displacement @@ -188,10 +203,10 @@ import Logging: Warn f, x0+disp, tv, _order, _abstol)) xv_disp2 = (@test_logs min_level=Logging.Warn taylorinteg( f, x0+disp, ta, _order, _abstol)) - @test xv_disp == xv_disp2 - xvT1_disp = xvT1.(disp) + @test xv_disp.x == xv_disp2.x + xvT1_disp = xvT1.x.(disp) @test norm(exactsol.(tv, x0+disp, tv[1])-xvT1_disp, Inf) < 1E-12 #analytical vs jet transport - @test norm(xv_disp-xvT1_disp, Inf) < 1E-12 # integration vs jet transport + @test norm(xv_disp.x-xvT1_disp, Inf) < 1E-12 # integration vs jet transport end y0 = 1.0 #"nominal" initial condition @@ -199,26 +214,26 @@ import Logging: Warn uv = 0.0:1/0.3:10/0.3 yvT1 = (@test_logs (Warn, max_iters_reached()) taylorinteg( g, y0T1, uv, _order, _abstol, maxsteps=1)) - @test size(yvT1) == (11,) + @test size(yvT1.x) == (11,) yvT1 = (@test_logs min_level=Logging.Warn taylorinteg( g, y0T1, uv, _order, _abstol)) yv = (@test_logs min_level=Logging.Warn taylorinteg( g, y0, uv, _order, _abstol)) exactsol_g(u, y0, u0) = y0*exp(0.3(u-u0)) - δsol_g = exactsol_g(uv[end], y0T1, uv[1])-yvT1[end] + δsol_g = exactsol_g(uv[end], y0T1, uv[1])-yvT1.x[end] δcoeffs_g = δsol_g.coeffs - @test length(uv) == length(yvT1) + @test length(uv) == length(yvT1.x) @test isapprox(δcoeffs_g, zeros(6), atol=1e-10, rtol=0) yv_analytical = exactsol_g.(uv, y0, uv[1]) - yvT1_0 = evaluate.(yvT1) + yvT1_0 = evaluate.(yvT1.x) @test isapprox(yv_analytical, yvT1_0, atol=1e-10, rtol=0) for i in 1:5 disp = 0.001*rand() #a small, random displacement - yv_disp = (@test_logs min_level=Logging.Warn taylorinteg( + sol_disp = (@test_logs min_level=Logging.Warn taylorinteg( g, y0+disp, uv, _order, _abstol)) - yvT1_disp = evaluate.(yvT1, disp) + yvT1_disp = evaluate.(yvT1.x, disp) @test norm(exactsol_g.(uv, y0+disp, uv[1])-yvT1_disp, Inf) < 1E-9 #analytical vs jet transport - @test norm(yv_disp-yvT1_disp, Inf) < 1E-9 # integration vs jet transport + @test norm(sol_disp.x-yvT1_disp, Inf) < 1E-9 # integration vs jet transport end end @@ -229,16 +244,17 @@ import Logging: Warn tv = 0.0:0.05:0.33 xvTN = (@test_logs (Warn, max_iters_reached()) taylorinteg( f, x0TN, tv, _order, _abstol, maxsteps=1)) - @test size(xvTN) == (7,) + @test size(xvTN.x) == (7,) ta = vec(tv) xvTN = (@test_logs (Warn, max_iters_reached()) taylorinteg( f, x0TN, ta, _order, _abstol, maxsteps=1)) - @test size(xvTN) == (7,) - xvTN = (@test_logs min_level=Logging.Warn taylorinteg( + @test size(xvTN.x) == (7,) + solTN = (@test_logs min_level=Logging.Warn taylorinteg( f, x0TN, tv, _order, _abstol)) - xvTN_ = (@test_logs min_level=Logging.Warn taylorinteg( + xvTN = solTN.x + solTN_ = (@test_logs min_level=Logging.Warn taylorinteg( f, x0TN, ta, _order, _abstol)) - @test xvTN == xvTN_ + @test xvTN == solTN_.x xv = (@test_logs min_level=Logging.Warn taylorinteg( f, x0, tv, _order, _abstol)) exactsol(t, x0, t0) = x0/(1.0-x0*(t-t0)) #the analytical solution @@ -252,24 +268,25 @@ import Logging: Warn for i in 1:5 disp = 0.001*rand() #a small, random displacement dv = map(x->[disp], tv) #a vector of identical displacements - xv_disp = (@test_logs min_level=Logging.Warn taylorinteg( + sol_disp = (@test_logs min_level=Logging.Warn taylorinteg( f, x0+disp, tv, _order, _abstol)) - xv_disp_ = (@test_logs min_level=Logging.Warn taylorinteg( + sol_disp_ = (@test_logs min_level=Logging.Warn taylorinteg( f, x0+disp, ta, _order, _abstol)) - @test xv_disp == xv_disp_ + @test sol_disp.x == sol_disp_.x xvTN_disp = evaluate.(xvTN, dv) @test norm(exactsol.(tv, x0+disp, tv[1])-xvTN_disp, Inf) < 1E-12 #analytical vs jet transport - @test norm(xv_disp-xvTN_disp, Inf) < 1E-12 # integration vs jet transport + @test norm(sol_disp.x-xvTN_disp, Inf) < 1E-12 # integration vs jet transport end y0 = 1.0 #"nominal" initial condition y0TN = y0 + p[1] #jet transport initial condition uv = 0.0:1/0.3:10/0.3 - yvTN = (@test_logs (Warn, max_iters_reached()) taylorinteg( + solTN = (@test_logs (Warn, max_iters_reached()) taylorinteg( g, y0TN, uv, _order, _abstol, maxsteps=1)) - @test size(yvTN) == (11,) - yvTN = (@test_logs min_level=Logging.Warn taylorinteg( + @test size(solTN.x) == (11,) + solTN = (@test_logs min_level=Logging.Warn taylorinteg( g, y0TN, uv, _order, _abstol)) + yvTN = solTN.x yv = (@test_logs min_level=Logging.Warn taylorinteg( g, y0, uv, _order, _abstol)) exactsol_g(u, y0, u0) = y0*exp(0.3(u-u0)) @@ -283,11 +300,11 @@ import Logging: Warn for i in 1:5 disp = 0.001*rand() #a small, random displacement dv = map(x->[disp], uv) #a vector of identical displacements - yv_disp = (@test_logs min_level=Logging.Warn taylorinteg( + sol_disp = (@test_logs min_level=Logging.Warn taylorinteg( g, y0+disp, uv, _order, _abstol)) yvTN_disp = evaluate.(yvTN, dv) @test norm(exactsol_g.(uv, y0+disp, uv[1])-yvTN_disp, Inf) < 1E-9 #analytical vs jet transport - @test norm(yv_disp-yvTN_disp, Inf) < 1E-9 # integration vs jet transport + @test norm(sol_disp.x-yvTN_disp, Inf) < 1E-9 # integration vs jet transport end end @@ -304,12 +321,14 @@ import Logging: Warn ω0 = 1.0 x0 = [0.0,ω0,ω0] x0T1 = x0+[0t,t,t] - tv1, xv1 = (@test_logs (Warn, max_iters_reached()) taylorinteg( + sol1 = (@test_logs (Warn, max_iters_reached()) taylorinteg( harmosc!, x0T1, 0.0, 100pi, _order, _abstol, maxsteps=1)) + tv1 = sol1.t; xv1 = sol1.x @test size(tv1) == (2,) @test size(xv1) == (2, 3) - tv1, xv1 = (@test_logs min_level=Logging.Warn taylorinteg( + sol1 = (@test_logs min_level=Logging.Warn taylorinteg( harmosc!, x0T1, 0.0, 100pi, _order, _abstol, maxsteps=2000)) + tv1 = sol1.t; xv1 = sol1.x y0 = evaluate.(xv1) x1(t,δω) = sin((ω0+δω)*t) x2(t,δω) = (ω0+δω)*cos((ω0+δω)*t) @@ -318,8 +337,9 @@ import Logging: Warn for i in 1:5 δω=0.001*rand() x0_disp = x0+[0.0,δω,δω] - tv, xv = (@test_logs min_level=Logging.Warn taylorinteg( + sol = (@test_logs min_level=Logging.Warn taylorinteg( harmosc!, x0_disp, 0.0, 100pi, _order, _abstol, maxsteps=2000)) + tv = sol.t; xv = sol.x y1 = evaluate.(xv1, δω) @test norm(y1[:,1]-x1.(tv1,δω),Inf) < 1E-11 @test norm(y1[:,2]-x2.(tv1,δω),Inf) < 1E-11 @@ -335,19 +355,23 @@ import Logging: Warn x0 = [0.0,ω0,ω0] x0T1 = x0+[0t,t,t] tv = 0.0:0.25*(2pi):100pi - xv1 = (@test_logs (Warn, max_iters_reached()) taylorinteg( + sol1 = (@test_logs (Warn, max_iters_reached()) taylorinteg( harmosc!, x0T1, tv, _order, _abstol, maxsteps=1)) + xv1 = sol1.x @test length(tv) == 201 @test size(xv1) == (201, 3) ta = vec(tv) - xv1 = (@test_logs (Warn, max_iters_reached()) taylorinteg( + sol1 = (@test_logs (Warn, max_iters_reached()) taylorinteg( harmosc!, x0T1, ta, _order, _abstol, maxsteps=1)) + xv1 = sol1.x @test length(ta) == 201 @test size(xv1) == (201, 3) - xv1 = (@test_logs min_level=Logging.Warn taylorinteg( + sol1 = (@test_logs min_level=Logging.Warn taylorinteg( harmosc!, x0T1, tv, _order, _abstol, maxsteps=2000)) - xv1_ = (@test_logs min_level=Logging.Warn taylorinteg( + xv1 = sol1.x + sol1_ = (@test_logs min_level=Logging.Warn taylorinteg( harmosc!, x0T1, ta, _order, _abstol, maxsteps=2000)) + xv1_ = sol1_.x @test xv1 == xv1_ y0 = evaluate.(xv1) x1(t,δω) = sin((ω0+δω)*t) @@ -357,10 +381,12 @@ import Logging: Warn for i in 1:5 δω=0.001*rand() x0_disp = x0+[0.0,δω,δω] - xv = (@test_logs min_level=Logging.Warn taylorinteg( + sol = (@test_logs min_level=Logging.Warn taylorinteg( harmosc!, x0_disp, tv, _order, _abstol, maxsteps=2000)) - xv_ = (@test_logs min_level=Logging.Warn taylorinteg( + xv = sol.x + sol_ = (@test_logs min_level=Logging.Warn taylorinteg( harmosc!, x0_disp, ta, _order, _abstol, maxsteps=2000)) + xv_ = sol_.x @test xv == xv_ y1 = evaluate.(xv1, δω) @test norm(y1[:,1]-x1.(tv,δω), Inf) < 1E-11 @@ -381,14 +407,17 @@ import Logging: Warn p = set_variables("ξ", numvars=2, order=5) x0 = [-1.0,0.45] x0TN = x0 + p - tvTN, xvTN = (@test_logs (Warn, max_iters_reached()) taylorinteg( + solTN = (@test_logs (Warn, max_iters_reached()) taylorinteg( harmosc!, x0TN, 0.0, 10pi, _order, _abstol, maxsteps=1)) + tvTN = solTN.t; xvTN = solTN.x @test length(tvTN) == 2 @test size(xvTN) == (length(tvTN), length(x0TN)) - tvTN, xvTN = (@test_logs min_level=Logging.Warn taylorinteg( + solTN = (@test_logs min_level=Logging.Warn taylorinteg( harmosc!, x0TN, 0.0, 10pi, _order, _abstol)) - tv , xv = (@test_logs min_level=Logging.Warn taylorinteg( + tvTN = solTN.t; xvTN = solTN.x + sol = (@test_logs min_level=Logging.Warn taylorinteg( harmosc!, x0 , 0.0, 10pi, _order, _abstol)) + tv = sol.t; xv = sol.x x_analyticsol(t,x0,p0) = p0*sin(t)+x0*cos(t) p_analyticsol(t,x0,p0) = p0*cos(t)-x0*sin(t) x_δsol = x_analyticsol(tvTN[end], x0TN[1], x0TN[2])-xvTN[end,1] @@ -428,17 +457,20 @@ import Logging: Warn #the time range tr = t0:integstep:tmax; #xv is the solution vector representing the propagation of the initial condition q0 propagated until time T - xv = (@test_logs min_level=Logging.Warn taylorinteg( + sol = (@test_logs min_level=Logging.Warn taylorinteg( pendulum!, q0, tr, _order, _abstol, maxsteps=100)) + xv = sol.x #xvTN is the solution vector representing the propagation of the initial condition q0 plus variations (ξ₁,ξ₂) propagated until time T #note that q0 is a Vector{Float64}, but q0TN is a Vector{TaylorN{Float64}} #but apart from that difference, we're calling `taylorinteg` essentially with the same parameters! #thus, jet transport is reduced to a beautiful application of Julia's multiple dispatch! - xvTN = (@test_logs (Warn, max_iters_reached()) taylorinteg( + solTN = (@test_logs (Warn, max_iters_reached()) taylorinteg( pendulum!, q0TN, tr, _order, _abstol, maxsteps=1)) + xvTN = solTN.x @test size(xvTN) == (5,2) - xvTN = (@test_logs min_level=Logging.Warn taylorinteg( + solTN = (@test_logs min_level=Logging.Warn taylorinteg( pendulum!, q0TN, tr, _order, _abstol, maxsteps=100)) + xvTN = solTN.x xvTN_0 = map( x->evaluate(x, [0.0, 0.0]), xvTN ) # the jet evaluated at the nominal solution @@ -446,8 +478,9 @@ import Logging: Warn @test isapprox(xv, xvTN_0) #nominal solution must coincide with jet evaluated at ξ=(0,0) #testing another jet transport method: - tv, xv = (@test_logs min_level=Logging.Warn taylorinteg( + sol = (@test_logs min_level=Logging.Warn taylorinteg( pendulum!, q0, t0, tmax, _order, _abstol, maxsteps=100)) + tv = sol.t; xv = sol.x @test isapprox(xv[1,:], xvTN_0[1,:]) #nominal solution must coincide with jet evaluated at ξ=(0,0) at initial time @test isapprox(xv[end,:], xvTN_0[end,:]) #nominal solution must coincide with jet evaluated at ξ=(0,0) at final time @@ -455,8 +488,9 @@ import Logging: Warn disp = 0.0001 #works even with 0.001, but we're giving it some margin #compare the jet solution evaluated at various variation vectors ξ, wrt to full solution, at each time evaluation point - xv_disp = (@test_logs (Warn, max_iters_reached()) taylorinteg( + sol_disp = (@test_logs (Warn, max_iters_reached()) taylorinteg( pendulum!, q0+[disp,0.0], tr, _order, _abstol, maxsteps=1)) + xv_disp = sol_disp.x @test size(xv_disp) == (5,2) for i in 1:10 # generate a random angle @@ -464,12 +498,12 @@ import Logging: Warn # generate a random point on a circumference of radius disp ξ = disp*[cos(ϕ), sin(ϕ)] #propagate in time full solution with initial condition q0+ξ - xv_disp = (@test_logs min_level=Logging.Warn taylorinteg( + sol_disp = (@test_logs min_level=Logging.Warn taylorinteg( pendulum!, q0+ξ, tr, _order, _abstol, maxsteps=100)) #evaluate jet at q0+ξ xvTN_disp = map( x->evaluate(x, ξ), xvTN ) #the propagated exact solution at q0+ξ should be approximately equal to the propagated jet solution evaluated at the same point: - @test isapprox( xvTN_disp, xv_disp ) + @test isapprox( xvTN_disp, sol_disp.x ) end end end diff --git a/test/lyapunov.jl b/test/lyapunov.jl index 8c8c534f0..b308c482f 100644 --- a/test/lyapunov.jl +++ b/test/lyapunov.jl @@ -118,15 +118,17 @@ import Logging: Warn xi = set_variables("δ", order=1, numvars=length(q0)) #Test lyap_taylorinteg with autodiff-computed Jacobian, maxsteps=5 - tv, xv, λv = (@test_logs (Warn, max_iters_reached()) lyap_taylorinteg( + sol = (@test_logs (Warn, max_iters_reached()) lyap_taylorinteg( lorenz!, q0, t0, tmax, _order, _abstol, nothing; maxsteps=5)) + tv, xv, λv = sol.t, sol.x, sol.λ @test size(tv) == (6,) @test size(xv) == (6,3) @test size(λv) == (6,3) #Test lyap_taylorinteg with user-defined Jacobian, maxsteps=5 - tv2, xv2, λv2 = (@test_logs (Warn, max_iters_reached()) lyap_taylorinteg( + sol2 = (@test_logs (Warn, max_iters_reached()) lyap_taylorinteg( lorenz!, q0, t0, tmax, _order, _abstol, nothing, lorenz_jac!; maxsteps=5)) + tv2, xv2, λv2 = sol2.t, sol2.x, sol2.λ @test size(tv2) == (6,) @test size(xv2) == (6,3) @test size(λv2) == (6,3) @@ -135,8 +137,9 @@ import Logging: Warn @test λv == λv2 #Test lyap_taylorinteg with autodiff-computed Jacobian, maxsteps=2000 - tv, xv, λv = (@test_logs min_level=Logging.Warn lyap_taylorinteg( + sol = (@test_logs min_level=Logging.Warn lyap_taylorinteg( lorenz!, q0, t0, tmax, _order, _abstol; maxsteps=2000)) + tv, xv, λv = sol.t, sol.x, sol.λ @test xv[1,:] == q0 @test tv[1] == t0 @test size(xv) == size(λv) @@ -147,18 +150,21 @@ import Logging: Warn @test isapprox(λv[end,2], -0.00830, rtol=mytol, atol=mytol) @test isapprox(λv[end,3], -22.46336, rtol=mytol, atol=mytol) # Check integration consistency (orbit should not depend on variational eqs) - t_, x_ = taylorinteg(lorenz!, q0, t0, tmax, _order, _abstol; maxsteps=2000) + sol_ = taylorinteg(lorenz!, q0, t0, tmax, _order, _abstol; maxsteps=2000) + t_, x_ = sol_.t, sol_.x @test t_ == tv @test x_ == xv # test backward integration - tvb, xvb, λvb = (@test_logs (Warn, max_iters_reached()) lyap_taylorinteg( + solb = (@test_logs (Warn, max_iters_reached()) lyap_taylorinteg( lorenz!, q0, t0, -tmax, _order, _abstol; maxsteps=2000)) + tvb, xvb, λvb = solb.t, solb.x, solb.λ @test isapprox(sum(λvb[1,:]), lorenztr) == false @test isapprox(sum(λvb[end,:]), lorenztr) #Test lyap_taylorinteg with user-defined Jacobian, maxsteps=2000 - tv2, xv2, λv2 = (@test_logs min_level=Logging.Warn lyap_taylorinteg( + sol2 = (@test_logs min_level=Logging.Warn lyap_taylorinteg( lorenz!, q0, t0, tmax, _order, _abstol, nothing, lorenz_jac!; maxsteps=2000)) + tv2, xv2, λv2 = sol2.t, sol2.x, sol2.λ @test xv2[1,:] == q0 @test tv2[1] == t0 @test size(xv2) == size(λv2) @@ -175,8 +181,9 @@ import Logging: Warn @test t_ == tv2 @test x_ == xv2 # test backward integration - tv2b, xv2b, λv2b = (@test_logs (Warn, max_iters_reached()) lyap_taylorinteg( + sol2b = (@test_logs (Warn, max_iters_reached()) lyap_taylorinteg( lorenz!, q0, t0, -tmax, _order, _abstol, nothing, lorenz_jac!; maxsteps=2000)) + tv2b, xv2b, λv2b = sol2b.t, sol2b.x, sol2b.λ @test isapprox(sum(λv2b[1,:]), lorenztr) == false @test isapprox(sum(λv2b[end,:]), lorenztr) end @@ -195,32 +202,40 @@ import Logging: Warn @test lorenztr == -(1+σ+β) trange = t0:1.0:tmax - xw, λw = (@test_logs (Warn, max_iters_reached()) lyap_taylorinteg( + solw = (@test_logs (Warn, max_iters_reached()) lyap_taylorinteg( lorenz!, q0, trange, _order, _abstol; maxsteps=5)) + xw, λw = solw.x, solw.λ + @test trange == solw.t @test xw[1,:] == q0 @test size(xw) == (length(trange), 3) @test λw[1,:] == zeros(3) @test size(λw) == size(xw) @test all(isnan.(xw[2:end,:])) @test all(isnan.(λw[2:end,:])) - xw2, λw2 = (@test_logs (Warn, max_iters_reached()) lyap_taylorinteg( + solw2 = (@test_logs (Warn, max_iters_reached()) lyap_taylorinteg( lorenz!, q0, vec(trange), _order, _abstol; maxsteps=5)) + xw2, λw2 = solw2.x, solw2.λ + @test trange == solw2.t @test xw[1, :] == xw2[1, :] @test all(isnan.(xw2[2:end,:])) @test λw2[1, :] == zeros(3) @test all(isnan.(λw2[2:end,:])) @test size(xw2) == (length(trange), 3) @test size(λw2) == size(xw2) - xw_, λw_ = (@test_logs (Warn, max_iters_reached()) lyap_taylorinteg( + solw_ = (@test_logs (Warn, max_iters_reached()) lyap_taylorinteg( lorenz!, q0, trange, _order, _abstol, nothing, lorenz_jac!; maxsteps=5)) + xw_, λw_ = solw_.x, solw_.λ + @test trange == solw_.t @test xw_[1,:] == q0 @test all(isnan.(xw_[2:end,:])) @test size(xw_) == (length(trange), 3) @test size(λw_) == size(xw_) @test λw_[1, :] == zeros(3) @test all(isnan.(λw_[2:end,:])) - xw2_, λw2_ = (@test_logs (Warn, max_iters_reached()) lyap_taylorinteg( + solw2_ = (@test_logs (Warn, max_iters_reached()) lyap_taylorinteg( lorenz!, q0, vec(trange), _order, _abstol, nothing, lorenz_jac!; maxsteps=5)) + xw2_, λw2_ = solw2_.x, solw2_.λ + @test vec(trange) == solw2_.t @test xw_[1, :] == xw2_[1, :] @test all(isnan.(xw2_[2:end,:])) @test λw2_[1, :] == zeros(3) @@ -228,15 +243,18 @@ import Logging: Warn @test size(xw2_) == (length(trange), 3) @test size(λw2_) == size(xw2_) - xw, λw = (@test_logs min_level=Logging.Warn lyap_taylorinteg( + solw = (@test_logs min_level=Logging.Warn lyap_taylorinteg( lorenz!, q0, trange, _order, _abstol; maxsteps=2000)) + xw, λw = solw.x, solw.λ + @test trange == solw.t @test xw[1,:] == q0 @test size(xw) == (length(trange), length(q0)) @test size(λw) == (length(trange), length(q0)) @test isapprox(sum(λw[1,:]), lorenztr) == false @test isapprox(sum(λw[end,:]), lorenztr) - tz, xz, λz = (@test_logs min_level=Logging.Warn lyap_taylorinteg( + solz = (@test_logs min_level=Logging.Warn lyap_taylorinteg( lorenz!, q0, trange[1], trange[end], _order, _abstol; maxsteps=2000)) + tz, xz, λz = solz.t, solz.x, solz.λ @test λw[end,:] == λz[end,:] @test xw[end,:] == xz[end,:] @test tz[end] == trange[end] @@ -244,26 +262,31 @@ import Logging: Warn @test isapprox(λw[end,1], 1.47167, rtol=mytol, atol=mytol) @test isapprox(λw[end,2], -0.00831, rtol=mytol, atol=mytol) @test isapprox(λw[end,3], -22.46336, rtol=mytol, atol=mytol) - xw_, λw_ = (@test_logs min_level=Logging.Warn lyap_taylorinteg( + solw_ = (@test_logs min_level=Logging.Warn lyap_taylorinteg( lorenz!, q0, trange, _order, _abstol, nothing, lorenz_jac!; maxsteps=2000)) + xw_, λw_ = solw_.x, solw_.λ @test xw == xw_ @test λw == λw_ - tz_, xz_, λz_ = (@test_logs min_level=Logging.Warn lyap_taylorinteg( + solz_ = (@test_logs min_level=Logging.Warn lyap_taylorinteg( lorenz!, q0, trange[1], trange[end], _order, _abstol, nothing, lorenz_jac!; maxsteps=2000)) + tz_, xz_, λz_ = solz_.t, solz_.x, solz_.λ @test λw_[end,:] == λz_[end,:] @test xw_[end,:] == xz_[end,:] @test tz_[end] == trange[end] - xw2, λw2 = (@test_logs min_level=Logging.Warn lyap_taylorinteg( + solw2 = (@test_logs min_level=Logging.Warn lyap_taylorinteg( lorenz!, q0, vec(trange), _order, _abstol; maxsteps=2000)) + xw2, λw2 = solw2.x, solw2.λ + @test vec(trange) == solw2.t @test xw2 == xw @test λw2 == λw @test xw2[1,:] == q0 @test size(xw) == (length(trange), length(q0)) @test isapprox(sum(λw2[1,:]), lorenztr) == false @test isapprox(sum(λw2[end,:]), lorenztr) - tz2, xz2, λz2 = (@test_logs min_level=Logging.Warn lyap_taylorinteg( + solz2 = (@test_logs min_level=Logging.Warn lyap_taylorinteg( lorenz!, q0, trange[1], trange[end], _order, _abstol; maxsteps=2000)) + tz2, xz2, λz2 = solz2.t, solz2.x, solz2.λ @test λw2[end,:] == λz2[end,:] @test xw2[end,:] == xz2[end,:] @test tz2[end] == trange[end] @@ -271,18 +294,21 @@ import Logging: Warn @test isapprox(λw2[end,1], 1.47167, rtol=mytol, atol=mytol) @test isapprox(λw2[end,2], -0.00831, rtol=mytol, atol=mytol) @test isapprox(λw2[end,3], -22.46336, rtol=mytol, atol=mytol) - xw2_, λw2_ = (@test_logs min_level=Logging.Warn lyap_taylorinteg( + solw2_ = (@test_logs min_level=Logging.Warn lyap_taylorinteg( lorenz!, q0, vec(trange), _order, _abstol, nothing, lorenz_jac!; maxsteps=2000)) + xw2_, λw2_ = solw2_.x, solw2.λ @test xw2 == xw2_ @test λw2 == λw2_ - tz2_, xz2_, λz2_ = (@test_logs min_level=Logging.Warn lyap_taylorinteg( + solz2_ = (@test_logs min_level=Logging.Warn lyap_taylorinteg( lorenz!, q0, trange[1], trange[end], _order, _abstol, nothing, lorenz_jac!; maxsteps=2000)) + tz2_, xz2_, λz2_ = solz2_.t, solz2_.x, solz2.λ @test λw2_[end,:] == λz2_[end,:] @test xw2_[end,:] == xz2_[end,:] @test tz2_[end] == trange[end] # Check integration consistency (orbit should not depend on variational eqs) - x_ = taylorinteg(lorenz!, q0, vec(trange), _order, _abstol; maxsteps=2000) + sol_ = taylorinteg(lorenz!, q0, vec(trange), _order, _abstol; maxsteps=2000) + x_ = sol_.x @test x_ == xw2 end diff --git a/test/many_ode.jl b/test/many_ode.jl index 728a37c45..902d2f6d8 100644 --- a/test/many_ode.jl +++ b/test/many_ode.jl @@ -34,16 +34,20 @@ import Logging: Warn δt = (_abstol/q0T[1].coeffs[end-1])^inv(_order-1) @test TaylorIntegration.stepsize(q0T, _abstol) == δt - tv, xv = (@test_logs (Warn, max_iters_reached()) taylorinteg( - eqs_mov!, q0, 0.0, 0.5, _order, _abstol, nothing)) + sol = @test_logs((Warn, max_iters_reached()), + @inferred(TaylorSolution{Float64, Float64, 2}, + taylorinteg(eqs_mov!, q0, 0.0, 0.5, _order, _abstol, nothing))) + @test_throws ErrorException sol(0.25) + tv = sol.t; xv = sol.x @test length(tv) == 501 @test isa(xv, SubArray) @test xv[1,:] == q0 @test tv[end] < 1/3 trange = 0.0:1/8:1.0 - xv = (@test_logs (Warn, max_iters_reached()) taylorinteg( + sol = (@test_logs (Warn, max_iters_reached()) taylorinteg( eqs_mov!, q0, trange, _order, _abstol)) + xv = sol.x @test size(xv) == (9,2) @test q0 == [3.0, 1.0] @test typeof(xv) == Transpose{Float64, Array{Float64,2}} @@ -54,8 +58,9 @@ import Logging: Warn @test abs(xv[2,1] - 4.8) ≤ eps(4.8) tarray = vec(trange) - xv2 = (@test_logs (Warn, max_iters_reached()) taylorinteg( + sol2 = (@test_logs (Warn, max_iters_reached()) taylorinteg( eqs_mov!, q0, tarray, _order, _abstol, nothing)) + xv2 = sol2.x @test xv[1:3,:] == xv2[1:3,:] @test xv2[1:3,:] ≈ xv[1:3,:] atol=eps() rtol=0.0 @test size(xv2) == (9,2) @@ -68,8 +73,9 @@ import Logging: Warn @test abs(xv2[2,1] - 4.8) ≤ eps(4.8) # Output includes Taylor polynomial solution - tv, xv, psol = (@test_logs (Warn, max_iters_reached()) taylorinteg( - eqs_mov!, q0, 0, 0.5, _order, _abstol, Val(true), nothing, maxsteps=2)) + sol = (@test_logs (Warn, max_iters_reached()) taylorinteg( + eqs_mov!, q0, 0, 0.5, _order, _abstol, nothing, dense=true, maxsteps=2)) + tv = sol.t; xv = sol.x; psol = sol.p @test size(psol) == (2, 2) @test xv[1,:] == q0 @test xv[2,:] == evaluate.(psol[1, :], tv[2]-tv[1]) @@ -88,8 +94,9 @@ import Logging: Warn q0 = [3.0, 3.0] tmax = 0.3 - tv, xv = (@test_logs min_level=Logging.Warn taylorinteg( + sol = (@test_logs min_level=Logging.Warn taylorinteg( eqs_mov!, q0, 0, tmax, _order, _abstol, nothing)) + tv = sol.t; xv = sol.x @test length(tv) < 501 @test length(xv[:,1]) < 501 @test length(xv[:,2]) < 501 @@ -104,8 +111,9 @@ import Logging: Warn @test abs(xv[end,2]-exactsol(tv[end], xv[1,2])) < 5e-14 tmax = 0.33 - tv, xv = (@test_logs min_level=Logging.Warn taylorinteg( + sol = (@test_logs min_level=Logging.Warn taylorinteg( eqs_mov!, [3, 3], 0.0, tmax, _order, _abstol)) + tv = sol.t; xv = sol.x @test length(tv) < 501 @test length(xv[:,1]) < 501 @test length(xv[:,2]) < 501 @@ -120,10 +128,12 @@ import Logging: Warn @test abs(xv[end,2]-exactsol(tv[end], xv[1,2])) < 1.0e-11 # Output includes Taylor polynomial solution - tv, xv, psol = taylorinteg(eqs_mov!, q0, 0, tmax, _order, _abstol, Val(true), nothing) + sol = taylorinteg(eqs_mov!, q0, 0, tmax, _order, _abstol, nothing, dense=true) + tv = sol.t; xv = sol.x; psol = sol.p @test size(psol) == ( size(xv, 1)-1, size(xv, 2) ) @test psol[1,:] == fill(Taylor1([3.0^i for i in 1:_order+1]), 2) @test xv[end,1] == evaluate.(psol[end, 1], tv[end]-tv[end-1]) + @test xv[end,:] == sol(tv[end]) @test psol[:,1] == psol[:,2] end @@ -139,8 +149,9 @@ import Logging: Warn order = 25 x0 = [t0, 0.0] #initial conditions such that x(t)=sin(t) - tv, xv = (@test_logs min_level=Logging.Warn taylorinteg( + sol = (@test_logs min_level=Logging.Warn taylorinteg( f!, x0, t0, tmax, order, abstol)) + tv = sol.t; xv = sol.x @test length(tv) < 501 @test length(xv[:,1]) < 501 @test length(xv[:,2]) < 501 @@ -150,8 +161,9 @@ import Logging: Warn @test abs(sin(tmax)-xv[end,2]) < 1e-14 # Backward integration - tb, xb = (@test_logs min_level=Logging.Warn taylorinteg( - f!, [tmax, sin(tmax)], tmax, t0, order, abstol)) + solb = (@test_logs min_level=Logging.Warn taylorinteg( + f!, [tmax, sin(tmax)], tmax, t0, order, abstol, dense=true)) + tb = solb.t; xb = solb.x @test length(tb) < 501 @test length(xb[:,1]) < 501 @test length(xb[:,2]) < 501 @@ -159,20 +171,27 @@ import Logging: Warn @test xb[1,1:end] == [tmax, sin(tmax)] @test tb[end] == xb[end,1] @test abs(sin(t0)-xb[end,2]) < 5e-14 + @test solb(tmax) == [tmax, sin(tmax)] + @test solb(x0[1]) == xb[end,:] + tmidT = (tmax+t0)/2 + Taylor1(order) + @test solb(tmidT)[1] == tmidT + @test norm(solb(tmidT)[2] - sin(tmidT), Inf) < 1e-13 # Tests with a range, for comparison with backward integration tmax = 15*(2pi) Δt = (tmax-t0)/1024 tspan = t0:Δt:tmax - xv = (@test_logs min_level=Logging.Warn taylorinteg( + sol = (@test_logs min_level=Logging.Warn taylorinteg( f!, x0, tspan, order, abstol, nothing)) + xv = sol.x @test xv[1,1:end] == x0 @test tmax == xv[end,1] @test abs(sin(tmax)-xv[end,2]) < 1e-14 # Backward integration - xback = (@test_logs min_level=Logging.Warn taylorinteg( + solback = (@test_logs min_level=Logging.Warn taylorinteg( f!, xv[end, :], reverse(tspan), order, abstol, nothing)) + xback = solback.x @test xback[1,:] == xv[end, :] @test abs(xback[end,1]-x0[1]) < 5.0e-14 @test abs(xback[end,2]-x0[2]) < 5.0e-14 @@ -195,8 +214,9 @@ import Logging: Warn abstol = 1e-20 order = 10 x0 = [10.0, 0.0] #initial conditions such that x(t)=sin(t) - tv, xv = (@test_logs min_level=Logging.Warn taylorinteg( + sol = (@test_logs min_level=Logging.Warn taylorinteg( fallball!, x0, t0, tmax, order, abstol)) + tv = sol.t; xv = sol.x @test length(tv) < 501 @test length(xv[:,1]) < 501 @test exactsol.(tv, x0[1], x0[2]) ≈ xv[:,1] diff --git a/test/one_ode.jl b/test/one_ode.jl index 6536cb9fa..4761d79e7 100644 --- a/test/one_ode.jl +++ b/test/one_ode.jl @@ -25,53 +25,66 @@ import Logging: Warn δt = _abstol^inv(_order-1) @test TaylorIntegration.stepsize(x0T, _abstol) == δt - tv, xv = (@test_logs (Warn, max_iters_reached()) taylorinteg( + sol = (@test_logs (Warn, max_iters_reached()) taylorinteg( eqs_mov, 1, 0.0, 1.0, _order, _abstol)) + @test sol isa TaylorSolution{Float64, Float64, 1} + tv = sol.t + xv = sol.x + @test isnothing(sol.p) @test length(tv) == 501 @test length(xv) == 501 @test xv[1] == x0 @test tv[end] < 1.0 - tv, xv = (@test_logs (Warn, max_iters_reached()) taylorinteg( + sol = (@test_logs (Warn, max_iters_reached()) taylorinteg( eqs_mov, x0, 0.0, 1.0, _order, _abstol, nothing)) + tv = sol.t; xv = sol.x @test length(tv) == 501 @test length(xv) == 501 @test xv[1] == x0 @test tv[end] < 1.0 trange = 0.0:1/8:1.0 - xv = (@test_logs (Warn, max_iters_reached()) taylorinteg( + sol = (@test_logs (Warn, max_iters_reached()) taylorinteg( eqs_mov, 1, trange, _order, _abstol)) + xv = sol.x + tv = sol.t @test length(xv) == length(trange) - @test typeof(xv) == Array{typeof(x0),1} + @test xv isa Vector{typeof(x0)} @test xv[1] == x0 @test isnan(xv[end]) @test abs(xv[5] - 2.0) ≤ eps(2.0) - tvr, xvr = (@test_logs min_level=Logging.Warn taylorinteg( + solr = (@test_logs min_level=Logging.Warn taylorinteg( eqs_mov, x0, trange[1], trange[end-1], _order, _abstol)) + tvr = solr.t; xvr = solr.x @test tvr[1] == trange[1] @test tvr[end] == trange[end-1] @test xvr[1] == xv[1] @test xvr[end] == xv[end-1] trange = 0.0:1/8:1.0 - xv = (@test_logs (Warn, max_iters_reached()) taylorinteg( + sol = (@test_logs (Warn, max_iters_reached()) taylorinteg( eqs_mov, x0, trange, _order, _abstol, nothing)) + xv = sol.x + @test trange == sol.t @test length(xv) == length(trange) - @test typeof(xv) == Array{typeof(x0),1} + @test xv isa Vector{typeof(x0)} @test xv[1] == x0 @test isnan(xv[end]) @test abs(xv[5] - 2.0) ≤ eps(2.0) - tvr, xvr = (@test_logs min_level=Logging.Warn taylorinteg( + solr = (@test_logs min_level=Logging.Warn taylorinteg( eqs_mov, x0, trange[1], trange[end-1], _order, _abstol, nothing)) + tvr = solr.t; xvr = solr.x @test tvr[1] == trange[1] @test tvr[end] == trange[end-1] @test xvr[1] == xv[1] @test xvr[end] == xv[end-1] tarray = vec(trange) - xv2 = (@test_logs (Warn, max_iters_reached()) taylorinteg( + sol2 = (@test_logs (Warn, max_iters_reached()) taylorinteg( eqs_mov, x0, tarray, _order, _abstol)) + xv2 = sol2.x + @test trange == sol2.t @test xv[1:end-1] == xv2[1:end-1] @test xv2[1:end-1] ≈ xv[1:end-1] atol=eps() rtol=0.0 @test length(xv2) == length(tarray) @@ -81,12 +94,17 @@ import Logging: Warn @test abs(xv2[5] - 2.0) ≤ eps(2.0) # Output includes Taylor polynomial solution - tv, xv, psol = (@test_logs (Warn, max_iters_reached()) taylorinteg( - eqs_mov, x0, 0, 0.5, _order, _abstol, Val(true), maxsteps=2)) + sol = (@test_logs (Warn, max_iters_reached()) taylorinteg( + eqs_mov, x0, 0, 0.5, _order, _abstol, dense=true, maxsteps=2)) + tv = sol.t; xv = sol.x + psol = sol.p @test length(psol) == 2 @test xv[1] == x0 @test psol[1] == Taylor1(ones(_order+1)) @test xv[2] == evaluate(psol[1], tv[2]-tv[1]) + @test xv[3] == sol(tv[3]) + @test psol[1] == sol(tv[1]+Taylor1(_order)) + @test psol[2] == sol(tv[2]+Taylor1(_order)) end @testset "Tests: dot{x}=x^2, x(0) = 3; nsteps <= maxsteps" begin @@ -102,8 +120,9 @@ import Logging: Warn δt = (_abstol/x0T.coeffs[end-1])^inv(_order-1) @test TaylorIntegration.stepsize(x0T, _abstol) == δt - tv, xv = (@test_logs min_level=Logging.Warn taylorinteg( + sol = (@test_logs min_level=Logging.Warn taylorinteg( eqs_mov, x0, 0, tmax, _order, _abstol)) + tv = sol.t; xv = sol.x @test length(tv) < 501 @test length(xv) < 501 @test length(tv) == 14 @@ -114,8 +133,9 @@ import Logging: Warn @test abs(xv[end]-exactsol(tv[end], xv[1])) < 5e-14 tmax = 0.33 - tv, xv = (@test_logs min_level=Logging.Warn taylorinteg( + sol = (@test_logs min_level=Logging.Warn taylorinteg( eqs_mov, x0, t0, tmax, _order, _abstol)) + tv = sol.t; xv = sol.x @test length(tv) < 501 @test length(xv) < 501 @test length(tv) == 28 @@ -134,8 +154,9 @@ import Logging: Warn order = 25 x0 = 0.0 #initial conditions such that x(t)=sin(t) - tv, xv = (@test_logs min_level=Logging.Warn taylorinteg( + sol = (@test_logs min_level=Logging.Warn taylorinteg( f, x0, t0, tmax, order, abstol)) + tv = sol.t; xv = sol.x @test length(tv) < 501 @test length(xv) < 501 @test xv[1] == x0 @@ -143,8 +164,9 @@ import Logging: Warn @test abs(sin(tmax)-xv[end]) < 1e-14 # Backward integration - tb, xb = (@test_logs min_level=Logging.Warn taylorinteg( + solb = (@test_logs min_level=Logging.Warn taylorinteg( f, sin(tmax), tmax, t0, order, abstol)) + tb = solb.t; xb = solb.x @test length(tb) < 501 @test length(xb) < 501 @test xb[1] == sin(tmax) @@ -155,14 +177,16 @@ import Logging: Warn tmax = 15*(2pi) Δt = (tmax-t0)/1024 tspan = t0:Δt:tmax - xv = (@test_logs min_level=Logging.Warn taylorinteg( + sol = (@test_logs min_level=Logging.Warn taylorinteg( f, x0, tspan, order, abstol)) + xv = sol.x @test xv[1] == x0 @test abs(sin(tmax)-xv[end]) < 1e-14 # Backward integration - xback = (@test_logs min_level=Logging.Warn taylorinteg( + solback = (@test_logs min_level=Logging.Warn taylorinteg( f, xv[end], reverse(tspan), order, abstol)) + xback = solback.x @test xback[1] == xv[end] @test abs(sin(t0)-xback[end]) < 5e-14 @test norm(xv[:]-xback[end:-1:1], Inf) < 5.0e-14 diff --git a/test/rootfinding.jl b/test/rootfinding.jl index 80cf46fa6..d8c4ece41 100644 --- a/test/rootfinding.jl +++ b/test/rootfinding.jl @@ -32,8 +32,9 @@ import Logging: Warn x01 = x0 + [ξ, ξ] #warm-up lap and preliminary tests - tv, xv, tvS, xvS, gvS = (@test_logs (Warn, max_iters_reached()) taylorinteg( + sol = (@test_logs (Warn, max_iters_reached()) taylorinteg( pendulum!, g, x0, t0, Tend, _order, _abstol, maxsteps=1)) + tv, xv, tvS, xvS, gvS = sol.t, sol.x, sol.tevents, sol.xevents, sol.gresids @test size(tv) == (2,) @test tv[1] == t0 @test size(xv) == (2,2) @@ -42,8 +43,9 @@ import Logging: Warn @test size(xvS) == (0,2) @test size(gvS) == (0,) - tvN, xvN, tvSN, xvSN, gvSN = (@test_logs (Warn, max_iters_reached()) taylorinteg( + solN = (@test_logs (Warn, max_iters_reached()) taylorinteg( pendulum!, g, x0N, t0, 3Tend, _order, _abstol, maxsteps=1)) + tvN, xvN, tvSN, xvSN, gvSN = solN.t, solN.x, solN.tevents, solN.xevents, solN.gresids @test eltype(tvN) == Float64 @test eltype(xvN) == TaylorN{Float64} @test eltype(tvSN) == TaylorN{Float64} @@ -57,8 +59,9 @@ import Logging: Warn @test size(xvSN) == (0,2) @test size(gvSN) == (0,) - tv1, xv1, tvS1, xvS1, gvS1 = (@test_logs (Warn, max_iters_reached()) taylorinteg( + sol1 = (@test_logs (Warn, max_iters_reached()) taylorinteg( pendulum!, g, x01, t0, 3Tend, _order, _abstol, maxsteps=1)) + tv1, xv1, tvS1, xvS1, gvS1 = sol1.t, sol1.x, sol1.tevents, sol1.xevents, sol1.gresids @test eltype(tv1) == Float64 @test eltype(xv1) == Taylor1{Float64} @test eltype(tvS1) == Taylor1{Float64} @@ -73,8 +76,9 @@ import Logging: Warn @test size(gvS1) == (0,) #testing 0-th order root-finding - tv, xv, tvS, xvS, gvS = (@test_logs min_level=Logging.Warn taylorinteg( + sol = (@test_logs min_level=Logging.Warn taylorinteg( pendulum!, g, x0, t0, 3Tend, _order, _abstol, maxsteps=1000)) + tv, xv, tvS, xvS, gvS = sol.t, sol.x, sol.tevents, sol.xevents, sol.gresids @test tv[1] == t0 @test xv[1,:] == x0 @test size(tvS) == (5,) @@ -82,8 +86,9 @@ import Logging: Warn @test norm(gvS,Inf) < eps() #testing 0-th order root-finding with dense output - tv, xv, psol, tvS, xvS, gvS = (@test_logs min_level=Logging.Warn taylorinteg( - pendulum!, g, x0, t0, 3Tend, _order, _abstol, Val(true), maxsteps=1000)) + sol = (@test_logs min_level=Logging.Warn taylorinteg( + pendulum!, g, x0, t0, 3Tend, _order, _abstol, maxsteps=1000, dense=true)) + tv, xv, psol, tvS, xvS, gvS = sol.t, sol.x, sol.p, sol.tevents, sol.xevents, sol.gresids @test tv[1] == t0 @test xv[1,:] == x0 @test size(tvS) == (5,) @@ -91,14 +96,16 @@ import Logging: Warn @test norm(gvS, Inf) < eps() for i in 1:length(tv)-1 @test norm(psol[i,:](tv[i+1]-tv[i]) - xv[i+1,:], Inf) < 1e-13 + @test norm(sol(tv[i+1]) - xv[i+1,:], Inf) < 1e-13 end #testing 0-th order root-finding with time ranges/vectors tvr = [t0, Tend/2, Tend, 3Tend/2, 2Tend, 5Tend/2, 3Tend] @test_throws AssertionError taylorinteg(pendulum!, g, x0, view(tvr, :), _order, _abstol, maxsteps=1000, eventorder=_order+1) - xvr, tvSr, xvSr, gvSr = (@test_logs min_level=Logging.Warn taylorinteg( + solr = (@test_logs min_level=Logging.Warn taylorinteg( pendulum!, g, x0, view(tvr, :), _order, _abstol, maxsteps=1000)) + xvr, tvSr, xvSr, gvSr = solr.x, solr.tevents, solr.xevents, solr.gresids @test xvr[1,:] == x0 @test size(tvSr) == (5,) @test size(tvSr) == size(tvr[2:end-1]) @@ -109,8 +116,9 @@ import Logging: Warn @test norm(tvS-tvSr, Inf) < 5E-15 #testing 0-th order root-finding + TaylorN jet transport - tvN, xvN, tvSN, xvSN, gvSN = (@test_logs min_level=Logging.Warn taylorinteg( + solN = (@test_logs min_level=Logging.Warn taylorinteg( pendulum!, g, x0N, t0, 3Tend, _order, _abstol, maxsteps=1000)) + tvN, xvN, tvSN, xvSN, gvSN = solN.t, solN.x, solN.tevents, solN.xevents, solN.gresids @test size(tvSN) == size(tvS) @test size(xvSN) == size(xvS) @test size(gvSN) == size(gvS) @@ -119,8 +127,9 @@ import Logging: Warn @test norm( xvSN()-xvS, Inf ) < 1E-14 #testing 0-th order root-finding + TaylorN jet transport + dense output - tvN, xvN, psolN, tvSN, xvSN, gvSN = (@test_logs min_level=Logging.Warn taylorinteg( - pendulum!, g, x0N, t0, 3Tend, _order, _abstol, Val(true), maxsteps=1000)) + solN = (@test_logs min_level=Logging.Warn taylorinteg( + pendulum!, g, x0N, t0, 3Tend, _order, _abstol, maxsteps=1000, dense=true)) + tvN, xvN, psolN, tvSN, xvSN, gvSN = solN.t, solN.x, solN.p, solN.tevents, solN.xevents, solN.gresids @test size(tvSN) == size(tvS) @test size(xvSN) == size(xvS) @test size(gvSN) == size(gvS) @@ -129,24 +138,27 @@ import Logging: Warn @test norm( xvSN()-xvS, Inf ) < 1E-14 for i in 1:length(tvN)-1 @test norm(psolN[i,:](tvN[i+1]-tvN[i]) - xvN[i+1,:], Inf) < 1e-12 + @test norm(solN(tvN[i+1]) - xvN[i+1,:], Inf) < 1e-12 end #testing 0-th root-finding + Taylor1 jet transport - tv1, xv1, tvS1, xvS1, gvS1 = (@test_logs min_level=Logging.Warn taylorinteg( + sol1 = (@test_logs min_level=Logging.Warn taylorinteg( pendulum!, g, x01, t0, 3Tend, _order, _abstol, maxsteps=1000)) + tv1, xv1, tvS1, xvS1, gvS1 = sol1.t, sol1.x, sol1.tevents, sol1.xevents, sol1.gresids @test size(tvS1) == size(tvS) @test size(xvS1) == size(xvS) @test size(gvS1) == size(gvS) @test norm(gvS1[:]) < 1E-14 - @test norm( tvS1()-tvr[2:end-1], Inf ) < 1E-13 + @test norm( tvS1()-tvS, Inf ) < 1E-13 @test norm( xvS1()-xvS, Inf ) < 1E-14 #testing surface higher order crossing detections and root-finding @test_throws AssertionError taylorinteg(pendulum!, g, x0, t0, 3Tend, _order, _abstol, maxsteps=1000, eventorder=_order+1, newtoniter=2) - tv, xv, tvS, xvS, gvS = (@test_logs (Warn, err_newton_raphson()) match_mode=:any taylorinteg( + sol = (@test_logs (Warn, err_newton_raphson()) match_mode=:any taylorinteg( pendulum!, g, x0, t0, 3Tend, _order, _abstol, maxsteps=1000, eventorder=2, newtoniter=2)) + tv, xv, tvS, xvS, gvS = sol.t, sol.x, sol.tevents, sol.xevents, sol.gresids @test tv[1] < tv[end] @test tv[1] == t0 @test xv[1,:] == x0 @@ -155,8 +167,9 @@ import Logging: Warn @test norm(gvS[:]) < 1E-15 # testing backward integrations - tvb, xvb, tvSb, xvSb, gvSb = (@test_logs (Warn, err_newton_raphson()) match_mode=:any taylorinteg( + solb = (@test_logs (Warn, err_newton_raphson()) match_mode=:any taylorinteg( pendulum!, g, xv[end,:], 3Tend, t0, _order, _abstol, maxsteps=1000, eventorder=2, newtoniter=2)) + tvb, xvb, tvSb, xvSb, gvSb = solb.t, solb.x, solb.tevents, solb.xevents, solb.gresids @test tvb[1] > tvb[end] @test tvSb[1] > tvSb[end] @test norm(gvSb[:]) < 1E-14 @@ -164,8 +177,9 @@ import Logging: Warn @test norm( xvSb[2:end,:] .- xvS[end:-1:1,:], Inf ) < 1E-13 #testing higher order root-finding + TaylorN jet transport - tvN, xvN, tvSN, xvSN, gvSN = (@test_logs min_level=Logging.Warn taylorinteg( + solN = (@test_logs min_level=Logging.Warn taylorinteg( pendulum!, g, x0N, t0, 3Tend, _order, _abstol, maxsteps=1000, eventorder=2)) + tvN, xvN, tvSN, xvSN, gvSN = solN.t, solN.x, solN.tevents, solN.xevents, solN.gresids @test size(tvSN) == size(tvS) @test size(xvSN) == size(xvS) @test size(gvSN) == size(gvS) @@ -174,8 +188,9 @@ import Logging: Warn @test norm(xvSN()-xvS, Inf) < 1E-14 #testing higher root-finding + Taylor1 jet transport - tv1, xv1, tvS1, xvS1, gvS1 = (@test_logs min_level=Logging.Warn taylorinteg( + sol1 = (@test_logs min_level=Logging.Warn taylorinteg( pendulum!, g, x01, t0, 3Tend, _order, _abstol, maxsteps=1000, eventorder=2)) + tv1, xv1, tvS1, xvS1, gvS1 = sol1.t, sol1.x, sol1.tevents, sol1.xevents, sol1.gresids @test size(tvS1) == size(tvS) @test size(xvS1) == size(xvS) @test size(gvS1) == size(gvS) @@ -186,10 +201,12 @@ import Logging: Warn # Tests if trange is properly sorted Δt = (3Tend-t0)/1000 tspan = t0:Δt:(3Tend-0.125) - xv1r, tvS1r, xvS1r, gvS1r = (@test_logs min_level=Logging.Warn taylorinteg( + sol1r = (@test_logs min_level=Logging.Warn taylorinteg( pendulum!, g, x01, tspan, _order, _abstol, maxsteps=1000, eventorder=2)) - xv1rb, tvS1rb, xvS1rb, gvS1rb = (@test_logs min_level=Logging.Warn taylorinteg( + xv1r, tvS1r, xvS1r, gvS1r = sol1r.x, sol1r.tevents, sol1r.xevents, sol1r.gresids + sol1rb = (@test_logs min_level=Logging.Warn taylorinteg( pendulum!, g, xv1r[end,:], reverse(tspan), _order, _abstol, maxsteps=1000, eventorder=2)) + xv1rb, tvS1rb, xvS1rb, gvS1rb = sol1rb.x, sol1rb.tevents, sol1rb.xevents, sol1rb.gresids @test size(xv1r) == size(xv1rb) @test size(tvS1r) == size(tvS1rb) @test size(xvS1r) == size(xvS1rb) diff --git a/test/runtests.jl b/test/runtests.jl index 0c2fc5e37..07e2ba3e7 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,6 +1,7 @@ # This file is part of the TaylorIntegration.jl package; MIT licensed testfiles = ( + "solution.jl", "one_ode.jl", "many_ode.jl", "complex.jl", diff --git a/test/solution.jl b/test/solution.jl new file mode 100644 index 000000000..4c74657e8 --- /dev/null +++ b/test/solution.jl @@ -0,0 +1,32 @@ +# This file is part of the TaylorIntegration.jl package; MIT licensed + +using TaylorIntegration +using Test +using Logging +import Logging: Warn + +@testset "Testing `solution.jl`" begin + tv = [1.,2] + xv = rand(2,2) + psol = Taylor1.(rand(2,1),2) + nsteps = 2 + sol2 = TaylorIntegration.build_solution(tv, xv, psol, nsteps) + @test sol2 isa TaylorSolution{Float64, Float64, 2} + @test string(sol2) == "tspan: (1.0, 2.0), x: 2 Float64 variables" + sol1 = TaylorIntegration.build_solution(tv, Vector(xv[1,:]), Vector(psol[1,:]), nsteps) + @test sol1 isa TaylorSolution{Float64, Float64, 1} + @test string(sol1) == "tspan: (1.0, 2.0), x: 1 Float64 variable" + tv = 0.1:0.1:1.1 + xv = rand(2, length(tv)) + sol = TaylorIntegration.build_solution(tv, xv, Taylor1.(xv, 2), 9) + t1 = 0.35 + Taylor1(get_variables()[1],2) + ind, δt = TaylorIntegration.timeindex(sol, t1) + @test ind == 3 + @test δt == t1 - sol.t[ind] + tv = collect((0:0.25:2)*pi) + xv = Matrix(hcat(sin.(tv), cos.(tv))') + psolv = Matrix(hcat(sin.(tv .+ Taylor1(25)), cos.(tv .+ Taylor1(25)))') + sol = TaylorIntegration.build_solution(tv, xv, psolv, length(tv)-2) + @test sol(sol.t[1]) == sol.x[1,:] + @test norm(sol(sol.t[end]) - sol.x[end,:], Inf) < 1e-14 +end \ No newline at end of file diff --git a/test/taylorize.jl b/test/taylorize.jl index 2dfffbcf0..d0dfcf947 100644 --- a/test/taylorize.jl +++ b/test/taylorize.jl @@ -33,9 +33,14 @@ import Logging: Warn @test (@isdefined xdot1) x0 = 1.0 - tv1, xv1 = taylorinteg( xdot1, x0, t0, tf, _order, _abstol, maxsteps=1000, parse_eqs=false) - tv1p, xv1p = (@test_logs min_level=Logging.Warn taylorinteg( - xdot1, x0, t0, tf, _order, _abstol, maxsteps=1000)) + sol1 = taylorinteg( xdot1, x0, t0, tf, _order, _abstol, maxsteps=1000, + dense=false, parse_eqs=false) + tv1 = sol1.t + xv1 = sol1.x + sol1p = (@test_logs min_level=Logging.Warn taylorinteg( + xdot1, x0, t0, tf, _order, _abstol, maxsteps=1000, dense=false)) + tv1p = sol1p.t + xv1p = sol1p.x @test length(tv1) == length(tv1p) @test iszero( norm(tv1-tv1p, Inf) ) @@ -45,9 +50,14 @@ import Logging: Warn @taylorize xdot2(x, p, t) = (local b2 = 3; b2-x^2) @test (@isdefined xdot2) - tv2, xv2 = taylorinteg( xdot2, x0, t0, tf, _order, _abstol, maxsteps=1000, parse_eqs=false) - tv2p, xv2p = (@test_logs min_level=Logging.Warn taylorinteg( - xdot2, x0, t0, tf, _order, _abstol, maxsteps=1000)) + sol2 = taylorinteg( xdot2, x0, t0, tf, _order, _abstol, maxsteps=1000, + dense=false, parse_eqs=false) + tv2 = sol2.t + xv2 = sol2.x + sol2p = (@test_logs min_level=Logging.Warn taylorinteg( + xdot2, x0, t0, tf, _order, _abstol, maxsteps=1000, dense=false)) + tv2p = sol2p.t + xv2p = sol2p.x @test length(tv2) == length(tv2p) @test iszero( norm(tv2-tv2p, Inf) ) @@ -57,10 +67,14 @@ import Logging: Warn @taylorize xdot3(x, p, t) = p-x^2 @test (@isdefined xdot3) - tv3, xv3 = taylorinteg( xdot3, x0, t0, tf, _order, _abstol, b1, maxsteps=1000, - parse_eqs=false) - tv3p, xv3p = (@test_logs min_level=Logging.Warn taylorinteg( - xdot3, x0, t0, tf, _order, _abstol, b1, maxsteps=1000)) + sol3 = taylorinteg( xdot3, x0, t0, tf, _order, _abstol, b1, maxsteps=1000, + dense=false, parse_eqs=false) + tv3 = sol3.t + xv3 = sol3.x + sol3p = (@test_logs min_level=Logging.Warn taylorinteg( + xdot3, x0, t0, tf, _order, _abstol, b1, maxsteps=1000, dense=false)) + tv3p = sol3p.t + xv3p = sol3p.x @test length(tv3) == length(tv3p) @test iszero( norm(tv3-tv3p, Inf) ) @@ -73,10 +87,12 @@ import Logging: Warn @test iszero( norm(xv1p-xv2p, Inf) ) @test iszero( norm(xv1p-xv3p, Inf) ) - xv4 = taylorinteg( xdot2, x0, t0:0.5:tf, _order, _abstol, maxsteps=1000, + sol4 = taylorinteg( xdot2, x0, t0:0.5:tf, _order, _abstol, maxsteps=1000, parse_eqs=false) - xv4p = (@test_logs min_level=Logging.Warn taylorinteg( + xv4 = sol4.x + sol4p = (@test_logs min_level=Logging.Warn taylorinteg( xdot2, x0, t0:0.5:tf, _order, _abstol, maxsteps=1000)) + xv4p = sol4p.x @test iszero( norm(xv4-xv4p, Inf) ) # Compare to exact solution @@ -120,8 +136,10 @@ import Logging: Warn @test (@isdefined xdot2_err) @test !isempty(methodswith(Val{xdot2_err}, TI.jetcoeffs!)) - tv2e, xv2e = (@test_logs (Warn, unable_to_parse(xdot2_err)) taylorinteg( - xdot2_err, x0, t0, tf, _order, _abstol, maxsteps=1000)) + sol2e = (@test_logs (Warn, unable_to_parse(xdot2_err)) taylorinteg( + xdot2_err, x0, t0, tf, _order, _abstol, maxsteps=1000, dense=false)) + tv2e = sol2e.t + xv2e = sol2e.x @test length(tv2) == length(tv2e) @test iszero( norm(tv2-tv2e, Inf) ) @test iszero( norm(xv2-xv2e, Inf) ) @@ -132,12 +150,16 @@ import Logging: Warn Val(xdot2_err), tT, xT, nothing, rv) # Output includes Taylor polynomial solution - tv3t, xv3t, psol3t = (@test_logs (Warn, max_iters_reached()) taylorinteg( - xdot3, x0, t0, tf, _order, _abstol, Val(true), 0.0, maxsteps=2)) + sol3 = (@test_logs (Warn, max_iters_reached()) taylorinteg( + xdot3, x0, t0, tf, _order, _abstol, 0.0, maxsteps=2)) + tv3t = sol3.t + xv3t = sol3.x + psol3t = sol3.p @test length(psol3t) == 2 @test xv3t[1] == x0 @test psol3t[1] == Taylor1([(-1.0)^i for i=0:_order]) @test xv3t[2] == evaluate(psol3t[1], tv3t[2]-tv3t[1]) + @test xv3t[2] == sol3(tv3t[2]) end @@ -147,9 +169,12 @@ import Logging: Warn @taylorize xdot1_parsed(x, p, t) = -10 # `zero(t)` can be avoided here ! @test (@isdefined xdot1_parsed) - tv1, xv1 = taylorinteg( xdot1, 10, 1, 20.0, _order, _abstol) - tv1p, xv1p = (@test_logs min_level=Logging.Warn taylorinteg( - xdot1_parsed, 10.0, 1.0, 20.0, _order, _abstol)) + sol1 = taylorinteg( xdot1, 10, 1, 20.0, _order, _abstol, dense=false, + parse_eqs=false) + tv1, xv1 = sol1.t, sol1.x + sol1p = (@test_logs min_level=Logging.Warn taylorinteg( + xdot1_parsed, 10.0, 1.0, 20.0, _order, _abstol, dense=false)) + tv1p, xv1p = sol1p.t, sol1p.x @test length(tv1) == length(tv1p) @test iszero( norm(tv1-tv1p, Inf) ) @@ -163,9 +188,12 @@ import Logging: Warn @test (@isdefined xdot2) - tv2, xv2 = taylorinteg( xdot2, 10.0, 1.0, 20.0, _order, _abstol, parse_eqs=false) - tv2p, xv2p = (@test_logs min_level=Logging.Warn taylorinteg( - xdot2, 10.0, 1, 20.0, _order, _abstol)) + sol2 = taylorinteg( xdot2, 10.0, 1.0, 20.0, _order, _abstol, dense=false, + parse_eqs=false) + tv2, xv2 = sol2.t, sol2.x + sol2p = (@test_logs min_level=Logging.Warn taylorinteg( + xdot2, 10.0, 1, 20.0, _order, _abstol, dense=false)) + tv2p, xv2p = sol2p.t, sol2p.x @test length(tv2) == length(tv2p) @test iszero( norm(tv2-tv2p, Inf) ) @@ -176,10 +204,12 @@ import Logging: Warn @test (@isdefined xdot3) - tv3, xv3 = taylorinteg( xdot3, 10.0, 1.0, 20.0, _order, _abstol, -10, - parse_eqs=false) - tv3p, xv3p = (@test_logs min_level=Logging.Warn taylorinteg( - xdot3, 10, 1.0, 20.0, _order, _abstol, -10)) + sol3 = taylorinteg( xdot3, 10.0, 1.0, 20.0, _order, _abstol, -10, + dense=false, parse_eqs=false) + tv3, xv3 = sol3.t, sol3.x + sol3p = (@test_logs min_level=Logging.Warn taylorinteg( + xdot3, 10, 1.0, 20.0, _order, _abstol, -10, dense=false)) + tv3p, xv3p = sol3p.t, sol3p.x @test length(tv3) == length(tv3p) @test iszero( norm(tv3-tv3p, Inf) ) @@ -231,10 +261,13 @@ import Logging: Warn @test (@isdefined pendulum!) q0 = [pi-0.001, 0.0] - tv2, xv2 = taylorinteg(pendulum!, q0, t0, tf, _order, _abstol, parse_eqs=false, - maxsteps=5000) - tv2p, xv2p = (@test_logs min_level=Logging.Warn taylorinteg( - pendulum!, q0, t0, tf, _order, _abstol, maxsteps=5000)) + sol2 = taylorinteg(pendulum!, q0, t0, tf, _order, _abstol, parse_eqs=false, + maxsteps=5000, dense=false) + tv2, xv2 = sol2.t, sol2.x + + sol2p = (@test_logs min_level=Logging.Warn taylorinteg( + pendulum!, q0, t0, tf, _order, _abstol, maxsteps=5000, dense=false)) + tv2p, xv2p = sol2p.t, sol2p.x @test length(tv2) == length(tv2p) @test iszero( norm(tv2-tv2p, Inf) ) @@ -266,10 +299,12 @@ import Logging: Warn @test (@isdefined eqscmplx) cx0 = complex(1.0, 0.0) - tv1, xv1 = taylorinteg(eqscmplx, cx0, t0, tf, _order, _abstol, maxsteps=1500, - parse_eqs=false) - tv1p, xv1p = (@test_logs min_level=Logging.Warn taylorinteg( - eqscmplx, cx0, t0, tf, _order, _abstol, maxsteps=1500)) + sol1 = taylorinteg(eqscmplx, cx0, t0, tf, _order, _abstol, maxsteps=1500, + parse_eqs=false, dense=false) + tv1, xv1 = sol1.t, sol1.x + sol1p = (@test_logs min_level=Logging.Warn taylorinteg( + eqscmplx, cx0, t0, tf, _order, _abstol, maxsteps=1500, dense=false)) + tv1p, xv1p = sol1p.t, sol1p.x @test length(tv1) == length(tv1p) @test iszero( norm(tv1-tv1p, Inf) ) @@ -279,10 +314,12 @@ import Logging: Warn @taylorize eqscmplx2(x, p, t) = (local cc1 = Complex(0.0,1.0); cc1*x) @test (@isdefined eqscmplx2) - tv2, xv2 = taylorinteg(eqscmplx2, cx0, t0, tf, _order, _abstol, maxsteps=1500, - parse_eqs=false) - tv2p, xv2p = (@test_logs min_level=Logging.Warn taylorinteg( - eqscmplx2, cx0, t0, tf, _order, _abstol, maxsteps=1500)) + sol2 = taylorinteg(eqscmplx2, cx0, t0, tf, _order, _abstol, maxsteps=1500, + parse_eqs=false, dense=false) + tv2, xv2 = sol2.t, sol2.x + sol2p = (@test_logs min_level=Logging.Warn taylorinteg( + eqscmplx2, cx0, t0, tf, _order, _abstol, maxsteps=1500, dense=false)) + tv2p, xv2p = sol2p.t, sol2p.x @test length(tv2) == length(tv2p) @test iszero( norm(tv2-tv2p, Inf) ) @@ -291,10 +328,12 @@ import Logging: Warn # Passing a parameter @taylorize eqscmplx3(x, p, t) = p*x @test (@isdefined eqscmplx3) - tv3, xv3 = taylorinteg(eqscmplx3, cx0, t0, tf, _order, _abstol, cc, maxsteps=1500, - parse_eqs=false) - tv3p, xv3p = (@test_logs min_level=Logging.Warn taylorinteg( - eqscmplx3, cx0, t0, tf, _order, _abstol, cc, maxsteps=1500)) + sol3 = taylorinteg(eqscmplx3, cx0, t0, tf, _order, _abstol, cc, maxsteps=1500, + parse_eqs=false, dense=false) + tv3, xv3 = sol3.t, sol3.x + sol3p = (@test_logs min_level=Logging.Warn taylorinteg( + eqscmplx3, cx0, t0, tf, _order, _abstol, cc, maxsteps=1500, dense=false)) + tv3p, xv3p = sol3p.t, sol3p.x @test length(tv3) == length(tv3p) @test iszero( norm(tv3-tv3p, Inf) ) @@ -327,12 +366,15 @@ import Logging: Warn z0 = complex(0.0, 1.0) zz0 = [z0, z0] ts = 0.0:pi:2pi + zsol = taylorinteg(eqs3!, zz0, ts, _order, _abstol, parse_eqs=false, maxsteps=10) + tz, xz = zsol.t, zsol.x zsolp = (@test_logs min_level=Logging.Warn taylorinteg( eqs3!, zz0, ts, _order, _abstol, maxsteps=10)) - @test length(tv3) == length(tv3p) - @test iszero( norm(tv3-tv3p, Inf) ) - @test iszero( norm(xv3-xv3p, Inf) ) + tzp, xzp = zsolp.t, zsolp.x + @test length(tz) == length(tzp) + @test iszero( norm(tz-tzp, Inf) ) + @test iszero( norm(xz-xzp, Inf) ) tT = t0 + Taylor1(_order) zT = zz0 .+ zero(tT) @@ -362,16 +404,22 @@ import Logging: Warn @test (@isdefined integ_cos1) @test (@isdefined integ_cos2) - tv11, xv11 = taylorinteg(integ_cos1, 0.0, 0.0, pi, _order, _abstol, parse_eqs=false) - tv12, xv12 = (@test_logs min_level=Logging.Warn taylorinteg( - integ_cos1, 0.0, 0.0, pi, _order, _abstol)) + sol11 = taylorinteg(integ_cos1, 0.0, 0.0, pi, _order, _abstol, parse_eqs=false, + dense=false) + tv11, xv11 = sol11.t, sol11.x + sol12 = (@test_logs min_level=Logging.Warn taylorinteg( + integ_cos1, 0.0, 0.0, pi, _order, _abstol, dense=false)) + tv12, xv12 = sol12.t, sol12.x @test length(tv11) == length(tv12) @test iszero( norm(tv11-tv12, Inf) ) @test iszero( norm(xv11-xv12, Inf) ) - tv21, xv21 = taylorinteg(integ_cos2, 0.0, 0.0, pi, _order, _abstol, parse_eqs=false) - tv22, xv22 = (@test_logs min_level=Logging.Warn taylorinteg( - integ_cos2, 0.0, 0.0, pi, _order, _abstol)) + sol21 = taylorinteg(integ_cos2, 0.0, 0.0, pi, _order, _abstol, parse_eqs=false, + dense=false) + tv21, xv21 = sol21.t, sol21.x + sol22 = (@test_logs min_level=Logging.Warn taylorinteg( + integ_cos2, 0.0, 0.0, pi, _order, _abstol, dense=false)) + tv22, xv22 = sol22.t, sol22.x @test length(tv21) == length(tv22) @test iszero( norm(tv21-tv22, Inf) ) @test iszero( norm(xv21-xv22, Inf) ) @@ -393,10 +441,13 @@ import Logging: Warn @test (@isdefined integ_vec) x0 = [0.0, 1.0] - tv11, xv11 = taylorinteg(integ_vec, x0, 0.0, pi, _order, _abstol, parse_eqs=false) - tv12, xv12 = (@test_logs min_level=Logging.Warn taylorinteg( - integ_vec, x0, 0.0, pi, _order, _abstol)) - @test length(tv11) == length(tv12) + sol11 = taylorinteg(integ_vec, x0, 0.0, pi, _order, _abstol, parse_eqs=false, + dense=false) + tv11, xv11 = sol11.t, sol11.x + sol12 = (@test_logs min_level=Logging.Warn taylorinteg( + integ_vec, x0, 0.0, pi, _order, _abstol, dense=false)) + tv12, xv12 = sol12.t, sol12.x + @test length(tv11) == length(tv12) @test iszero( norm(tv11-tv12, Inf) ) @test iszero( norm(xv11-xv12, Inf) ) @@ -418,9 +469,13 @@ import Logging: Warn return du end - tvF, xvF = taylorinteg(ggg!, [1.0, 0.0], 0.0, 10.0, 20, 1.0e-20, parse_eqs=false) - tv1, xv1 = taylorinteg(ggg!, [1.0, 0.0], 0.0, 10.0, 20, 1.0e-20) - tv2, xv2 = taylorinteg(fff!, [1.0, 0.0], 0.0, 10.0, 20, 1.0e-20) + solF = taylorinteg(ggg!, [1.0, 0.0], 0.0, 10.0, 20, 1.0e-20, parse_eqs=false, + dense=false) + tvF, xvF = solF.t, solF.x + sol1 = taylorinteg(ggg!, [1.0, 0.0], 0.0, 10.0, 20, 1.0e-20) + tv1, xv1 = sol1.t, sol1.x + sol2 = taylorinteg(fff!, [1.0, 0.0], 0.0, 10.0, 20, 1.0e-20) + tv2, xv2 = sol2.t, sol2.x @test tv1 == tv2 == tvF @test xv1 == xv2 == xvF @@ -441,10 +496,12 @@ import Logging: Warn q0 = [1.0, 0.0] p = [2.0] local tf = 300.0 - tv1, xv1 = taylorinteg(harm_osc!, q0, t0, tf, _order, _abstol, - p, maxsteps=1000, parse_eqs=false) - tv1p, xv1p = (@test_logs min_level=Logging.Warn taylorinteg( - harm_osc!, q0, t0, tf, _order, _abstol, p, maxsteps=1000)) + sol1 = taylorinteg(harm_osc!, q0, t0, tf, _order, _abstol, + p, maxsteps=1000, parse_eqs=false, dense=false) + tv1, xv1 = sol1.t, sol1.x + sol1p = (@test_logs min_level=Logging.Warn taylorinteg( + harm_osc!, q0, t0, tf, _order, _abstol, p, maxsteps=1000, dense=false)) + tv1p, xv1p = sol1p.t, sol1p.x @test length(tv1) == length(tv1p) @test iszero( norm(tv1-tv1p, Inf) ) @@ -477,9 +534,11 @@ import Logging: Warn end @test !isempty(methodswith(Val{harm_osc_error!}, TI.jetcoeffs!)) - tv2e, xv2e = ( + sol2e = ( @test_logs (Warn, unable_to_parse(harm_osc_error!)) taylorinteg( - harm_osc_error!, q0, t0, tf, _order, _abstol, p, maxsteps=1000, parse_eqs=true)) + harm_osc_error!, q0, t0, tf, _order, _abstol, p, maxsteps=1000, + parse_eqs=true, dense=false)) + tv2e, xv2e = sol2e.t, sol2e.x @test length(tv1) == length(tv1p) @test iszero( norm(tv1-tv2e, Inf) ) @test iszero( norm(xv1-xv2e, Inf) ) @@ -507,10 +566,13 @@ import Logging: Warn @test (@isdefined multpendula1!) q0 = [pi-0.001, 0.0, pi-0.001, 0.0, pi-0.001, 0.0] - tv1, xv1 = taylorinteg(multpendula1!, q0, t0, tf, _order, _abstol, - [NN, nnrange], maxsteps=1000, parse_eqs=false) - tv1p, xv1p = (@test_logs min_level=Logging.Warn taylorinteg( - multpendula1!, q0, t0, tf, _order, _abstol, [NN, nnrange], maxsteps=1000)) + sol1 = taylorinteg(multpendula1!, q0, t0, tf, _order, _abstol, + [NN, nnrange], maxsteps=1000, parse_eqs=false, dense=false) + tv1, xv1 = sol1.t, sol1.x + sol1p = (@test_logs min_level=Logging.Warn taylorinteg( + multpendula1!, q0, t0, tf, _order, _abstol, [NN, nnrange], maxsteps=1000, + dense=false)) + tv1p, xv1p = sol1p.t, sol1p.x @test length(tv1) == length(tv1p) @test iszero( norm(tv1-tv1p, Inf) ) @@ -527,10 +589,12 @@ import Logging: Warn end @test (@isdefined multpendula2!) - tv2, xv2 = taylorinteg(multpendula2!, q0, t0, tf, _order, _abstol, - maxsteps=1000, parse_eqs=false) - tv2p, xv2p = (@test_logs min_level=Logging.Warn taylorinteg( - multpendula2!, q0, t0, tf, _order, _abstol, maxsteps=1000)) + sol2 = taylorinteg(multpendula2!, q0, t0, tf, _order, _abstol, + maxsteps=1000, parse_eqs=false, dense=false) + tv2, xv2 = sol2.t, sol2.x + sol2p = (@test_logs min_level=Logging.Warn taylorinteg( + multpendula2!, q0, t0, tf, _order, _abstol, maxsteps=1000, dense=false)) + tv2p, xv2p = sol2p.t, sol2p.x @test length(tv2) == length(tv2p) @test iszero( norm(tv2-tv2p, Inf) ) @@ -547,10 +611,13 @@ import Logging: Warn end @test (@isdefined multpendula3!) - tv3, xv3 = taylorinteg(multpendula3!, q0, t0, tf, _order, _abstol, - [NN, nnrange], maxsteps=1000, parse_eqs=false) - tv3p, xv3p = (@test_logs min_level=Logging.Warn taylorinteg( - multpendula3!, q0, t0, tf, _order, _abstol, [NN, nnrange], maxsteps=1000)) + sol3 = taylorinteg(multpendula3!, q0, t0, tf, _order, _abstol, + [NN, nnrange], maxsteps=1000, parse_eqs=false, dense=false) + tv3, xv3 = sol3.t, sol3.x + sol3p = (@test_logs min_level=Logging.Warn taylorinteg( + multpendula3!, q0, t0, tf, _order, _abstol, [NN, nnrange], maxsteps=1000, + dense=false)) + tv3p, xv3p = sol3p.t, sol3p.x # Comparing integrations @test length(tv1) == length(tv2) == length(tv3) @@ -619,15 +686,19 @@ import Logging: Warn @test (@isdefined kepler2!) pars = (2, -1.0) q0 = [0.2, 0.0, 0.0, 3.0] - tv1, xv1 = taylorinteg(kepler1!, q0, t0, tf, _order, _abstol, -1.0, - maxsteps=500000, parse_eqs=false) - tv1p, xv1p = (@test_logs min_level=Logging.Warn taylorinteg(kepler1!, q0, t0, tf, _order, _abstol, -1.0, - maxsteps=500000)) - - tv6p, xv6p = taylorinteg(kepler2!, q0, t0, tf, _order, _abstol, pars, - maxsteps=500000, parse_eqs=false) - tv7p, xv7p = (@test_logs min_level=Logging.Warn taylorinteg(kepler2!, q0, t0, tf, _order, _abstol, pars, - maxsteps=500000)) + sol1 = taylorinteg(kepler1!, q0, t0, tf, _order, _abstol, -1.0, + maxsteps=500000, parse_eqs=false, dense=false) + tv1, xv1 = sol1.t, sol1.x + sol1p = (@test_logs min_level=Logging.Warn taylorinteg(kepler1!, q0, t0, tf, _order, + _abstol, -1.0, maxsteps=500000, dense=false)) + tv1p, xv1p = sol1p.t, sol1p.x + + sol6p = taylorinteg(kepler2!, q0, t0, tf, _order, _abstol, pars, + maxsteps=500000, parse_eqs=false, dense=false) + tv6p, xv6p = sol6p.t, sol6p.x + sol7p = (@test_logs min_level=Logging.Warn taylorinteg(kepler2!, q0, t0, tf, _order, + _abstol, pars, maxsteps=500000, dense=false)) + tv7p, xv7p = sol7p.t, sol7p.x @test length(tv1) == length(tv1p) == length(tv6p) @test iszero( norm(tv1-tv1p, Inf) ) @@ -668,15 +739,19 @@ import Logging: Warn @test (@isdefined kepler3!) @test (@isdefined kepler4!) - tv3, xv3 = taylorinteg(kepler3!, q0, t0, tf, _order, _abstol, - maxsteps=500000, parse_eqs=false) - tv3p, xv3p = (@test_logs min_level=Logging.Warn taylorinteg(kepler3!, q0, t0, tf, _order, _abstol, - maxsteps=500000)) - - tv4, xv4 = taylorinteg(kepler4!, q0, t0, tf, _order, _abstol, - maxsteps=500000, parse_eqs=false) - tv4p, xv4p = (@test_logs min_level=Logging.Warn taylorinteg(kepler4!, q0, t0, tf, _order, _abstol, - maxsteps=500000)) + sol3 = taylorinteg(kepler3!, q0, t0, tf, _order, _abstol, + maxsteps=500000, parse_eqs=false, dense=false) + tv3, xv3 = sol3.t, sol3.x + sol3p = (@test_logs min_level=Logging.Warn taylorinteg(kepler3!, q0, t0, tf, _order, + _abstol, maxsteps=500000, dense=false)) + tv3p, xv3p = sol3p.t, sol3p.x + + sol4 = taylorinteg(kepler4!, q0, t0, tf, _order, _abstol, + maxsteps=500000, parse_eqs=false, dense=false) + tv4, xv4 = sol4.t, sol4.x + sol4p = (@test_logs min_level=Logging.Warn taylorinteg(kepler4!, q0, t0, tf, _order, + _abstol, maxsteps=500000, dense=false)) + tv4p, xv4p = sol4p.t, sol4p.x @test length(tv3) == length(tv3p) == length(tv4) @test iszero( norm(tv3-tv3p, Inf) ) @@ -758,15 +833,19 @@ import Logging: Warn @test (@isdefined kepler1!) @test (@isdefined kepler2!) q0 = [0.2, 0.0, 0.0, 3.0] - tv1, xv1 = taylorinteg(kepler1!, q0, t0, tf, _order, _abstol, - maxsteps=500000, parse_eqs=false) - tv1p, xv1p = (@test_logs min_level=Logging.Warn taylorinteg( - kepler1!, q0, t0, tf, _order, _abstol, maxsteps=500000)) - - tv2, xv2 = taylorinteg(kepler2!, q0, t0, tf, _order, _abstol, -1.0, - maxsteps=500000, parse_eqs=false) - tv2p, xv2p = (@test_logs min_level=Logging.Warn taylorinteg( - kepler2!, q0, t0, tf, _order, _abstol, -1.0, maxsteps=500000)) + sol1 = taylorinteg(kepler1!, q0, t0, tf, _order, _abstol, + maxsteps=500000, parse_eqs=false, dense=false) + tv1, xv1 = sol1.t, sol1.x + sol1p = (@test_logs min_level=Logging.Warn taylorinteg( + kepler1!, q0, t0, tf, _order, _abstol, maxsteps=500000, dense=false)) + tv1p, xv1p = sol1p.t, sol1p.x + + sol2 = taylorinteg(kepler2!, q0, t0, tf, _order, _abstol, -1.0, + maxsteps=500000, parse_eqs=false, dense=false) + tv2, xv2 = sol2.t, sol2.x + sol2p = (@test_logs min_level=Logging.Warn taylorinteg( + kepler2!, q0, t0, tf, _order, _abstol, -1.0, maxsteps=500000, dense=false)) + tv2p, xv2p = sol2p.t, sol2p.x @test length(tv1) == length(tv1p) == length(tv2) @test iszero( norm(tv1-tv1p, Inf) ) @@ -810,15 +889,19 @@ import Logging: Warn @test (@isdefined kepler3!) @test (@isdefined kepler4!) - tv3, xv3 = taylorinteg(kepler3!, q0, t0, tf, _order, _abstol, -1.0, - maxsteps=500000, parse_eqs=false) - tv3p, xv3p = (@test_logs min_level=Logging.Warn taylorinteg( - kepler3!, q0, t0, tf, _order, _abstol, -1.0, maxsteps=500000)) - - tv4, xv4 = taylorinteg(kepler4!, q0, t0, tf, _order, _abstol, - maxsteps=500000, parse_eqs=false) - tv4p, xv4p = (@test_logs min_level=Logging.Warn taylorinteg( - kepler4!, q0, t0, tf, _order, _abstol, maxsteps=500000)) + sol3 = taylorinteg(kepler3!, q0, t0, tf, _order, _abstol, -1.0, + maxsteps=500000, parse_eqs=false, dense=false) + tv3, xv3 = sol3.t, sol3.x + sol3p = (@test_logs min_level=Logging.Warn taylorinteg( + kepler3!, q0, t0, tf, _order, _abstol, -1.0, maxsteps=500000, dense=false)) + tv3p, xv3p = sol3p.t, sol3p.x + + sol4 = taylorinteg(kepler4!, q0, t0, tf, _order, _abstol, + maxsteps=500000, parse_eqs=false, dense=false) + tv4, xv4 = sol4.t, sol4.x + sol4p = (@test_logs min_level=Logging.Warn taylorinteg( + kepler4!, q0, t0, tf, _order, _abstol, maxsteps=500000, dense=false)) + tv4p, xv4p = sol4p.t, sol4p.x @test length(tv3) == length(tv3p) == length(tv4) @test iszero( norm(tv3-tv3p, Inf) ) @@ -896,31 +979,35 @@ import Logging: Warn q0 = [19.0, 20.0, 50.0] #the initial condition xi = set_variables("δ", order=1, numvars=length(q0)) - tv1, lv1, xv1 = lyap_taylorinteg(lorenz1!, q0, t0, tf, _order, _abstol, + sol1 = lyap_taylorinteg(lorenz1!, q0, t0, tf, _order, _abstol, params, maxsteps=2000, parse_eqs=false) + tv1, xv1, lv1 = sol1.t, sol1.x, sol1.λ - tv1p, lv1p, xv1p = (@test_logs min_level=Logging.Warn lyap_taylorinteg( + solp = (@test_logs min_level=Logging.Warn lyap_taylorinteg( lorenz1!, q0, t0, tf, _order, _abstol, params, maxsteps=2000)) + tv1p, xv1p, lv1p = solp.t, solp.x, solp.λ @test tv1 == tv1p - @test lv1 == lv1p @test xv1 == xv1p + @test lv1 == lv1p - tv2, lv2, xv2 = lyap_taylorinteg(lorenz1!, q0, t0, tf, _order, _abstol, + sol2 = lyap_taylorinteg(lorenz1!, q0, t0, tf, _order, _abstol, params, lorenz1_jac!, maxsteps=2000, parse_eqs=false) + tv2, xv2, lv2 = sol2.t, sol2.x, sol2.λ - tv2p, lv2p, xv2p = (@test_logs min_level=Logging.Warn lyap_taylorinteg( + sol2p = (@test_logs min_level=Logging.Warn lyap_taylorinteg( lorenz1!, q0, t0, tf, _order, _abstol, params, lorenz1_jac!, maxsteps=2000, parse_eqs=true)) + tv2p, xv2p, lv2p = sol2p.t, sol2p.x, sol2p.λ @test tv2 == tv2p - @test lv2 == lv2p @test xv2 == xv2p + @test lv2 == lv2p # Comparing both integrations (lorenz1) @test tv1 == tv2 - @test lv1 == lv2 @test xv1 == xv2 + @test lv1 == lv2 #Lorenz system ODE: @taylorize function lorenz2!(dx, x, p, t) @@ -954,53 +1041,61 @@ import Logging: Warn nothing end - tv3, lv3, xv3 = lyap_taylorinteg(lorenz2!, q0, t0, tf, _order, _abstol, + sol3 = lyap_taylorinteg(lorenz2!, q0, t0, tf, _order, _abstol, maxsteps=2000, parse_eqs=false) + tv3, xv3, lv3 = sol3.t, sol3.x, sol3.λ - tv3p, lv3p, xv3p = (@test_logs min_level=Logging.Warn lyap_taylorinteg( + sol3p = (@test_logs min_level=Logging.Warn lyap_taylorinteg( lorenz2!, q0, t0, tf, _order, _abstol, maxsteps=2000)) + tv3p, xv3p, lv3p = sol3p.t, sol3p.x, sol3p.λ @test tv3 == tv3p - @test lv3 == lv3p @test xv3 == xv3p + @test lv3 == lv3p - tv4, lv4, xv4 = lyap_taylorinteg(lorenz2!, q0, t0, tf, _order, _abstol, + sol4 = lyap_taylorinteg(lorenz2!, q0, t0, tf, _order, _abstol, lorenz2_jac!, maxsteps=2000, parse_eqs=false) + tv4, xv4, lv4 = sol4.t, sol4.x, sol4.λ - tv4p, lv4p, xv4p = (@test_logs min_level=Logging.Warn lyap_taylorinteg( + sol4p = (@test_logs min_level=Logging.Warn lyap_taylorinteg( lorenz2!, q0, t0, tf, _order, _abstol, lorenz2_jac!, maxsteps=2000, parse_eqs=true)) + tv4p, xv4p, lv4p = sol4p.t, sol4p.x, sol4p.λ @test tv4 == tv4p - @test lv4 == lv4p @test xv4 == xv4p + @test lv4 == lv4p # Comparing both integrations (lorenz2) @test tv3 == tv4 - @test lv3 == lv4 @test xv3 == xv4 + @test lv3 == lv4 # Comparing both integrations (lorenz1! vs lorenz2!) @test tv1 == tv3 - @test lv1 == lv3 @test xv1 == xv3 + @test lv1 == lv3 # Using ranges - lv5, xv5 = lyap_taylorinteg(lorenz2!, q0, t0:0.125:tf, _order, _abstol, + sol5 = lyap_taylorinteg(lorenz2!, q0, t0:0.125:tf, _order, _abstol, maxsteps=2000, parse_eqs=false) + lv5, xv5 = sol5.x, sol5.λ - lv5p, xv5p = (@test_logs min_level=Logging.Warn lyap_taylorinteg( + sol5p = (@test_logs min_level=Logging.Warn lyap_taylorinteg( lorenz2!, q0, t0:0.125:tf, _order, _abstol, maxsteps=2000, parse_eqs=true)) + lv5p, xv5p = sol5p.x, sol5p.λ @test lv5 == lv5p @test xv5 == xv5p - lv6, xv6 = lyap_taylorinteg(lorenz2!, q0, t0:0.125:tf, _order, _abstol, + sol6 = lyap_taylorinteg(lorenz2!, q0, t0:0.125:tf, _order, _abstol, lorenz2_jac!, maxsteps=2000, parse_eqs=false) + lv6, xv6 = sol6.x, sol6.λ - lv6p, xv6p = (@test_logs min_level=Logging.Warn lyap_taylorinteg( + sol6p = (@test_logs min_level=Logging.Warn lyap_taylorinteg( lorenz2!, q0, t0:0.125:tf, _order, _abstol, lorenz2_jac!, maxsteps=2000, parse_eqs=true)) + lv6p, xv6p = sol6p.x, sol6p.λ @test lv6 == lv6p @test xv6 == xv6p @@ -1133,41 +1228,44 @@ import Logging: Warn #the time range tr = t0:integstep:T; #note that as called below, taylorinteg uses the parsed jetcoeffs! method by default - xvp = (@test_logs min_level=Logging.Warn taylorinteg( + solp = (@test_logs min_level=Logging.Warn taylorinteg( pendulum!, q0, tr, _order, _abstol, maxsteps=100)) + @test tr == solp.t + xvp = solp.x # "warmup" for jet transport integration - xvTN = (@test_logs (Warn, max_iters_reached()) taylorinteg( + solTN = (@test_logs (Warn, max_iters_reached()) @inferred taylorinteg( pendulum!, q0TN, tr, _order, _abstol, maxsteps=1, parse_eqs=false)) - xvTN = (@test_logs (Warn, max_iters_reached()) taylorinteg( - pendulum!, q0TN, tr, _order, _abstol, maxsteps=1)) - @test size(xvTN) == (5,2) + @test size(solTN.x) == (5,2) #jet transport integration with parsed jetcoeffs! - xvTNp = (@test_logs min_level=Logging.Warn taylorinteg( + solTNp = (@test_logs min_level=Logging.Warn taylorinteg( pendulum!, q0TN, tr, _order, _abstol, maxsteps=100)) #jet transport integration with non-parsed jetcoeffs! - xvTN = (@test_logs min_level=Logging.Warn taylorinteg( + solTN = (@test_logs min_level=Logging.Warn taylorinteg( pendulum!, q0TN, tr, _order, _abstol, maxsteps=100, parse_eqs=false)) - @test xvTN == xvTNp - @test norm(xvTNp[:,:]() - xvp, Inf) < 1e-15 + @test solTN.x == solTNp.x + @test norm(solTNp.x[:,:]() - xvp, Inf) < 1e-15 dq = 0.0001rand(2) q1 = q0 + dq - yv = (@test_logs min_level=Logging.Warn taylorinteg( + y = (@test_logs min_level=Logging.Warn taylorinteg( pendulum!, q1, tr, _order, _abstol, maxsteps=100)) - yv_jt = xvTNp[:,:](dq) - @test norm(yv-yv_jt, Inf) < 1e-11 + y_jt = solTNp.x[:,:](dq) + @test norm(y.x-y_jt, Inf) < 1e-11 dq = 0.001 t = Taylor1([0.0, 1.0], 10) x0T1 = q0+[0t,t] q1 = q0+[0.0,dq] - tv, xv = (@test_logs min_level=Logging.Warn taylorinteg( - pendulum!, q1, t0, 2T, _order, _abstol)) - tvT1, xvT1 = (@test_logs min_level=Logging.Warn taylorinteg( - pendulum!, x0T1, t0, 2T, _order, _abstol, parse_eqs=false)) - tvT1p, xvT1p = (@test_logs min_level=Logging.Warn taylorinteg( - pendulum!, x0T1, t0, 2T, _order, _abstol)) + sol = (@test_logs min_level=Logging.Warn taylorinteg( + pendulum!, q1, t0, 2T, _order, _abstol, dense=false)) + tv, xv = sol.t, sol.x + solT1 = (@test_logs min_level=Logging.Warn taylorinteg( + pendulum!, x0T1, t0, 2T, _order, _abstol, parse_eqs=false, dense=false)) + tvT1, xvT1 = solT1.t, solT1.x + solT1p = (@test_logs min_level=Logging.Warn taylorinteg( + pendulum!, x0T1, t0, 2T, _order, _abstol, dense=false)) + tvT1p, xvT1p = solT1p.t, solT1p.x @test tvT1 == tvT1p @test xvT1 == xvT1p xv_jt = xvT1p[:,:](dq) @@ -1175,12 +1273,10 @@ import Logging: Warn @taylorize function kepler1!(dq, q, p, t) local μ = -1.0 - x1s = q[3] - y1s = q[4] - r_p2 = ((x1s^2)+(y1s^2)) + r_p2 = ((q[1]^2)+(q[2]^2)) r_p3d2 = r_p2^1.5 - dq[1] = x1s - dq[2] = y1s + dq[1] = q[3] + dq[2] = q[4] newtonianCoeff = μ / r_p3d2 dq[3] = q[1] * newtonianCoeff dq[4] = q[2] * newtonianCoeff @@ -1192,23 +1288,68 @@ import Logging: Warn q0 = [0.2, 0.0, 0.0, 3.0] q0TN = q0 + p # JT initial condition - tv, xv, psol = (@test_logs (Warn, max_iters_reached()) taylorinteg( - kepler1!, q0TN, t0, tf, _order, _abstol, Val(true), maxsteps=2, parse_eqs=false)) - tvp, xvp, psolp = (@test_logs (Warn, max_iters_reached()) taylorinteg( - kepler1!, q0TN, t0, tf, _order, _abstol, Val(true), maxsteps=2)) + sol = (@test_logs (Warn, max_iters_reached()) taylorinteg( + kepler1!, q0TN, t0, tf, _order, _abstol, maxsteps=2, parse_eqs=false, dense=false)) + tv, xv, psol = sol.t, sol.x, sol.p + solp = (@test_logs (Warn, max_iters_reached()) taylorinteg( + kepler1!, q0TN, t0, tf, _order, _abstol, maxsteps=2, dense=false)) + tvp, xvp, psolp = sol.t, sol.x, sol.p @test tv == tvp @test xv == xvp @test psol == psolp - tv, xv, psol = (@test_logs min_level=Logging.Warn taylorinteg( - kepler1!, q0TN, t0, tf, _order, _abstol, Val(true), maxsteps=3000, parse_eqs=false)) - tvp, xvp, psolp = (@test_logs min_level=Logging.Warn taylorinteg( - kepler1!, q0TN, t0, tf, _order, _abstol, Val(true), maxsteps=3000)) + sol = (@test_logs min_level=Logging.Warn taylorinteg( + kepler1!, q0TN, t0, tf, _order, _abstol, maxsteps=3000, parse_eqs=false)) + solp = (@test_logs min_level=Logging.Warn taylorinteg( + kepler1!, q0TN, t0, tf, _order, _abstol, maxsteps=3000)) + + @test length(sol.t) == length(solp.t) + @test sol.t == solp.t + @test solp.t[end] == tf + @test iszero( norm(sol.x-solp.x, Inf) ) + @test iszero( norm(sol.p-solp.p, Inf) ) + + # Keplerian energy + function visviva(x) + v2 = x[3]^2 + x[4]^2 + r = sqrt(x[1]^2 + x[2]^2) + return 0.5v2 - 1/r + end - @test length(tv) == length(tvp) - @test tv == tvp - @test iszero( norm(xv-xvp, Inf) ) - @test iszero( norm(psol-psolp, Inf) ) + # initial energy + E0 = visviva(solp(t0)) + # final energy + Ef = visviva(solp(tf)) + @test norm( E0() - Ef(), Inf ) < 1e-12 + + @testset "Test _taylorinteg function barrier" begin + _t0 = t0 + _tmax = tf + _q0 = q0TN + _params = nothing + _maxsteps = 3000 + _T = eltype(_t0) + _U = eltype(_q0) + # Initialize the vector of Taylor1 expansions + _dof = length(_q0) + _t = _t0 + Taylor1( _T, _order ) + _x = Array{Taylor1{_U}}(undef, _dof) + _dx = Array{Taylor1{_U}}(undef, _dof) + @inbounds for i in eachindex(_q0) + _x[i] = Taylor1( _q0[i], _order ) + _dx[i] = Taylor1( zero(_q0[i]), _order ) + end + # Determine if specialized jetcoeffs! method exists + __parse_eqs, __rv = TaylorIntegration._determine_parsing!(true, kepler1!, _t, _x, _dx, _params); + # Re-initialize the Taylor1 expansions + _t = _t0 + Taylor1( _T, _order ) + _x .= Taylor1.( _q0, _order ) + _dx .= Taylor1.( zero.(_q0), _order) + solTN = @inferred TaylorIntegration._taylorinteg!(Val(true), kepler1!, _t, _x, _dx, _q0, _t0, _tmax, _abstol, __rv, _params; parse_eqs=__parse_eqs, maxsteps=_maxsteps) + solTN2 = @inferred TaylorIntegration._taylorinteg!(Val(false), kepler1!, _t, _x, _dx, _q0, _t0, _tmax, _abstol, __rv, _params; parse_eqs=__parse_eqs, maxsteps=_maxsteps) + @test solTN isa TaylorSolution{typeof(_t0), eltype(_q0), ndims(solTN.x), typeof(solTN.t), typeof(solTN.x), typeof(solTN.p), Nothing, Nothing, Nothing} + @test solTN2 isa TaylorSolution{typeof(_t0), eltype(_q0), ndims(solTN.x), typeof(solTN.t), typeof(solTN.x), Nothing, Nothing, Nothing, Nothing} + end end @@ -1240,8 +1381,9 @@ import Logging: Warn pendulum!, g, x0, t0, Tend, _order, _abstol, Val(false), maxsteps=1, eventorder=_order+1) #testing 0-th order root-finding - tv, xv, tvS, xvS, gvS = (@test_logs min_level=Logging.Warn taylorinteg( + sol = (@test_logs min_level=Logging.Warn taylorinteg( pendulum!, g, x0, t0, 3Tend, _order, _abstol, Val(false), maxsteps=1000)) + tv, xv, tvS, xvS, gvS = sol.t, sol.x, sol.tevents, sol.xevents, sol.gresids @test tv[1] == t0 @test xv[1,:] == x0 @test size(tvS) == (2,) @@ -1254,8 +1396,9 @@ import Logging: Warn pendulum!, g, x0, view(tvr, :), _order, _abstol, maxsteps=1) @test_throws AssertionError taylorinteg( pendulum!, g, x0, view(tvr, :), _order, _abstol, maxsteps=1, eventorder=_order+1) - xvr, tvSr, xvSr, gvSr = (@test_logs min_level=Logging.Warn taylorinteg( + solr = (@test_logs min_level=Logging.Warn taylorinteg( pendulum!, g, x0, view(tvr, :), _order, _abstol, maxsteps=1000)) + xvr, tvSr, xvSr, gvSr = solr.x, solr.tevents, solr.xevents, solr.gresids @test xvr[1,:] == x0 @test size(tvSr) == (2,) @test norm(tvSr-tvr[3:2:end-1], Inf) < 1e-13 @@ -1503,10 +1646,12 @@ import Logging: Warn @test x == x_ @test dx == dx_ - tv, xv = (@test_logs min_level=Warn taylorinteg( + sol = (@test_logs min_level=Warn taylorinteg( harmosc1dchain!, x0, t0, 10000.0, _order, _abstol, μ)) - tv_, xv_ = (@test_logs min_level=Warn taylorinteg( + tv, xv = sol.t, sol.x + sol_ = (@test_logs min_level=Warn taylorinteg( harmosc1dchain_threads!, x0, t0, 10000.0, _order, _abstol, μ)) + tv_, xv_ = sol_.t, sol_.x @test tv == tv_ @test xv == xv_