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

```
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

*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

```
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 [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 [ ]:
```