<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Variational Equations
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.
Step1:  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
Step2:  We'll use this function to create a list of true final positions to which we later compare our results.
Step3:  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$.
Step4:  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$.
Step5:  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$.
<ASSISTANT_TASK:>
Python Code:
import rebound
import numpy as np
%matplotlib inline
import matplotlib;
import matplotlib.pyplot as plt
Explanation: Variational Equations
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.
End of explanation
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
Explanation: 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:
End of explanation
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)
Explanation: We'll use this function to create a list of true final positions to which we later compare our results.
End of explanation
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
Explanation: 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$.
End of explanation
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
Explanation: 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$.
End of explanation
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);
Explanation: 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$.
End of explanation
<END_TASK>