Wrong lyapunov exponent of unstable systems when the solver works adaptively

[mapi1]

Describe the bug I was testing the lyapunov() function for the simple exponential function dx/dt = ax as I know it has to be a in this case. Tough, just using the function gives some interesting results whenT > 1 is set:

Debugging a little I found that the problem might be related to the adaptiveness of the solver leading to large time steps and therefore wrong results after rescaling took place. Setting diffeq = (alg = Tsit5(), adaptive = false, dt = 1.0) gives correct results:

Maybe it is not so relevant for more complex problems, but would it be reasonable to reset the step size after rescaling took place?

Minimal Working Example

using DynamicalSystems
using ChaosTools
using OrdinaryDiffEq
using GLMakie

function sys(du, u, p, t)
    a = p[1]
    du[1] = a * u[1]
end

function simple(u0=[0.0]; a=0.0)
    # return CoupledODEs(sys, u0, [a]) # wrong
    return CoupledODEs(sys, u0, [a], diffeq = (alg = Tsit5(), adaptive = false, dt = 1.0)) # correct
end

ds = simple(; a = -1.0)

a_values = -1.0:0.025:1.0
λs = zeros(length(a_values), 4)

for i in eachindex(a_values)
    system = deepcopy(ds)
    set_parameter!(system, 1, a_values[i])
    λs[i, 1] = lyapunov(system, 1; Ttr = 0, u0 = 0.0)
    λs[i, 2] = lyapunov(system, 2; Ttr = 0, u0 = 0.0)
    λs[i, 3] = lyapunov(system, 50; Ttr = 0, u0 = 0.0)
    λs[i, 4] = lyapunov(system, 500; Ttr = 0, u0 = 0.0)
end

fig = Figure(fontsize = 25, dpi = 400)
ax = Axis(fig[1,1]; xlabel = L"a", ylabel = L"\lambda", title = L"\dot{x} = ax")
lins = [lines!(ax, a_values, λs[:, j], linewidth = 3) for j in 1:size(λs, 2)]
Legend(fig[1, 2], lins, ["T = $j" for j in [1,2,50,500]])
fig

Package versions

With Julia v1.9.0:

  [608a59af] ChaosTools v3.0.2
⌃ [61744808] DynamicalSystems v3.0.0
⌃ [e9467ef8] GLMakie v0.8.2
⌃ [1dea7af3] OrdinaryDiffEq v6.41.0

[Datseris]

Maybe it is not so relevant for more complex problems, but would it be reasonable to reset the step size after rescaling took place?

That sounds like a great suggestion, and in fact it sounds to me that the step size should be rescaled any time a new state is set to the integrator.

[Datseris]

I've known of the problem you quote here for unstable systems for quite some time, but haven't done anything about it yet...

[mapi1]

So, I tried to look deeper into this and have a try at fixing this by rescaling the step size. There are several ways to achieve this, but I settled on calling reinit!() instead of set_state!(). The thing is, it did not bring a significant improvement, sticking with the above example:

I added some debug printing and found that step!() rapidly increased dt, with much larger increments than when I solve the same problem via solve(ODEProblem(parallel_rule(ds, [0.0, 1e-9])..., (1,111), [a]), Tsit5()). The increments are actually the same as for solving solve(ODEProblem(sys, [0.0], (1,111), [a]), Tsit5()) which has the trivial solution 0 for all t. So it seems like the adaptivity is only decided on the first variable of the ParallelDynamicalSystem which I could further trace down to the DynamicalSystemsBase.matrixnorm that is used as internal norm. Removing it from the ParallelDynamicalSystem(ds::CoreDynamicalSystem, states) definition and letting DE use its default gave substantially better results:

The remaining error is related to error control and adjusting that finally gave correct results, even if the dt reset is removed and the original lyapunov() function is used:

I did not yet figure out what's wrong with the norm, but I wanted to write down my debugging results somewhere. Maybe this issue can be closed then and a new one in DynamicalSystemsBase should be opened