In the documentation it is specified the following about the Jacobian:

Jacobian

The keyword J provides the Jacobian function. It must be a Julia function in the same form as f, the dynamic_rule. Specifically, J(u, p, n) -> M::SMatrix for the out-of-place version or J(M, u, p, n) for the in-place version acting in-place on M. in both cases M is a matrix whose columns are the deviation vectors.

The bold part does not look consistent with what I can see in the source code and in the examples. Function J seems to return the Jacobian matrix J_f(x). The time evolution of the deviation vectors is then obtained multiplying the previous vectors by this matrix:

\frac{dY}{dt} = J_f(x)\cdot Y

If J was what the documentation says, then it would be possible to provide the evolution in tangent space as any type of function. I was planning to exploit this to speed up my computations. However, I donâ€™t think that this tangent evolution works as it is explained in the docs.

Could anyone clarify whether my interpretation is right? If so, the docs should be modifiedâ€¦

What do you mean? What do you see in the docs that is inconsistent?

I donâ€™t understand what you are proposing here. Can you give an example? What does â€śany type of functionâ€ť mean? You hope to provide a function f(Y, u, p, t) that returns directly dY/dt ? If so, why?

Evolution in tangent space works as explained in the docs, which is the equation you state in your post. I am very confident on that. Why do you think it is wrong?

Oh sh*t you are right, the docstring doesnâ€™t make sense. Of course the Jacobian function should return the Jacobian matrix, not a matrix with columns the deviation vectorsâ€¦ The jacobian matrix multiplies the latter matrix. that docstirng statement is indeed dead wrong Do you mind doing a quick PR to fix it?

coincidentally, two weeks ago, I was nearly about to re-write the continuous time version for Tangent systems but I decided I better focus on my job insteadâ€¦ But a rework may also allow more flexibility in how to specify the Jacobian.

Just did it. I hope I did it well, it is my first PR everâ€¦

I donâ€™t understand what you are proposing here. Can you give an example? What does â€śany type of functionâ€ť mean? You hope to provide a function f(Y, u, p, t) that returns directly dY/dt ? If so, why?

Exactly, I am thinking on providing a function that gives the evolution of Y directly, like what you said. The idea would be to speed up the computations. I am computing the Lyapunov spectra of a very large system in which the Jacobian is very sparse, but sparse does not improve the performance. I was hoping that something hard coded would help, as I can exploit some properties of the system and avoid multiplying by so many zeros.

I understand that this is a very particular case, and the performance might end up not improving. Anyhow, I will try to do this by redefining TangentDynamicalSystem so that tangentrule is my new function.

Thatâ€™s a good point, I did not benchmark in detail before. I just ran a simulation using Juliaâ€™s profiler and it does look like step! is the bottleneck (I attach a screenshot of ProfView).

Ok, for the moment I will try to speed things up by tweaking the implementation in tangent space, but in the future I might want to have a deeper look.