Fix in Lyapunov spectrum caches

PerezHz · Jan 6, 2025 · d7d4769 · d7d4769
1 parent 617cafb
commit d7d4769
Show file tree

Hide file tree

Showing 2 changed files with 13 additions and 3 deletions.
diff --git a/src/integrator/cache.jl b/src/integrator/cache.jl
@@ -196,7 +196,9 @@ function init_cache_lyap(t0::T, q0::Vector{U}, maxsteps::Int, order::Int) where
     # Initialize the vector of Taylor1 expansions
     dof = length(q0)
     jt = Matrix{U}(I, dof, dof)
-    x0 = vcat(q0, reshape(jt, dof*dof))
+    x0 = similar(q0, dof + dof*dof)
+    x0[1:dof] .= q0
+    @views x0[dof+1:end] .= jt[:]
     t, x, dx = init_expansions(t0, x0, order)
     # Initialize cache
     nx0 = length(x0)
@@ -233,7 +235,9 @@ function init_cache_lyap(trange::AbstractVector{T}, q0::Vector{U}, maxsteps::Int
     t0 = trange[1]
     dof = length(q0)
     jt = Matrix{U}(I, dof, dof)
-    x0 = vcat(q0, reshape(jt, dof*dof))
+    x0 = similar(q0, dof + dof*dof)
+    x0[1:dof] .= q0
+    @views x0[dof+1:end] .= jt[:]
     t, x, dx = init_expansions(t0, x0, order)
     # Initialize cache
     nn = length(trange)
@@ -244,7 +248,7 @@ function init_cache_lyap(trange::AbstractVector{T}, q0::Vector{U}, maxsteps::Int
         Array{U}(undef, dof, nn),
         nothing,
         Array{Taylor1{U}}(undef, nx0),
-        similar(x0),
+        x0,
         similar(q0),
         Array{U}(undef, dof, nn),
         similar(q0),

diff --git a/src/lyapunovspectrum.jl b/src/lyapunovspectrum.jl
@@ -285,6 +285,9 @@ function lyap_taylorinteg!(f!,
 
     # Initial conditions
     sign_tstep = copysign(1, tmax-t0)
+    x0[1:dof] .= q0
+    @views x0[dof+1:end] .= jt[:]
+    update!(cache, t0, x0)
     t00 = t0
     tspan = zero(T)
     @inbounds tv[1] = t0
@@ -374,6 +377,9 @@ function lyap_taylorinteg!(f!,
     # Initial conditions
     @inbounds t0, t1, tmax = trange[1], trange[2], trange[end]
     sign_tstep = copysign(1, tmax-t0)
+    x0[1:dof] .= q0
+    @views x0[dof+1:end] .= jt[:]
+    update!(cache, t0, x0)
     t00 = t0
     tspan = zero(T)
     @inbounds for ind in eachindex(q0)