Collocation for finding periodic orbits of ODEs

[JonasKoziorek]

This blogpost shows working collocation algorithm for detection of periodic orbits of ODEs. It could be implemented in this package as well.

[Datseris]

cc @dawbarton who wrote the blog post and cc @rveltz because I know BifurcationKit.jl has the method (or a similar enough) implemented.

[rveltz]

I would say it is mostly done in BK and could be sent to NonlinearSolve. [I am quite busy to beginning of October]

[dawbarton]

The Fourier-based method in the post is easy to implement but it’s not great numerically. In certain systems (particularly forced mechanical systems) you can get numerical resonance effects because of the equispaced grid, which kills the convergence. Non-equispaced grids (e.g. Chebyshev) or piecewise polynomials with mesh adaptation (like AUTO or COCO use) are much better and more stable. I do have codes for these as well but someone would need to put a nicer interface on them. (That said, if the Fourier approach does work, the accuracy is very high.)

Ideally, someone would implement AUTO’s condensation method as that gives you stability information as well very cheaply (the parallel version that Bart Oldeman did would be even better).

What I would like to do in the future is generate composable discretisations, such that you can ask for things like a periodic orbit with a connecting orbit to an equilibrium (similar to COCO but not limited by MATLAB constraints). I did try this in ModelingToolkit.jl a few years ago, but found that the compilation times are excessively long.

[rveltz]

A few additionnal comments for collocation:

• condensation of parameters COP is implemented in BK (the serial one) as well as the Floquet exponents computation. Actually, I just noticed that COP is not mentioned in BK docs :(
• you need to have the analytical jacobian if you want the performances of Auto because the COP is basically half the time of the newton step. So MTK is used, it has to identify the block structure of the problem.
• having analytical jacobian allows you to compute periodic orbits in large dimensional systems but you need specific linear bordered solvers for this to be possible
• mesh adaptation also helps, especially close to homoclinic orbits ; this is also in BK. It is nice to have for continuing PO but I also wanted to use it to locate PO as well, like during a single newton computation. I did not do it.

All this (cop, analytical jacobian and mesh adaptation) is very tricky to implement. One tiny mistake and it allocates and kill performance. Or it is numerically unstable, or does not converge.

[Datseris]

How well does this method work for detecting unstable periodic orbits in chaotic systems?

[dawbarton]

How well does this method work for detecting unstable periodic orbits in chaotic systems?

If you’ve got a half decent initial guess for the orbit (or you’re tracking orbits through continuation, e.g., continuing through a sequence of period doubling bifurcations) then they are very good. However, if you don’t have prior knowledge and just want to find nearby orbits, then they aren’t so good. Sometimes you can get good results by dropping the accuracy right down (using just a handful of mesh points usually increases the radius of convergence; once converged you can refine the mesh again) but this is problem dependent.

mesh adaptation also helps, especially close to homoclinic orbits ; this is also in BK. It is nice to have for continuing PO but I also wanted to use it to locate PO as well, like during a single newton computation. I did not do it.

Mesh adaptation during a Newton iteration isn’t always recommended; it can create cycles in the iteration that stop convergence. I seem to remember that this is particularly a problem for PDE discretisations but it’s been ca. 15 years since I last looked at these details so my memory is hazy.

It’s great to hear that COP is in BK - I’ll have to take a look.

When you talk about analytical Jacobians, do you mean of the collocation scheme or of the underlying vector field? (I’m assuming the former since AUTO usually just uses finite differences for the vector field - most people don’t bother to provide derivatives.)

[rveltz]

I just mean that if yiu dont use the band structure of the PO jacobian you miss a lot of perfs.

I am sometimes shocked by the robustness and perfs of Auto given the choice it makes ( FD jacobians, Newton chord at iteration 3,...) even for very stiff problems.

I think a lot of perfs of BK is left on the table even if we are close to Auto (although it is very difficult to check). I maybe need a fresh eye ;P

[rveltz]

just mean that if yiu dont use the band structure of the PO jacobian you miss a lot of perfs.

Ok i understand the mistunderstanding. I meant that if you use MTK for encoding the PO probkem it would have to identify the band structure. But it likely will.

I do not rely on MTK for that because I want to treat problems that MTK cannot handle yet. But it wiyld be nice to see how it compares to my analytical jacobian. I bet MTK would be able to beat BK if it can use parallel jacobian filling

[dawbarton]

Ah, ok, I see what you mean. When I last tried, MTK exploits some of the structure - more than I expected but not as much as COP or a bordered solve. (Though I was focusing on flexible problem definition rather than performance at that point.)