For a complete introduction to variational equations, please read the paper by Rein and Tamayo (2016).
For this tutorial, we work with a two planet system. We vary the initial semi-major axis $a$ of the outer planet. Because the planets interact with each other, the final $x$-position of the inner planet at the end of the simulation will depend on the initial semi-major axis of the outer planet. We run the simulation once for a fixed $a_0$ and then use first and second order variational equations to predict the final position of the outer planet for different $a$s in a neighbourhood of $a_0$.
To do that, let us first import REBOUND, numpy and matplotlib.
In [1]:
import rebound
import numpy as np
%matplotlib inline
import matplotlib;
import matplotlib.pyplot as plt
Before using variational equations, let us define a function that calculates the final position of the inner planet as a function of $a$ in the brute-force way:
In [2]:
def run_sim(a):
sim = rebound.Simulation()
sim.add(m=1.)
sim.add(primary=sim.particles[0],m=1e-3, a=1)
sim.add(primary=sim.particles[0],m=1e-3, a=a)
sim.integrate(2.*np.pi*10.)
return sim.particles[1].x
We'll use this function to create a list of true final positions to which we later compare our results.
In [3]:
N=400
x_exact = np.zeros((N))
a_grid = np.linspace(1.4,1.7,N)
for i,a in enumerate(a_grid):
x_exact[i] = run_sim(a)
Running a simulation with variational equations is very easy. We start by creating a simulation and add the three particles (the star and two planets) just as before. Note that the vary convenience function we use below only accepts heliocentric coordinates, so we explicitly tell REBOUND that the star is the primary when adding particles to the simulation.
We then add variational particles to the simulation. We vary one paramter ($a$) and thus need only one set of first order variational equaitions. The second order variational equations depend on the first order ones. Thus, when initializing them, one has to pass the set of first order variational equations using the 'first_order' parameter.
After adding a variation, one must always initialize it. We do this below with REBOUND's vary() convenience function, which makes varying orbital parameters particularly easy. Alternatively, one can also initialize the variational particles directly, e.g. using var_da.particles[1].x = 1. Note that variations are implemented as particles, but you they really represent derivatives of a particle's corrdinates with respect to some initial parameter. For more details, see Rein and Tamayo (2016).
The function below does all that and returns the final position of the inner planet, as well as the first and second derivatives of the position with respect to $a$.
In [4]:
def run_sim_var(a):
sim = rebound.Simulation()
sim.add(m=1.)
sim.add(primary=sim.particles[0],m=1e-3, a=1)
sim.add(primary=sim.particles[0],m=1e-3, a=a)
var_da = sim.add_variation()
var_dda = sim.add_variation(order=2, first_order=var_da)
var_da.vary(2, "a")
var_dda.vary(2, "a")
sim.integrate(2.*np.pi*10.)
return sim.particles[1].x, var_da.particles[1].x, var_dda.particles[1].x
We can now use the variational equations to predict the final position of the inner particle. Note that we only run one simulation, at $a_0=1.56$.
In [5]:
a_0 = 1.56
x, dxda, ddxdda = run_sim_var(a_0)
x_1st_order = np.zeros(N)
x_2nd_order = np.zeros(N)
for i,a in enumerate(a_grid):
x_1st_order[i] = x + (a-a_0)*dxda
x_2nd_order[i] = x + (a-a_0)*dxda + 0.5*(a-a_0)*(a-a_0)*ddxdda
In the figure below, we plot the final position as a function of the initial semi-major axis. The black line corresponds to the true final position as calculated by the brute-force approach. The dashed and dotted lines correspond to the approximations using first and second order variational equations. As one can see, the second order approximation is very accurate within a neighbourhood of $a_0$.
In [6]:
fig = plt.figure(figsize=(6,4))
ax = plt.subplot(111)
ax.set_xlim(a_grid[0],a_grid[-1])
ax.set_ylim(np.min(x_exact),np.max(x_exact)*1.01)
ax.set_xlabel("initial semi-major axis of the outer planet")
ax.set_ylabel("$x$ position of inner planet after 10 orbits")
ax.plot(a_grid, x_exact, "-", color="black", lw=2)
ax.plot(a_grid, x_1st_order, "--", color="green")
ax.plot(a_grid, x_2nd_order, ":", color="blue")
ax.plot(a_0, x, "ro",ms=10);
In [ ]: