Implementation of Bennetin's algorithm for Lyapunov exponents

I want to calculate the Lyapunov indices for a system - problem is that I can’t get DynamicalSystems.jl to work (described here).

Therefore, I decided to do it a much more crude way by writing out the steps myself. Now, if I understand it right, the algorithm will compare the systems (with initial distance d_0) at a certain time. We note the distance in the 2 trajectories here as d_1, move the perturbed system back to d_0 distance from the evolved unperturbed trajectory, run the evolution for the next interval, repeat.

In order to do this, I decided to solve my original system for some tspan. Now, getting d_1 is easy enough - solve the same system with different initial conditions and subtract the 2 timeseries. So, for d_1, lets say I evolve the system for time, t_j. I am not sure what to do to get d_2, but I thought of the following way - solve the same system again but this time I use,

u_0 = [ x = sol[1, j] + d_0 ]
tspan = (t_j, t_k)

So, I choose the initial value of my variable this time as the value of the unperturbed variable at t_j plus the deviation, d_0. Now, I evolve this system from t_j to some other time t_k. Repeat the process.

I have not yet tried this (I will as soon as I get the chance), but I would like to know what the people here, who are familiar with Lyapunov exponents and so on, think of the validity and the viability of the process.

Thanks in advance!!

1 Like

Instead of reimplementing the well working Bennetin algorithm, from ChaosTools, better post the natural form of your system (not as you tried to define it here) to get help in redefining it as a right ds. Julia discourse renders a LaTeX text, hence you can write it down with all details.