MEGNO time evolution
alelatt opened this issue · 17 comments
Hi, I found a strange behaviour in the time evolution of MEGNO and I'd like to check if it should be expected.
What you're looking at is the (backwards) time evolution of MEGNO for a three-planet system with different starting conditions taken from within uncertainties.
These planets are integrated for 10^7 orbits of the outer planet (much as in Gajdos, Vanko, 2022). IAS15 is used as the integrator.
Plots 0, 2, 3, 6, 7 display chaotic behaviour. My question lies in the other ones: is the deviation from 2 to be expected?
This behaviour is also present when integrating forward in time, so I'd think that this has nothing to do with the time direction of choice.
The only information I managed to gather in the literature is Cincotta, Giordano, Theory and applications of the Mean Exponential Growth Factor of Nearby Orbits(MEGNO) method, Lecture Notes in Physics, 915, 93 (2016) in which it's briefly said (first half of page 7) that in the limit that the orbits are close to a stable periodic orbit the mean of Y should go to zero, but I haven't found any other information.
It's worth noting that I have done the same for another system and the same behaviour appears. What puzzles me is that I haven't found even a single realization (after about 100 different realizations across both systems) in which I observe a flat MEGNO~2.
Should I look for a bug in my code or is this behaviour in some way expected when integrating for so long?
Thank you in advance
This is not a direct answer to your question, but if you look at this example, you'll see that there is a region where MEGNO roughly approaches 2. But you'll also see the uncertainty (the color is not a uniform green). In general, I find it hard to interpret the results by looking at the MEGNO versus time curves. I think you'll probably gain more insight if you plot a slice of the parameter space like in the example.
Thank you for the answer.
The plot in the example is more or less what I was aiming for as an end result.
I tried plotting MEGNO vs time only after seeing several results with MEGNO < 2 to see what was happening and found that unexpected behaviour, so the plots I've shown you are definitely not the aim for my project.
Do you think that integrating for so many orbits might be part of the problem?
Yes. You need to choose the integration time with care. I don't have a mathematical proof, but my intuition is that you roughly want to integrate for a Laypunov timescale, or a couple of hundred secular/resonant timescales. If the integration time is many order of magnitude shorter, then chaos won't show up because the motion is dominated by periodic orbits. Anything many orders of magnitude longer then you'll notice that any system you're integrating is probably chaotic in the strictest sense (if not, then at some point finite floating point precision will be an issue). So from a practical point of view, you don't really want to measure if the system is chaotic or not (we're not mathematicians), but how chaotic it is compared to another. Hence the importance of choosing the right integration time...
Thank you very much, this is a very clear explaination.
I'll mark this question as closed for now. Should I find anything of relevance I'll reopen it.
Thanks again for the help
You're welcome!
Hi, sorry to reopen the question but I have some additional information.
I'm backwards integrating a much simpler and more stable system (2 coplanar planets with nearly no eccentricity), where I'd expect MEGNO ~ 2.
Integrating with IAS15 I get the same behaviour I've shown with the other system:
I then tried integrating with WHFast and I get exactly what I expect:
Please note that, while I am sampling from within the uncertainties for each orbital parameter to get the initial conditions (so no two systems are exactly the same in these examples), I'm only showing you a small subset of all the realizations I have integrated, so I feel that the two behaviours should be statistical and not due to good or bad luck in the sampling. All the integrations done with IAS15 show the same behaviour and same for all the integrations done using WHFast.
My question now is: could there be an issue in the way MEGNO is implemented in IAS15?
What are the units of the x axis?
Years, I've integrated backwards for 10^7 orbits of the outer planet.
Thanks. Hm. I'm not sure what's going on. It could be that the variable timestep leads to the discrepancy. If so, I'm not sure if this is a bug or maybe expected behaviour. In any case you could try out if it makes a difference if you set the IAS15 timestep to some reasonable value (maybe the same as WHFast) and turn off adaptive timesteps. Set sim.ri_ias15.epsilon = 0
Sure, I'll try. I'll get back to you tomorrow, thanks again!
I have tried using IAS15 with the same timestep as WHFast. The run hasn't completed since it exceeded the computation time I had available, but it's gone far enough to notice the same behaviour as shown previously, so I'd say that it shouldn't be due to the adaptive timestep.
Thanks. I'll see if I can reproduce the issue.
Here's a short notebook integrating a regular system with both IAS15 and WHFast. MEGNO is going towards 2 in both cases as expected. I'm not sure what could be different in your case. You'd need to post a similar short example that reproduces the issue if you want me to further investigate.
Thank you for the reply. I have looked at your plot and you're right, these seem to be converging to 2. The issue though is that you have integrated for ~10^5 orbits, and on that scale even my values are around 2 for both IAS and WHFast. You start seeing the wierd behaviour after at least 10^6 orbits for both systems I've looked at. I'll try to integrate the system in your notebook for longer and see what happens, otherwise I'll try to give you a simpler working example.
Hi, sorry for not replying for long.
To check that the results I have are independent of my code (which is a bit more convoluted than something as in you notebook) I have copied the initial conditions of a run in this way:
import rebound
import numpy as np
sim = rebound.Simulation()
sim.units = ('AU', 'MEarth', 'yr')
init_conf = np.loadtxt('init.txt')
sim.add(m = init_conf[0,0], x = init_conf[0,1], y = init_conf[0,2], z = init_conf[0,3], vx = init_conf[0,4], vy = init_conf[0,5], vz = init_conf[0,6])
sim.add(m = init_conf[1,0], x = init_conf[1,1], y = init_conf[1,2], z = init_conf[1,3], vx = init_conf[1,4], vy = init_conf[1,5], vz = init_conf[1,6])
sim.add(m = init_conf[2,0], x = init_conf[2,1], y = init_conf[2,2], z = init_conf[2,3], vx = init_conf[2,4], vy = init_conf[2,5], vz = init_conf[2,6])
sim.dt = -sim.particles[1].P/720
sim.init_megno()
times = np.linspace(0,-690000,1000)
megno = np.zeros(len(times))
for i, t in enumerate(times):
sim.integrate(t)
megno[i] = sim.megno()
res = np.array([times,megno])
np.savetxt(f'./ias.txt', res.T)
and same with WHFast
The initial conditions are in
init.txt
Would you say that replicating the initial conditions like this is enough for the simulation to be independent of the original one? And in that case, can you run this code and check?
Thanks for the code. However, this simulation integrates to almost 1e8 orbits (and takes about 8 hours to run). This is too long, both for debugging any potential issues, and for interpreting MEGNO.
My conclusion remains the same. I don't think there is any bug. WHFast uses a symplectic tangent map which is probably better at the long term integration of variational equations used within MEGNO (in terms of numerical errorS). However, the main issue here is the long integration time. This is not what a fast chaos indicator like MEGNO was designed to do. I'd advise to clarify what important timescales are present in your problem (orbital period, resonant, secular, Lyapunov, ..). Instead of using a fast chaos indicator, you might also just want to run an ensemble of simulations and analyze the statistical differences after long integrations. This is what people typically do for this kind of simulation. See e.g. https://arxiv.org/abs/2303.05567
Thank you very much for your answer and for the reference!