From 4a01df0e7938e9228d9342c15c74994684267888 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jorge=20A=2E=20P=C3=A9rez-Hern=C3=A1ndez?= Date: Fri, 20 Dec 2024 12:24:52 -0600 Subject: [PATCH] Add Taylor expansions to TaylorIntegration cache types --- src/integrator/cache.jl | 236 ++++++++++++++++++++++------------ src/integrator/taylorinteg.jl | 59 +++------ src/lyapunovspectrum.jl | 95 ++++++-------- src/rootfinding.jl | 39 ++---- test/taylorize.jl | 25 +--- 5 files changed, 232 insertions(+), 222 deletions(-) diff --git a/src/integrator/cache.jl b/src/integrator/cache.jl index 0f2151a8..ecabc67a 100644 --- a/src/integrator/cache.jl +++ b/src/integrator/cache.jl @@ -24,10 +24,108 @@ struct ScalarCache{TV, XV, PSOL, T, X} <: AbstractTaylorIntegrationCache x::X end -function init_cache(cachetype::Type{ScalarCache}, dense::Val{D}, t0::T, x0::U, maxsteps::Int, order::Int) where {D, U, T} - # Initialize the Taylor1 expansions +# VectorCache + +struct VectorCache{TV, XV, PSOL, XAUX, T, X, DX} <: AbstractVectorCache + tv::TV + xv::XV + psol::PSOL + xaux::XAUX + t::T + x::X + dx::DX +end + +# VectorTRangeCache + +struct VectorTRangeCache{TV, XV, PSOL, XAUX, X0, X1, T, X, DX} <: AbstractVectorCache + tv::TV + xv::XV + psol::PSOL + xaux::XAUX + x0::X0 + x1::X1 + t::T + x::X + dx::DX +end + +# LyapunovSpectrumCache + +struct LyapunovSpectrumCache{TV, XV, PSOL, XAUX, X0, Λ, ΛTSUM, ΔX, DΔX, JAC, VARSAUX, QQH, RRH, AJ, QI, VJ, T, X, DX, JT, DVARS} <: AbstractVectorCache + tv::TV + xv::XV + psol::PSOL + xaux::XAUX + x0::X0 + λ::Λ + λtsum::ΛTSUM + δx::ΔX + dδx::DΔX + jac::JAC + varsaux::VARSAUX + QH::QQH + RH::RRH + aⱼ::AJ + qᵢ::QI + vⱼ::VJ + t::T + x::X + dx::DX + jt::JT + dvars::DVARS +end + +# LyapunovSpectrumTRangeCache + +struct LyapunovSpectrumTRangeCache{TV, XV, PSOL, XAUX, X0, Q1, Λ, ΛTSUM, ΔX, DΔX, JAC, VARSAUX, QQH, RRH, AJ, QI, VJ, T, X, DX, JT, DVARS} <: AbstractVectorCache + tv::TV + xv::XV + psol::PSOL + xaux::XAUX + x0::X0 + q1::Q1 + λ::Λ + λtsum::ΛTSUM + δx::ΔX + dδx::DΔX + jac::JAC + varsaux::VARSAUX + QH::QQH + RH::RRH + aⱼ::AJ + qᵢ::QI + vⱼ::VJ + t::T + x::X + dx::DX + jt::JT + dvars::DVARS +end + +# init_expansions + +function init_expansions(t0::T, x0::U, order::Int) where {T, U} t = t0 + Taylor1( T, order ) x = Taylor1( x0, order ) + return t, x +end + +function init_expansions(t0::T, q0::Vector{U}, order::Int) where {T, U} + dof = length(q0) + 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 ) + return t, x, dx +end + +# init_cache + +function init_cache(cachetype::Type{ScalarCache}, dense::Val{D}, t0::T, x0::U, maxsteps::Int, order::Int) where {D, U, T} + # Initialize the Taylor1 expansions + t, x = init_expansions(t0, x0, order) # Initialize cache return cachetype( Array{T}(undef, maxsteps + 1), @@ -40,8 +138,7 @@ end function init_cache(cachetype::Type{ScalarCache}, ::Val{false}, trange::AbstractVector{T}, x0::U, maxsteps::Int, order::Int) where {U, T} # Initialize the Taylor1 expansions t0 = trange[1] - t = t0 + Taylor1( T, order ) - x = Taylor1( x0, order ) + t, x = init_expansions(t0, x0, order) # Initialize cache nn = length(trange) cache = cachetype( @@ -54,45 +151,38 @@ function init_cache(cachetype::Type{ScalarCache}, ::Val{false}, trange::Abstract return cache end -# VectorCache - -struct VectorCache{TV, XV, PSOL, XAUX} <: AbstractVectorCache - tv::TV - xv::XV - psol::PSOL - xaux::XAUX -end - -function init_cache(cachetype::Type{VectorCache}, dense::Val{D}, t0::T, x::Vector{Taylor1{U}}, maxsteps::Int) where {D, U, T} - dof = length(x) +function init_cache(cachetype::Type{VectorCache}, dense::Val{D}, t0::T, q0::Vector{U}, maxsteps::Int, order::Int) where {D, U, T} + # Initialize the vector of Taylor1 expansions + t, x, dx = init_expansions(t0, q0, order) + # Initialize cache + dof = length(q0) return cachetype( Array{T}(undef, maxsteps + 1), Array{U}(undef, dof, maxsteps + 1), init_psol(dense, maxsteps, dof, x), - Array{Taylor1{U}}(undef, dof)) -end - -# VectorTRangeCache - -struct VectorTRangeCache{TV, XV, PSOL, XAUX, X0, X1} <: AbstractVectorCache - tv::TV - xv::XV - psol::PSOL - xaux::XAUX - x0::X0 - x1::X1 + Array{Taylor1{U}}(undef, dof), + t, + x, + dx) end -function init_cache(cachetype::Type{VectorTRangeCache}, ::Val{false}, trange::AbstractVector{T}, x::Vector{Taylor1{U}}, maxsteps::Int) where {U, T} +function init_cache(cachetype::Type{VectorTRangeCache}, ::Val{false}, trange::AbstractVector{T}, q0::Vector{U}, maxsteps::Int, order::Int) where {U, T} + # Initialize the vector of Taylor1 expansions + t0 = trange[1] + t, x, dx = init_expansions(t0, q0, order) + # Initialize cache nn = length(trange) - dof = length(x) + dof = length(q0) cache = cachetype( trange, Array{U}(undef, dof, nn), init_psol(Val(false), maxsteps, dof, x), Array{Taylor1{U}}(undef, dof), - similar(constant_term.(x)), - similar(constant_term.(x))) + similar(q0), + similar(q0), + t, + x, + dx) fill!(cache.x0, T(NaN)) for ind in 1:nn @inbounds cache.xv[:,ind] .= cache.x0 @@ -100,36 +190,21 @@ function init_cache(cachetype::Type{VectorTRangeCache}, ::Val{false}, trange::Ab return cache end -# LyapunovSpectrumCache - -struct LyapunovSpectrumCache{TV, XV, PSOL, XAUX, X0, Λ, ΛTSUM, ΔX, DΔX, JAC, VARSAUX, QQH, RRH, AJ, QI, VJ} <: AbstractVectorCache - tv::TV - xv::XV - psol::PSOL - xaux::XAUX - x0::X0 - λ::Λ - λtsum::ΛTSUM - δx::ΔX - dδx::DΔX - jac::JAC - varsaux::VARSAUX - QH::QQH - RH::RRH - aⱼ::AJ - qᵢ::QI - vⱼ::VJ -end - -function init_cache(cachetype::Type{LyapunovSpectrumCache}, dense, t0::T, x::Vector{Taylor1{U}}, maxsteps::Int) where {U, T} - nx0 = length(x) # equals dof + dof^2 - dof = Int(sqrt(nx0 + 1/4) - 1/2) +function init_cache(cachetype::Type{LyapunovSpectrumCache}, dense, t0::T, q0::Vector{U}, maxsteps::Int, order::Int) where {U, T} + # Initialize the vector of Taylor1 expansions + dof = length(q0) + jt = Matrix{U}(I, dof, dof) + x0 = vcat(q0, reshape(jt, dof*dof)) + t, x, dx = init_expansions(t0, x0, order) + # Initialize cache + nx0 = length(x0) + dvars = Array{TaylorN{Taylor1{U}}}(undef, dof) cache = cachetype( Array{T}(undef, maxsteps+1), Array{U}(undef, dof, maxsteps+1), nothing, Array{Taylor1{U}}(undef, nx0), - getcoeff.(x, 0), + x0, Array{U}(undef, dof, maxsteps+1), similar(constant_term.(x[1:dof])), Array{TaylorN{Taylor1{U}}}(undef, dof), @@ -140,38 +215,28 @@ function init_cache(cachetype::Type{LyapunovSpectrumCache}, dense, t0::T, x::Vec Array{U}(undef, dof, dof), Array{U}(undef, dof), Array{U}(undef, dof), - Array{U}(undef, dof) + Array{U}(undef, dof), + t, + x, + dx, + jt, + dvars ) fill!(cache.jac, zero(x[1])) return cache end -# LyapunovSpectrumTRangeCache - -struct LyapunovSpectrumTRangeCache{TV, XV, PSOL, XAUX, X0, Q1, Λ, ΛTSUM, ΔX, DΔX, JAC, VARSAUX, QQH, RRH, AJ, QI, VJ} <: AbstractVectorCache - tv::TV - xv::XV - psol::PSOL - xaux::XAUX - x0::X0 - q1::Q1 - λ::Λ - λtsum::ΛTSUM - δx::ΔX - dδx::DΔX - jac::JAC - varsaux::VARSAUX - QH::QQH - RH::RRH - aⱼ::AJ - qᵢ::QI - vⱼ::VJ -end - -function init_cache(cachetype::Type{LyapunovSpectrumTRangeCache}, dense, trange::AbstractVector{T}, x::Vector{Taylor1{U}}, maxsteps::Int) where {U, T} - nx0 = length(x) # equals dof + dof^2 - dof = Int(sqrt(nx0 + 1/4) - 1/2) +function init_cache(cachetype::Type{LyapunovSpectrumTRangeCache}, dense, trange::AbstractVector{T}, q0::Vector{U}, maxsteps::Int, order::Int) where {U, T} + # Initialize the vector of Taylor1 expansions + t0 = trange[1] + dof = length(q0) + jt = Matrix{U}(I, dof, dof) + x0 = vcat(q0, reshape(jt, dof*dof)) + t, x, dx = init_expansions(t0, x0, order) + # Initialize cache nn = length(trange) + nx0 = length(x0) + dvars = Array{TaylorN{Taylor1{U}}}(undef, dof) cache = cachetype( trange, Array{U}(undef, dof, nn), @@ -189,7 +254,12 @@ function init_cache(cachetype::Type{LyapunovSpectrumTRangeCache}, dense, trange: Array{U}(undef, dof, dof), Array{U}(undef, dof), Array{U}(undef, dof), - Array{U}(undef, dof) + Array{U}(undef, dof), + t, + x, + dx, + jt, + dvars ) fill!(cache.xv, U(NaN)) fill!(cache.λ, U(NaN)) diff --git a/src/integrator/taylorinteg.jl b/src/integrator/taylorinteg.jl index 9d679b4d..87e96600 100644 --- a/src/integrator/taylorinteg.jl +++ b/src/integrator/taylorinteg.jl @@ -57,14 +57,14 @@ function taylorinteg(f, x0::U, t0::T, tmax::T, order::Int, abstol::T, params = n # Determine if specialized jetcoeffs! method exists parse_eqs, rv = _determine_parsing!(parse_eqs, f, cache.t, cache.x, params) - return _taylorinteg!(Val(dense), f, x0, t0, tmax, abstol, rv, cache, params; parse_eqs, maxsteps) + return taylorinteg!(Val(dense), f, x0, t0, tmax, abstol, rv, cache, params; parse_eqs, maxsteps) end -function _taylorinteg!(dense::Val{D}, f, +function taylorinteg!(dense::Val{D}, f, x0::U, t0::T, tmax::T, abstol::T, rv::RetAlloc{Taylor1{U}}, cache::ScalarCache, params; parse_eqs::Bool=true, maxsteps::Int=500) where {T<:Real,U<:Number,D} - @unpack t, x, tv, xv, psol = cache + @unpack tv, xv, psol, t, x = cache # Initial conditions nsteps = 1 @@ -98,32 +98,24 @@ function _taylorinteg!(dense::Val{D}, f, end -function taylorinteg(f!, q0::Array{U,1}, t0::T, tmax::T, order::Int, abstol::T, params = nothing; +function taylorinteg(f!, q0::Vector{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 vector of Taylor1 expansions - dof = length(q0) - 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 ) - # Allocation - cache = init_cache(VectorCache, Val(dense), t0, x, maxsteps) + cache = init_cache(VectorCache, Val(dense), t0, q0, maxsteps, order) # Determine if specialized jetcoeffs! method exists - parse_eqs, rv = _determine_parsing!(parse_eqs, f!, t, x, dx, params) + parse_eqs, rv = _determine_parsing!(parse_eqs, f!, cache.t, cache.x, cache.dx, params) - return _taylorinteg!(Val(dense), f!, t, x, dx, q0, t0, tmax, abstol, rv, + return taylorinteg!(Val(dense), f!, q0, t0, tmax, abstol, rv, cache, params; parse_eqs, maxsteps) end -function _taylorinteg!(dense::Val{D}, f!, t::Taylor1{T}, x::Array{Taylor1{U},1}, dx::Array{Taylor1{U},1}, +function taylorinteg!(dense::Val{D}, f!, q0::Array{U,1}, t0::T, tmax::T, abstol::T, rv::RetAlloc{Taylor1{U}}, cache::VectorCache, params; parse_eqs::Bool=true, maxsteps::Int=500) where {T<:Real,U<:Number,D} - @unpack tv, xv, psol, xaux = cache + @unpack tv, xv, psol, xaux, t, x, dx = cache # Initial conditions @inbounds t[0] = t0 @@ -257,7 +249,7 @@ function taylorinteg(f, x0::U, trange::AbstractVector{T}, # Check if trange is increasingly or decreasingly sorted @assert (issorted(trange) || - issorted(reverse(trange))) "`trange` or `reverse(trange)` must be sorted" + issorted(trange, rev=true)) "`trange` or `reverse(trange)` must be sorted" # Allocation cache = init_cache(ScalarCache, Val(false), trange, x0, maxsteps, order) @@ -265,13 +257,13 @@ function taylorinteg(f, x0::U, trange::AbstractVector{T}, # Determine if specialized jetcoeffs! method exists parse_eqs, rv = _determine_parsing!(parse_eqs, f, cache.t, cache.x, params) - return _taylorinteg!(f, x0, trange, abstol, rv, cache, params; parse_eqs, maxsteps) + return taylorinteg!(f, x0, trange, abstol, rv, cache, params; parse_eqs, maxsteps) end -function _taylorinteg!(f, x0::U, trange::AbstractVector{T}, +function taylorinteg!(f, x0::U, trange::AbstractVector{T}, abstol::T, rv::RetAlloc{Taylor1{U}}, cache::ScalarCache, params; parse_eqs::Bool=true, maxsteps::Int=500) where {T<:Real, U<:Number} - @unpack t, x, xv = cache + @unpack xv, t, x = cache # Initial conditions @inbounds t0, t1, tmax = trange[1], trange[2], trange[end] @@ -313,38 +305,29 @@ function _taylorinteg!(f, x0::U, trange::AbstractVector{T}, return build_solution(trange, xv) end -function taylorinteg(f!, q0::Array{U,1}, trange::AbstractVector{T}, +function taylorinteg(f!, q0::Vector{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 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 ) + issorted(trange, rev=true)) "`trange` or `reverse(trange)` must be sorted" # Allocation - cache = init_cache(VectorTRangeCache, Val(false), trange, x, maxsteps) + cache = init_cache(VectorTRangeCache, Val(false), trange, q0, maxsteps, order) # Determine if specialized jetcoeffs! method exists - parse_eqs, rv = _determine_parsing!(parse_eqs, f!, t, x, dx, params) + parse_eqs, rv = _determine_parsing!(parse_eqs, f!, cache.t, cache.x, cache.dx, params) - return _taylorinteg!(f!, t, x, dx, q0, trange, abstol, rv, + return taylorinteg!(f!, q0, trange, abstol, rv, cache, 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}}, cache::VectorTRangeCache, params; +function taylorinteg!(f!, + q0::Vector{U}, trange::AbstractVector{T}, abstol::T, rv::RetAlloc{Taylor1{U}}, cache::VectorTRangeCache, params; parse_eqs::Bool=true, maxsteps::Int=500) where {T<:Real, U<:Number} - @unpack xv, xaux, x0, x1 = cache + @unpack xv, xaux, x0, x1, t, x, dx = cache # Initial conditions @inbounds t[0] = trange[1] diff --git a/src/lyapunovspectrum.jl b/src/lyapunovspectrum.jl index d0cc0edb..f656b3d9 100644 --- a/src/lyapunovspectrum.jl +++ b/src/lyapunovspectrum.jl @@ -36,11 +36,11 @@ build_lyap_solution(t::AbstractVector{T}, # stabilitymatrix! """ - stabilitymatrix!(eqsdiff!, t, x, δx, dδx, jac, _δv, params[, jacobianfunc!=nothing]) + stabilitymatrix!(eqsdiff!, t, x, δx, dδx, jac, dvars, params[, jacobianfunc!=nothing]) Updates the matrix `jac::Matrix{Taylor1{U}}` (linearized equations of motion) computed from the equations of motion (`eqsdiff!`), at time `t` at `x`; `x` is -of type `Vector{Taylor1{U}}`, where `U<:Number`. `δx`, `dδx` and `_δv` are +of type `Vector{Taylor1{U}}`, where `U<:Number`. `δx`, `dδx` and `dvars` are auxiliary arrays of type `Vector{TaylorN{Taylor1{U}}}` to avoid allocations. Optionally, the user may provide a Jacobian function `jacobianfunc!` to compute `jac`. Otherwise, `jac` is computed via automatic differentiation using @@ -49,11 +49,11 @@ 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, + jac::Matrix{Taylor1{U}}, dvars::Vector{TaylorN{Taylor1{U}}}, params, ::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] + @inbounds δx[ind] = x[ind] + dvars[ind] end # Equations of motion eqsdiff!(dδx, δx, params, t) @@ -62,13 +62,13 @@ function stabilitymatrix!(eqsdiff!, t::Taylor1{T}, x::Vector{Taylor1{U}}, 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, + jac::Matrix{Taylor1{U}}, dvars::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) +stabilitymatrix!(eqsdiff!, t, x, δx, dδx, jac, dvars, params) = + stabilitymatrix!(eqsdiff!, t, x, δx, dδx, jac, dvars, params, nothing) # Modified from `cgs` and `mgs`, obtained from: # http://nbviewer.jupyter.org/url/math.mit.edu/~stevenj/18.335/Gram-Schmidt.ipynb @@ -194,12 +194,12 @@ function lyap_jetcoeffs!(t::Taylor1{T}, x::AbstractVector{Taylor1{S}}, end """ - lyap_taylorstep!(::Val{V}, f!, t, x, dx, xaux, δx, dδx, jac, abstol, _δv, varsaux, params, rv[, jacobianfunc!]) + lyap_taylorstep!(::Val{V}, f!, t, x, dx, xaux, δx, dδx, jac, abstol, dvars, 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 linearization of the equations of motion, i.e, the Jacobian. `xaux`, `δx`, `dδx`, -`varsaux` and `_δv` are auxiliary vectors, and `params` define the parameters +`varsaux` and `dvars` are auxiliary vectors, and `params` define the parameters of the ODEs. Optionally, the user may provide a Jacobian function `jacobianfunc!` to compute `jac`. Otherwise, `jac` is computed via automatic differentiation using `TaylorSeries.jl`. @@ -208,7 +208,7 @@ via automatic differentiation using `TaylorSeries.jl`. 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}}}, + jac::Array{Taylor1{U},2}, abstol::T, dvars::Vector{TaylorN{Taylor1{U}}}, varsaux::Array{Taylor1{U},3}, params, rv::RetAlloc{Taylor1{U}}, jacobianfunc! =nothing) where {T<:Real, U<:Number, V} @@ -220,7 +220,7 @@ function lyap_taylorstep!(::Val{V}, f!, t::Taylor1{T}, x::Vector{Taylor1{U}}, dx __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!) + stabilitymatrix!(f!, t, x, δx, dδx, jac, dvars, 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) @@ -252,42 +252,33 @@ function lyap_taylorinteg(f!, q0::Array{U,1}, t0::T, tmax::T, order::Int, abstol::T, params = nothing, jacobianfunc! =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 ) - jt = Matrix{U}(I, dof, dof) - x0 = vcat(q0, reshape(jt, dof*dof)) - nx0 = length(x0) - x = Array{Taylor1{U}}(undef, nx0) - dx = Array{Taylor1{U}}(undef, nx0) - x .= Taylor1.( x0, order ) - _dv = Array{TaylorN{Taylor1{U}}}(undef, dof) # Allocation - cache = init_cache(LyapunovSpectrumCache, Val(false), t0, x, maxsteps) + cache = init_cache(LyapunovSpectrumCache, Val(false), t0, q0, maxsteps, order) # If user does not provide Jacobian, check number of TaylorN variables and initialize _dv if isa(jacobianfunc!, Nothing) @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) + cache.dvars[ind] = one(cache.x[1])*TaylorN(Taylor1{U}, ind, order=1) end end # Determine if specialized jetcoeffs! method exists - parse_eqs, rv = _determine_parsing!(parse_eqs, f!, t, - view(x, 1:dof), view(dx, 1:dof), params) + parse_eqs, rv = _determine_parsing!(parse_eqs, f!, cache.t, + view(cache.x, 1:dof), view(cache.dx, 1:dof), params) - return _lyap_taylorinteg!(f!, t, x, dx, q0, t0, tmax, abstol, jt, _dv, rv, + return lyap_taylorinteg!(f!, q0, t0, tmax, abstol, rv, cache, 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}, +function lyap_taylorinteg!(f!, + q0::Array{U,1}, t0::T, tmax::T, abstol::T, rv::RetAlloc{Taylor1{U}}, cache::LyapunovSpectrumCache, params, jacobianfunc!; parse_eqs::Bool=true, maxsteps::Int=500) where {T<:Real, U<:Number} - @unpack tv, xv, xaux, x0, λ, λtsum, δx, dδx, jac, varsaux, QH, RH, aⱼ, qᵢ, vⱼ = cache + @unpack tv, xv, xaux, x0, λ, λtsum, δx, dδx, jac, varsaux, QH, RH, aⱼ, qᵢ, vⱼ, t, x, dx, jt, dvars = cache # Allocation order = get_order(t) @@ -308,7 +299,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!(Val(parse_eqs), f!, t, x, dx, xaux, δx, dδx, jac, abstol, _δv, + δt = lyap_taylorstep!(Val(parse_eqs), f!, t, x, dx, xaux, δx, dδx, jac, abstol, dvars, 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)) @@ -330,7 +321,10 @@ function _lyap_taylorinteg!(f!, t::Taylor1{T}, x::Array{Taylor1{U},1}, dx::Array for ind in eachindex(QH) @inbounds x0[dof+ind] = QH[ind] end - x .= Taylor1.( x0, order ) + @inbounds for i in eachindex(x0) + x[i][0] = x0[i] + TaylorSeries.zero!(dx[i], 0) + end if nsteps > maxsteps @warn(""" Maximum number of integration steps reached; exiting. @@ -347,48 +341,38 @@ function lyap_taylorinteg(f!, q0::Array{U,1}, trange::AbstractVector{T}, order::Int, abstol::T, params = nothing, jacobianfunc! = 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) - t = trange[1] + Taylor1( T, order ) - jt = Matrix{U}(I, dof, dof) - t = Taylor1(T, order) - x0 = vcat(q0, reshape(jt, dof*dof)) - nx0 = length(x0) - x = Array{Taylor1{U}}(undef, nx0) - dx = Array{Taylor1{U}}(undef, nx0) - x .= Taylor1.( x0, order ) - _dv = Array{TaylorN{Taylor1{U}}}(undef, dof) + + # Check if trange is increasingly or decreasingly sorted + @assert (issorted(trange) || + issorted(trange, rev=true)) "`trange` or `reverse(trange)` must be sorted" # Allocation - cache = init_cache(LyapunovSpectrumTRangeCache, Val(false), trange, x, maxsteps) + cache = init_cache(LyapunovSpectrumTRangeCache, Val(false), trange, q0, maxsteps, order) # If user does not provide Jacobian, check number of TaylorN variables and initialize _dv 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) + cache.dvars[ind] = one(cache.x[1])*TaylorN(Taylor1{U}, ind, order=1) end end # Determine if specialized jetcoeffs! method exists - parse_eqs, rv = _determine_parsing!(parse_eqs, f!, t, - view(x, 1:dof), view(dx, 1:dof), params) + parse_eqs, rv = _determine_parsing!(parse_eqs, f!, cache.t, + view(cache.x, 1:dof), view(cache.dx, 1:dof), params) - return _lyap_taylorinteg!(f!, t, x, dx, q0, trange, abstol, jt, _dv, rv, cache, + return lyap_taylorinteg!(f!, q0, trange, abstol, rv, cache, params, jacobianfunc!; parse_eqs, maxsteps) 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}, +function lyap_taylorinteg!(f!, + q0::Array{U,1}, trange::AbstractVector{T}, abstol::T, rv::RetAlloc{Taylor1{U}}, cache::LyapunovSpectrumTRangeCache, params, jacobianfunc!; parse_eqs::Bool=true, maxsteps::Int=500) where {T<:Real, U<:Number} - @unpack xv, xaux, x0, q1, λ, λtsum, δx, dδx, jac, varsaux, QH, RH, aⱼ, qᵢ, vⱼ = cache + @unpack xv, xaux, x0, q1, λ, λtsum, δx, dδx, jac, varsaux, QH, RH, aⱼ, qᵢ, vⱼ, t, x, dx, jt, dvars = cache # Auxiliary variables order = get_order(t) @@ -410,7 +394,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!(Val(parse_eqs), f!, t, x, dx, xaux, δx, dδx, jac, abstol, _δv, + δt = lyap_taylorstep!(Val(parse_eqs), f!, t, x, dx, xaux, δx, dδx, jac, abstol, dvars, 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)) @@ -458,7 +442,10 @@ function _lyap_taylorinteg!(f!, t::Taylor1{T}, x::Array{Taylor1{U},1}, dx::Array for ind in eachindex(QH) @inbounds x0[dof+ind] = QH[ind] end - x .= Taylor1.( x0, order ) + @inbounds for i in eachindex(x0) + x[i][0] = x0[i] + TaylorSeries.zero!(dx[i], 0) + end if nsteps > maxsteps @warn(""" Maximum number of integration steps reached; exiting. diff --git a/src/rootfinding.jl b/src/rootfinding.jl index e8e1c4f3..482ffba3 100644 --- a/src/rootfinding.jl +++ b/src/rootfinding.jl @@ -212,30 +212,22 @@ function taylorinteg(f!, g, q0::Array{U,1}, t0::T, tmax::T, @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) - x .= Taylor1.( q0, order ) - dx .= Taylor1.( zero.(q0), order ) - # Allocation - cache = init_cache(VectorCache, Val(dense), t0, x, maxsteps) + cache = init_cache(VectorCache, Val(dense), t0, q0, maxsteps, order) # Determine if specialized jetcoeffs! method exists - parse_eqs, rv = _determine_parsing!(parse_eqs, f!, t, x, dx, params) + parse_eqs, rv = _determine_parsing!(parse_eqs, f!, cache.t, cache.x, cache.dx, params) - return _taylorinteg!(Val(dense), f!, g, t, x, dx, q0, t0, tmax, abstol, rv, cache, params; + return taylorinteg!(Val(dense), f!, g, q0, t0, tmax, abstol, rv, cache, params; parse_eqs, maxsteps, eventorder, newtoniter, nrabstol) end -function _taylorinteg!(dense::Val{D}, f!, g, t::Taylor1{T}, x::Array{Taylor1{U},1}, dx::Array{Taylor1{U},1}, +function taylorinteg!(dense::Val{D}, f!, g, q0::Array{U,1}, t0::T, tmax::T, abstol::T, rv::RetAlloc{Taylor1{U}}, cache::VectorCache, params; parse_eqs::Bool=true, maxsteps::Int=500, eventorder::Int=0, newtoniter::Int=10, nrabstol::T=eps(T)) where {T <: Real,U <: Number, D} - @unpack tv, xv, psol, xaux = cache + @unpack tv, xv, psol, xaux, t, x, dx = cache # Initial conditions order = get_order(t) @@ -305,32 +297,23 @@ function taylorinteg(f!, g, q0::Array{U,1}, trange::AbstractVector{T}, # 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) - x .= Taylor1.( q0, order ) - dx .= Taylor1.( zero.(q0), order ) + issorted(trange, rev=true)) "`trange` or `reverse(trange)` must be sorted" # Allocation - cache = init_cache(VectorTRangeCache, Val(false), trange, x, maxsteps) + cache = init_cache(VectorTRangeCache, Val(false), trange, q0, maxsteps, order) # Determine if specialized jetcoeffs! method exists - parse_eqs, rv = _determine_parsing!(parse_eqs, f!, t, x, dx, params) + parse_eqs, rv = _determine_parsing!(parse_eqs, f!, cache.t, cache.x, cache.dx, params) - return _taylorinteg!(f!, g, t, x, dx, q0, trange, abstol, rv, cache, params; + return taylorinteg!(f!, g, q0, trange, abstol, rv, cache, 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}, +function taylorinteg!(f!, g, q0::Array{U,1}, trange::AbstractVector{T}, abstol::T, rv::RetAlloc{Taylor1{U}}, cache::VectorTRangeCache, params; parse_eqs::Bool=true, maxsteps::Int=500, eventorder::Int=0, newtoniter::Int=10, nrabstol::T=eps(T)) where {T <: Real,U <: Number} - @unpack tv, xv, xaux, x0, x1 = cache + @unpack tv, xv, xaux, x0, x1, t, x, dx = cache # Initial conditions @inbounds t0, t1, tmax = trange[1], trange[2], trange[end] diff --git a/test/taylorize.jl b/test/taylorize.jl index 393dbdba..5b96bfac 100644 --- a/test/taylorize.jl +++ b/test/taylorize.jl @@ -1330,26 +1330,13 @@ import Logging: Warn _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) # Allocation - _cache_true = TaylorIntegration.init_cache(TaylorIntegration.VectorCache, Val(true), _t0, _x, _maxsteps) - _cache_false = TaylorIntegration.init_cache(TaylorIntegration.VectorCache, Val(false), _t0, _x, _maxsteps) - solTN = @inferred TaylorIntegration._taylorinteg!(Val(true), kepler1!, _t, _x, _dx, _q0, _t0, _tmax, _abstol, __rv, _cache_true, _params; parse_eqs=__parse_eqs, maxsteps=_maxsteps) - solTN2 = @inferred TaylorIntegration._taylorinteg!(Val(false), kepler1!, _t, _x, _dx, _q0, _t0, _tmax, _abstol, __rv, _cache_false, _params; parse_eqs=__parse_eqs, maxsteps=_maxsteps) + _cache_true = TaylorIntegration.init_cache(TaylorIntegration.VectorCache, Val(true), _t0, _q0, _maxsteps, _order) + _cache_false = TaylorIntegration.init_cache(TaylorIntegration.VectorCache, Val(false), _t0, _q0, _maxsteps, _order) + # Determine if specialized jetcoeffs! method exists + __parse_eqs, __rv = TaylorIntegration._determine_parsing!(true, kepler1!, _cache_true.t, _cache_true.x, _cache_true.dx, _params); + solTN = @inferred TaylorIntegration.taylorinteg!(Val(true), kepler1!, _q0, _t0, _tmax, _abstol, __rv, _cache_true, _params; parse_eqs=__parse_eqs, maxsteps=_maxsteps) + solTN2 = @inferred TaylorIntegration.taylorinteg!(Val(false), kepler1!, _q0, _t0, _tmax, _abstol, __rv, _cache_false, _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