Nonlinear vector field

[ghost]

Hi, I am working on a system of delay differential equations where the vector field is quite nonlinear. My main interest is calculating the largest Lyapunov Exponent as function of vector field parameters. For this purpose, I can use the algorithm itself because it automatically calculates the Jacobian matrix. Unfortunately, this Jacobian matrix is again nonlinear due to nonlinear character of the vector field and I obtain incorrect values considering the expected Lyapunov Exponents. My question is: how could I set the Jacobian matrix before hand instead of being dependent on the algorithm?

Thank you, Edmilson Roque.

[Wrzlprmft]

Right now, there is no interface to provide the Jacobians yourself and I guess that it would be rather tedious to implement. The easiest way to manually supply the Jacobians right now would probably be to monkey-patch the functions jitcdde._jac or jitcdde.tangent_vector_f (if needed). (Look at the source code to see how they are working and what kind of output they supply.)

That being said, non-linear Jacobians are no reason for the automatic routines to fail (in fact, this is the most common case), and even if they were, I would like to know why this happens so I can fix it. Another possibility you have to consider is that the mistake is located somewhere else entirely. Can you share your problem with me (either here or via mail) as well as what you think the Jacobians (for each delay) should look like, so I can have a look at it? If you do not want to, you can use the aforementioned function jitcdde._jac to obtain the Jacobians in JiTCDDE’s way and compare them with your results.

[ghost]

Well I should say your answer helped me to realize some way to fix my problem. On the other hand, it could be nice to discuss the problem which I am pursuing to solve.

Problem: somehow I am interested to measure stability of a given sync manifold. There is an option in jitcdde, jitcdde_restricted_lyap, which automatically provides such linear stability analysis to an entire coupled system. However, I did not understand properly how it works. Then I calculated by hand (similar procedure (canonical) developed here: https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.80.2109).

In the end of the day, after some calculations, we end up with a reduced system of differential equations. For sure, the paper that I'm sending you is ordinary and not delay differential equations, but for instance look at here: https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.92.074104 (Maybe everything is trivial to you but I needed to make my point).

I tried to solve such reduced system of delayed differential equations (part of the entire vector field is the Jacobian matrix of the individual dynamics for each oscillator) and I failed when I increased the coupling among the oscillators. I got frustrated. My question/suggestion would be: could you describe deeper how jitcdde_restricted_lyap works? I can not see how we guarantee transversal vectors of a given sync state which a priori we do not know. And a tiny suggestion: maybe implement the aforementioned procedure could optimize the code to someone. Then, what do you think to include a comment about it to people that do not know the procedure and still need to implement really large systems?

Anyway, Thank you so much. Edmilson Roque.

[Wrzlprmft]

I only have time for a brief answer. I will probably extend this tomorrow.

your answer helped me to realize some way to fix my problem.

I am glad to help. But please tell me if there is any problem with JiTCDDE.

However, I did not understand properly how it works.

I recently added a test for this that may also help to illustrate how to use this feature.

I can not see how we guarantee transversal vectors of a given sync state which a priori we do not know.

Well, you know what the completely synchronised state looks like: All corresponding dynamical variables have identical values. That makes it easy to write down a basis of the respective hyperplane containing all these states (the synchronisation manifold) and staying transversal to that is just a question of subtracting projections to that hyperplane. (Note that to do all of this with JiTCDDE you have to set your dynamics to be on that hyperplane.)

[Wrzlprmft]

Through your comments I got some ideas that (if proven reasonable) may cause me to redesign this feature soon. That being said, here is how it works as of now:

For simplicity, suppose first that you have three, one-dimensional oscillators and an ODE. Then the synchronisation manifold is described by y(0) == y(1) == y(2). A base vector comprising these states is (1,1,1) and thus for the transversal Lyapunov exponents, we are interested in the growth of tangent vectors (or separation functions) orthogonal to this. What we can thus do is to evolve the tangent vectors regularly and regularly remove any projections to (1,1,1) (this is the vector you would give as an argument to jitcode_lyap). The growth of the resulting tangent vector is thus described by the maximum transverse Lyapunov exponent. Note that for all of this to work, we have to be on the synchronisation manifold, which we can achieve by initialising our differential equations with identical states. Also note that if your oscillators have more dimensions, you get additional base vectors for each: for two dimensions, your base vectors would be either [ [1,0,1,0], [0,1,0,1] ] (if your primary order of dynamical variables is by oscillator) or [ [1,1,0,0], [0,0,1,1] ] (if your primary order of dynamical variables is by component).

Now, for the method used by JiTCDDE, states are not only identified by the value of dynamical variables, but also by their derivative. Thus in the above example you would have to specify separate basis vectors for the state and derivative, namely the vectors argument of JiTCDDE would be [ ( [1,1,1], [0,0,0] ), ( [0,0,0], [1,1,1] ) ]. Moreover, as the states also include the past, the scalar products needed for removing the projection become quite expensive. This can be addressed by transforming your differential equations such that your base vectors (of the synchronisation manifold) only have one non-zero component – JiTCDDE automatically recognises this and is much quicker in removing this component. You can do this by extending your basis of the synchronisation manifold to an orthogonal basis of the entire state space (without delay) and transform your differential equations to that basis. This is more or less what is done in the second paper you referenced.

I am aware that there are many inefficiencies in this process, e.g., as you use identical initial conditions for the oscillators, you integrate the same differential equation multiple times. I am now thinking of a way to account for this.

[Wrzlprmft]

I implemented a more easy, robust, and efficient way to calculate transversal Lyapunov exponents to synchronisation manifolds for JiTCODE. I will translate this feature to DDEs soon (and don’t expect any major problems), but until then you can take a look at the ODE version.

[Wrzlprmft]

I implemented a considerably better tool for transversal Lyapunov exponents in the above commit. Note that, for now, you need to install directly from GitHub for this and need the GitHub version of Symengine.py and JiTC*DE. For the documentation see docstrings. Here is an example for how to use this feature for ODEs, which translates directly to DDEs.

I would be very interested in your feedback.