Wrong lyapunov exponent of unstable systems when the solver works adaptively about chaostools.jl HOT 3 OPEN

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.

from chaostools.jl.

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

from chaostools.jl.

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

from chaostools.jl.