ias15 and Megno

[mstalport]

Hi,

I have some issues with the integrator ias15 when it is coupled with Megno. Indeed, there are some behaviors I don’t understand when I launch a series of identical simulations (i.e. same integration setups and same initial conditions) of an unstable system.

For some given initial conditions, I noticed different outcomes between the several simulations if the Megno calculation is included or not. Without the Megno, all simulations are stopped at the same output time, when the system gets unstable. When the Megno calculation is included however, the system’s survival time can be variable from one simulation to another. Would you have an idea about to which extent could the Megno calculation modify the integration?

Here below, I show a case example: a 2-planets unstable system.

import os 
import math
import rebound

NBpl = 2 # The number of planets in the system 
T = 250000 # The integration time 

for simu in range(5): # 5 is the number of identical simulations we will sequentially launch. 
    sim = rebound.Simulation()
    sim.add(m=1.0)
    sim.add(m=1e-3, a=0.5, e=0.15, omega=1.0)
    sim.add(primary=sim.particles[0], m=1e-3, a=0.8, e=0.2) 

    sim.move_to_com()
    dt = 1.0*2*math.pi / 100. # The integration timestep 
    sim.dt = dt
    sim.integrator = "ias15" 
    sim.min_dt = dt/10. 
    sim.exit_max_distance = 200. * 0.8
    sim.init_megno() # Comment this in order to skip Megno calculation 

    try:
        out_file = open('OutputTime_%01d.dat' % simu,'a') 
        for j in range(62500):
            t_output = 400*j*dt
            sim.integrate(t_output) 
            out_file.write("{0:.18e}".format(t_output/(2*math.pi)) + '\n') # Output time in years (every 4 years)  

        out_file.close()

############################ Comment that part to skip Megno calculation                              
        LyTime = sim.calculate_lyapunov()
        LyTime = math.log10(abs(1/LyTime))

        LyTime_file = open('T-LyTime_%01d.dat' % simu,'a') 
        LyTime_file.write("{0:.2e}".format(T) + ' ' + "{0:.18e}".format(LyTime) + '\n') 
        LyTime_file.close()
############################

    except rebound.Escape:
        pass

    del sim 

There is another behavior I don’t understand, probably linked to the former, regarding this coupling between ias15 and Megno. To see this, one needs to modify the initial conditions, replacing the value of ac (semi-major axis of the outer planet) from 0.8 to 0.7 AU. Without Megno, the system holds for around 28000 years. Including the Megno calculation however, some simulations will stay blocked. That means, the output file does not update anymore, but the CPU is still requisitioned at a rate of 100%. It is like if the integration froze. This problem never shows up when the Megno is not calculated, or when the integrated system is stable. Would you have an idea on why such a freezing is observed, and how to avoid it?

Many thanks in advance for your help!

Manu

[hannorein]

Hi Manu,

You're integrating a chaotic systems and tiny changes to the initial conditions can lead to a dramatic difference after a finite time.

If you include MEGNO, then it might change the integration because IAS15 may take different timesteps. The difference is only at the level of machine precision at each step, but the chaotic systems amplifies this exponentially.

So I don't think this is a bug, but expected. Try varying your initial conditions a tiny bit and look at the distribution of the results, not a single simulation. (If there is still a major difference between simulations which use MEGNO and those who don't, then I would conclude that there is a bug.)

If the simulation gets stuck, then most likely a close encounter happened which it is trying to resolve. To avoid this, set sim.exit_min_distance to some small distance. REBOUND then throws an exception if two particles come any closer. A reasonable value might be the physical size of planets.

I hope this helps!

Hanno

[mstalport]

Hi Hanno,

Many thanks for your helpful answers!

I didn't think about Megno impacting the choice of the timestep, but it makes sense to me actually (noting also that I hadn't such issues with whfast). I totally understand that this could therefore change the results of the integration (for an unstable system). It sounds good, I will try integrating nearby systems with and without the Megno, and verify that the distributions of the results are similar.

For the second issue, it was indeed a close encounter problem. The exception is raised when I include sim.exit_min_distance. I didn't think my simulations could get stuck because of a close encounter, because I had set sim.min_dt to a reasonable value. Therefore I was expecting that, in case of a close encounter, the simulation would not resolve well such an event but would not spend too much time on it neither. Anyway, many thanks again for your help!

Manu

[hannorein]

Note that you would want to set sim.ri_ias15.min_dt (it's an IAS15 specific setting).

[mstalport]

Great, thank you for the tip!