\documentclass[11pt]{article}
\usepackage[margin=1in]{geometry}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{graphicx}
\usepackage{booktabs}
\usepackage{tikz}
\usepackage{pgfplots}
\pgfplotsset{compat=1.18}
\usepackage[colorlinks=true, linkcolor=black, citecolor=black, urlcolor=blue]{hyperref}

\definecolor{iedgreen}{HTML}{1C4A34}
\definecolor{iedmid}{HTML}{2D7A53}
\definecolor{iedblue}{HTML}{2F5F96}
\definecolor{iedred}{HTML}{B04A2F}
\definecolor{iedgrey}{HTML}{7C8A82}

\title{How to Lose a Planet by Choosing the Wrong Integration Method}
\author{The Science Journaling Club}
\date{Volume 3, Issue 1 $\cdot$ Fall 2026 $\cdot$ Methods}

\begin{document}
\maketitle

\begin{abstract}
Simulating gravity is one equation applied over and over: acceleration equals
$-GM$ times the position vector, divided by the cube of the distance. We
integrated a single two-body Kepler orbit (eccentricity $0.6$, one hundred
steps per orbit, a resolution that looks perfectly reasonable) with three
standard stepping rules, forward Euler, velocity Verlet, and fourth-order
Runge--Kutta (RK4), and checked every trajectory against the exact solution
of Kepler's equation rather than against each other. At that resolution,
Euler's energy flips sign within the first orbit, from $-0.500$ to $+0.083$,
and never recovers. Velocity Verlet's energy error oscillates in a bounded
band between roughly $0\%$ and $2.99\%$ for the full $300$ orbits run, with
a secular slope consistent with zero to $0.0$ of its own standard error,
exactly as symplectic theory predicts, and its angular momentum is conserved
to $2.4\times10^{-14}$ relative error, $107$ machine epsilons, indistinguishable
from floating-point roundoff. RK4 drifts too, far more slowly and in the
opposite direction: it crosses $1\%$ energy loss only after $47.6$ orbits and
$5\%$ after $234.2$ orbits, a secular trend $21{,}915$ standard errors from
zero that quietly contracts the orbit rather than expanding it. Fitting the
measured convergence order across a range of timesteps gives $0.952$ for
Euler, $1.998$ for Verlet and $4.166$ for RK4 against textbook values of $1$,
$2$ and $4$; repeating the fit across $40$ Monte Carlo trials at random
eccentricities on $[0.10,0.60]$ gives $0.9866\pm0.0021$, $2.0015\pm0.0002$
and $4.1183\pm0.0051$, each a small, real, explainable disagreement with the
textbook exponent reported here rather than rounded away. Giving Euler four
times the timestep resolution of RK4, so both cost the same number of force
evaluations per orbit, does not rescue it: at $e=0.6$ it still loses $1\%$ of
the orbital energy within about a thousandth of an orbit, a stability
failure a bigger budget cannot buy back. We also flag one place our own
title oversells the result: at this configuration Euler's planet does not
fall into the sun. It gains energy and flies apart, the better-documented
and opposite failure mode; ``spiral into the sun'' describes RK4's slow
contraction far better than it describes Euler's fast escape. Master seed
20260901. Total run time 65.0 seconds.
\end{abstract}

\section{One equation, applied over and over}

Gravity between two bodies is not a complicated law: one mass pulls the
other straight toward itself, weakening as the square of the distance.
``Integrating'' the orbit means exactly what it sounds like: take the
current position and velocity, apply a small correction, repeat tens of
thousands of times. The correction rule is not unique. Different rules
disagree, sometimes wildly, about how much energy the orbit keeps by the
end, and a rule chosen badly makes a simulated planet do something no real
planet does: gain or lose energy for free, for reasons living entirely in
the arithmetic. The solar system itself is chaotic on timescales of tens of
millions of years, so real long-term integrations are run at very fine
resolution and cross-checked between independent codes \cite{laskar1989,
quinn1991}. Nothing in the two-body problem studied here is chaotic; the
failures below are a wrong, repeatable answer from a fixed method at a fixed
resolution, not sensitivity to initial conditions. This is a methods study.
The two-body problem has a closed-form solution, and every number below is
checked against that exact solution, never against another simulation.
Nothing here was observed; the computation is the whole experiment.

\section{Three ways to take a step}

Units are chosen so $GM=1$ and the semi-major axis $a=1$, making the orbital
period exactly $T=2\pi$. State is $(x,y,v_x,v_y)$ and acceleration is always
$\mathbf{a}=-\mathbf{r}/r^3$.

\textbf{Forward Euler.} $x_{n+1}=x_n+v_n\Delta t$, $v_{n+1}=v_n+a(x_n)\Delta
t$. One force evaluation per step. First order.

\textbf{Velocity Verlet.} A half-kick, a full drift, a second half-kick:
$v_{n+1/2}=v_n+\tfrac12 a(x_n)\Delta t$, $x_{n+1}=x_n+v_{n+1/2}\Delta t$,
$v_{n+1}=v_{n+1/2}+\tfrac12 a(x_{n+1})\Delta t$. One force evaluation per
step (the final acceleration is reused as the next step's first). Second
order and symplectic, meaning the discrete map preserves phase-space volume
exactly \cite{channell1990}; higher-order symplectic integrators exist by
composing several such substeps with carefully chosen weights
\cite{yoshida1990}. Algebraically identical, in this kick-drift-kick form, to
the leapfrog scheme of plasma and $N$-body physics \cite{hut1995}, and to
Verlet's original 1967 scheme \cite{verlet1967}; velocity Verlet itself was
named in the molecular-dynamics literature \cite{swope1982}.

\textbf{Fourth-order Runge--Kutta (RK4).} Four intermediate force
evaluations per step combined with weights $1,2,2,1$ give fourth-order
accuracy \cite{butcher2016}. Four times the cost of Euler or Verlet, step for
step, and, despite the higher order, not symplectic \cite{sanzserna1992}.

Cost is charged in force evaluations throughout: one per step for Euler and
Verlet, four per step for RK4.

\section{Checking the machine before trusting it}

Two validations ran before any orbit was trusted. \textbf{A:} Kepler's
equation $M=E-e\sin E$ was solved for the eccentric anomaly $E$ by Newton's
method \cite{danby1987}, checked over 740 combinations of mean anomaly and
eccentricity; worst residual $8.882\times10^{-16}$, at most 9 iterations.
\textbf{B:} the closed-form solution was sampled at 4{,}000 points across 7
orbits at $e=0.6$ and its energy and angular momentum recomputed
independently; worst relative deviation from the theoretical constants
($E_0=-GM/2a=-0.5$, $L_0=\sqrt{GMa(1-e^2)}=0.8$) was $4.441\times10^{-15}$
for energy and $5.551\times10^{-16}$ for angular momentum, both at the
double-precision floor. Both checks passed before a single stepping rule was
run.

\section{The headline number}

Start energy $E_0=-0.500000000000000$, exact by construction. After one
orbit at $\Delta t=T/100$, Euler's energy is $E_0+1.166583\times|E_0| =
-0.5+0.583=+0.083$: positive, unbound, and unrecovered 300 orbits later
($+0.0828$). Verlet's energy after orbit 1 is $-0.4865$ and after orbit 300
is $-0.4855$, oscillating between roughly $0\%$ and $2.99\%$ relative error
the entire time. RK4's energy creeps from $-0.500105$ after orbit 1 to
$-0.532241$ after orbit 300, a small steady contraction. Table
\ref{tab:headline} gives the orbit at which each method first crosses a
fixed energy-loss threshold, and whether the crossing is permanent or a
transient excursion that later recovers.

\begin{table}[htbp]
\centering
\caption{Orbit of first crossing of a relative energy-loss threshold,
$e=0.6$, $\Delta t=T/100$, 300 orbits run. ``Permanent'' means the trace
never falls back under half the threshold again; ``transient'' means it
does.}
\label{tab:headline}
\begin{tabular}{lrlrlr}
\toprule
method & orbit @ 1\% & character & orbit @ 5\% & character & orbit @ 10\% \\
\midrule
forward Euler & 0.00 & permanent & 0.00 & permanent & 0.00 (permanent) \\
velocity Verlet & 0.02 & transient & not reached & --- & not reached \\
RK4 & 47.64 & permanent & 234.19 & permanent & not reached \\
\bottomrule
\end{tabular}
\end{table}

\section{How bad does it get}

\begin{figure}[htbp]
\centering
\begin{tikzpicture}
\begin{axis}[
  width=13cm, height=6.5cm,
  xlabel={orbit number}, ylabel={relative energy error, velocity Verlet and RK4 (\%)},
  xmin=0, xmax=300, ymin=-8, ymax=4,
  grid=both, grid style={gray!18},
  legend pos=south west, legend cell align=left,
  legend style={font=\footnotesize, draw=gray!40},
  tick label style={font=\footnotesize},
  label style={font=\small}
]
\addplot[iedmid, thick, mark=*, mark size=1.6pt,
         error bars/.cd, y dir=both, y explicit]
 coordinates {
 (1,1.4946) +- (0,1.4946) (21,1.4946) +- (0,1.4946)
 (41,1.5370) +- (0,1.4507) (61,1.4967) +- (0,1.4924)
 (81,1.5206) +- (0,1.4677) (101,1.5029) +- (0,1.4860)
 (121,1.5081) +- (0,1.4807) (141,1.6092) +- (0,1.3796)
 (161,1.5132) +- (0,1.4755) (181,1.4996) +- (0,1.4894)
 (201,1.5274) +- (0,1.4607) (221,1.4952) +- (0,1.4939)
 (241,1.5455) +- (0,1.4420) (261,1.4950) +- (0,1.4942)
 (281,1.4950) +- (0,1.4942) (300,1.4990) +- (0,1.4901)
 };
\addlegendentry{velocity Verlet, min--max band per orbit}
\addplot[iedblue, thick, no markers] coordinates {
 (1,-0.0107) (21,-0.4470) (41,-0.8628) (61,-1.2774) (81,-1.7080)
 (101,-2.1241) (121,-2.5562) (141,-2.9871) (161,-3.4170) (181,-3.8459)
 (201,-4.2736) (221,-4.7004) (241,-5.1463) (261,-5.5709) (281,-6.0147)
 (300,-6.4372)
};
\addlegendentry{RK4}
\addplot[gray, dashed, thick, no markers] coordinates {(0,0) (300,0)};
\end{axis}
\end{tikzpicture}
\caption{Relative energy error against orbit number, $e=0.6$, $\Delta t=T/100$,
300 orbits, velocity Verlet and RK4 only (forward Euler is already at
$+116.6\%$ by orbit 1 and is omitted from this axis; see Table
\ref{tab:headline}). Verlet's error bars show the true minimum-to-maximum
range within each orbit, plotted at the orbit's midpoint value, so the band
width is the real within-orbit oscillation rather than an artifact of
sampling phase. RK4 is a single steadily descending line: the
ordinary-least-squares slope is $-2.146\times10^{-4}$ per orbit against a
standard error of $9.79\times10^{-9}$, 21{,}915 of its own standard errors
from zero.}
\label{fig:drift}
\end{figure}

Symplectic theory makes a second, sharper prediction for Verlet specifically:
its angular momentum should be conserved up to roundoff alone, because a
kick-drift-kick step applied to a central force leaves $L=xv_y-yv_x$
algebraically unchanged at every substep (the kick changes velocity by a
vector parallel to position, so $\mathbf{r}\times\Delta\mathbf{v}=0$; the
drift changes position by a vector parallel to velocity, so
$\Delta\mathbf{r}\times\mathbf{v}=0$) \cite{hairer2006,benettin1994}. Measured
over the same 300 orbits, Verlet's angular momentum drifted by at most
$2.387\times10^{-14}$ relative to its start, 107 machine epsilons, against a
naive roundoff random-walk estimate of $3.85\times10^{-14}$ over 30{,}000
steps, a ratio of $0.62$: consistent with pure floating-point noise.
Forward Euler's angular momentum drifted $1.48\times10^{16}$ times further
over the identical run, and RK4's $3.06\times10^{11}$ times further, both
from truncation error rather than roundoff.

\section{Does the error actually shrink the way theory says}

\begin{figure}[htbp]
\centering
\begin{tikzpicture}
\begin{axis}[
  width=13cm, height=8cm,
  xmode=log, ymode=log,
  xlabel={timestep $\Delta t$}, ylabel={position error after one orbit},
  grid=both, grid style={gray!18},
  legend pos=north west, legend cell align=left,
  legend style={font=\footnotesize, draw=gray!40},
  tick label style={font=\footnotesize},
  label style={font=\small}
]
\addplot[iedred, thick, mark=*, mark size=1.6pt] coordinates {
 (0.003142,1.054000) (0.001571,0.610100) (0.000785,0.326400)
 (0.000393,0.166700) (0.000196,0.083740) (0.000098,0.041900)
 (0.000049,0.020950)
};
\addlegendentry{Euler, slope 0.952 (theory 1)}
\addplot[iedmid, thick, mark=square*, mark size=1.6pt] coordinates {
 (0.062832,0.527000) (0.031416,0.136100) (0.015708,0.033820)
 (0.007854,0.008433) (0.003927,0.002107) (0.001963,0.000527)
 (0.000982,0.000132)
};
\addlegendentry{Verlet, slope 1.998 (theory 2)}
\addplot[iedblue, thick, mark=triangle*, mark size=2pt] coordinates {
 (0.062832,0.004006) (0.031416,0.0001824) (0.015708,0.0000093)
 (0.007854,0.00000052) (0.003927,0.0000000306) (0.001963,0.00000000185)
 (0.000982,0.00000000011)
};
\addlegendentry{RK4, slope 4.166 (theory 4)}
\end{axis}
\end{tikzpicture}
\caption{Position error after exactly one orbit against timestep, $e=0.6$,
log-log axes, least-squares slope quoted per series. Each method uses its
own range of step counts, chosen so the fit sits in that method's asymptotic
regime rather than a shared range that would understate Euler's order
(\S\ref{sec:asymptotic}).}
\label{fig:convergence}
\end{figure}

\label{sec:asymptotic}
Each method needs its own range of step counts before the theoretical
scaling of Table \ref{tab:conv} shows up cleanly, and that range differs by
two orders of magnitude between methods: Verlet and RK4 settle into slopes
of 2 and 4 by a few hundred steps per orbit, but Euler's locally fitted
slope at a few thousand steps per orbit is still only about 0.8, reaching
1.000 only in the tens of thousands of steps per orbit, several thousand
times finer than the 100-step configuration used everywhere else in this
article. Butcher's order conditions \cite{butcher2016} describe the limit
$\Delta t\to0$; they say nothing about how close a given problem must get to
that limit before it is a useful description.

\begin{table}[htbp]
\centering
\caption{Fitted convergence order, single orbit at $e=0.6$ versus 40 Monte
Carlo trials at eccentricities drawn uniformly on $[0.10,0.60]$, against
theory.}
\label{tab:conv}
\begin{tabular}{lrrrrr}
\toprule
method & theory & single-orbit fit & difference & Monte Carlo mean $\pm$ SE ($n=40$) & $z$ \\
\midrule
forward Euler & 1 & 0.952 & $-0.048$ & $0.9866\pm0.0021$ & $-6.31$ \\
velocity Verlet & 2 & 1.998 & $-0.002$ & $2.0015\pm0.0002$ & $+7.72$ \\
RK4 & 4 & 4.166 & $+0.166$ & $4.1183\pm0.0051$ & $+23.09$ \\
\bottomrule
\end{tabular}
\end{table}

\section{Is this just one lucky orbit?}

The single-eccentricity fit was repeated 40 times at eccentricities drawn
independently and uniformly from $[0.10,0.60]$, seeded once from the master
seed via numpy's SeedSequence. The Monte Carlo means (Table \ref{tab:conv})
agree with the single-orbit numbers to within two percent for every method,
so $e=0.6$ was not a fluke in either direction; the larger sample instead
buys a standard error small enough to make RK4's residual disagreement with
theory, $4.118$ against $4$, impossible to shrug off at $23.1$ of its own
standard errors. We attribute this to a genuine, small, higher-order
correction term in RK4's local truncation error that has not fully died away
even at the finest step counts used, the same phenomenon documented for
Euler above at much larger scale, and we report the gap rather than narrow
the step-count range until it disappears.

\section{Buying back stability}

Comparing methods at the same timestep compares them at different costs,
since RK4 spends four force evaluations per step against Euler and Verlet's
one. Table \ref{tab:cost} instead holds cost fixed: RK4 at 60 steps per
orbit (240 force evaluations per orbit) against Euler and Verlet at 240
steps per orbit each, four times finer a timestep, for the identical spend.

\begin{table}[htbp]
\centering
\caption{Orbits survived before a 1\% energy loss at matched computational
cost (240 force evaluations per orbit for all three methods), across three
eccentricities.}
\label{tab:cost}
\begin{tabular}{lllrl}
\toprule
$e$ & method & steps/orbit & orbit @ 1\% loss & character \\
\midrule
0.60 & forward Euler & 240 & 0.001 & permanent \\
0.60 & velocity Verlet & 240 & $>$1500 & not reached \\
0.60 & RK4 & 60 & 3.950 & permanent \\
0.20 & forward Euler & 240 & 0.011 & permanent \\
0.20 & velocity Verlet & 240 & $>$1500 & not reached \\
0.20 & RK4 & 60 & 1480.181 & permanent \\
0.10 & forward Euler & 240 & 0.019 & permanent \\
0.10 & velocity Verlet & 240 & $>$1500 & not reached \\
0.10 & RK4 & 60 & $>$1500 & not reached \\
\bottomrule
\end{tabular}
\end{table}

Four times the timestep resolution does not buy Euler a longer life at
$e=0.6$ or $e=0.2$; at matched timestep (all three at 60 steps per orbit,
$e=0.6$) Euler crosses 1\% at orbit 0.0002, and at matched cost it crosses at
orbit 0.001, five times longer measured in a unit where five times longer
than instantaneous is still instantaneous for every practical purpose. This
is an accuracy-versus-stability distinction: a bigger budget fixes the
first kind of problem and not the second, and at this eccentricity Euler's
failure is squarely the second kind. The eccentricity dependence in Table
\ref{tab:cost} is not incidental: at $e=0.1$, RK4 does not reach 1\% loss
inside 1{,}500 orbits at the same cost that gets it there in under 4 orbits
at $e=0.6$, because perihelion speed, not orbit count, sets the resolution a
fixed-step method actually needs.

\section{The strongest objection we could make}

Nobody serious runs forward Euler on an eccentric Kepler orbit, so what does
this demonstrate beyond a straw man? The answer is not that Euler is
secretly fine; it is that plenty of serious codes do not use velocity
Verlet either, for reasons unrelated to the failure documented here. The
Hermite integrator, fourth order and time-symmetric like RK4 but built
around continuously estimated higher derivatives, is standard in
direct-summation planet-formation codes because it tolerates close
encounters and a wide range of dynamical timescales that a fixed-step
symplectic method handles badly \cite{kokubo1998}; it is not symplectic
either, and by our own \S6 logic should show RK4-like secular drift at a
rate its designers judged an acceptable trade. Large cosmological $N$-body
codes such as GADGET-2 default to exactly the kick-drift-kick scheme used
here for its bounded-error property \cite{springel2005}, and modern
solar-system integrators such as REBOUND ship a Wisdom--Holman symplectic
mapping as their recommended default \cite{rein2012}, following the original
demonstration that symplectic mapping holds up over solar-system-length
integrations where naive methods do not \cite{wisdom1991}. The honest
picture is not ``symplectic good, everything else bad'' but ``know which
failure mode your problem can tolerate,'' since every fast method accepts
one.

\section{Where a different modelling choice would change the answer}

Three choices made above would each change a headline number if made
differently. The eccentricity, $0.6$ by our own pedagogical pick, sets how
fast failure arrives (Table \ref{tab:cost}: RK4 survives 3.95 orbits at
$e=0.6$ and over 1{,}500 at $e=0.1$ for the identical cost budget), so
quoting a crossing orbit without the eccentricity attached quotes a number
that depends on an undefended choice. The step-count range used for the
convergence fit (\S\ref{sec:asymptotic}) would have reported Euler's order
as $0.8$ rather than $1.0$ had the fit used the same few-hundred-step range
that works for Verlet and RK4, a natural shortcut that would silently
understate Euler's asymptotic order without any buggy code. And defining
cost in force evaluations rather than wall-clock seconds sidesteps
implementation quality entirely; a poorly vectorised RK4 against a
hand-tuned Euler could reverse every ranking in Table \ref{tab:cost} on
wall-clock time alone, which is exactly why cost was charged the way it was.

\section{Reproducing this}

Everything in this article comes from one file and one command:
\begin{verbatim}
python integrator-energy-drift.py > integrator-energy-drift-output.txt
\end{verbatim}
Python 3.12 and numpy are required; no other package is imported. Our run
took $65.0$ seconds on numpy 2.4.2 and Python 3.12.3, well inside a
ten-minute budget. Master seed 20260901 governs the Monte Carlo draw in
\S7 and nowhere else; every other number in this article is fully
deterministic and reproduces exactly on a rerun. An interactive version of
the same three integrators, at the same default eccentricity and timestep,
is available alongside this article.

\begin{thebibliography}{99}

\bibitem{verlet1967}
Verlet, L. (1967). Computer ``experiments'' on classical fluids. I.
Thermodynamical properties of Lennard-Jones molecules. \emph{Physical
Review} \textbf{159}, 98--103.
\href{https://doi.org/10.1103/PhysRev.159.98}{doi:10.1103/PhysRev.159.98}

\bibitem{swope1982}
Swope, W. C., Andersen, H. C., Berens, P. H. and Wilson, K. R. (1982). A
computer simulation method for the calculation of equilibrium constants for
the formation of physical clusters of molecules: application to small water
clusters. \emph{Journal of Chemical Physics} \textbf{76}, 637--649.
\href{https://doi.org/10.1063/1.442716}{doi:10.1063/1.442716}

\bibitem{wisdom1991}
Wisdom, J. and Holman, M. (1991). Symplectic maps for the $n$-body problem.
\emph{The Astronomical Journal} \textbf{102}, 1528--1538.
\href{https://doi.org/10.1086/115978}{doi:10.1086/115978}

\bibitem{yoshida1990}
Yoshida, H. (1990). Construction of higher order symplectic integrators.
\emph{Physics Letters A} \textbf{150}, 262--268.
\href{https://doi.org/10.1016/0375-9601(90)90092-3}{doi:10.1016/0375-9601(90)90092-3}

\bibitem{hut1995}
Hut, P., Makino, J. and McMillan, S. (1995). Building a better leapfrog.
\emph{The Astrophysical Journal} \textbf{443}, L93--L96.
\href{https://doi.org/10.1086/187844}{doi:10.1086/187844}

\bibitem{quinn1991}
Quinn, T. R., Tremaine, S. and Duncan, M. (1991). A three million year
integration of the earth's orbit. \emph{The Astronomical Journal}
\textbf{101}, 2287--2305.
\href{https://doi.org/10.1086/115850}{doi:10.1086/115850}

\bibitem{channell1990}
Channell, P. J. and Scovel, C. (1990). Symplectic integration of Hamiltonian
systems. \emph{Nonlinearity} \textbf{3}, 231--259.
\href{https://doi.org/10.1088/0951-7715/3/2/001}{doi:10.1088/0951-7715/3/2/001}

\bibitem{sanzserna1992}
Sanz-Serna, J. M. (1992). Symplectic integrators for Hamiltonian problems:
an overview. \emph{Acta Numerica} \textbf{1}, 243--286.
\href{https://doi.org/10.1017/S0962492900002282}{doi:10.1017/S0962492900002282}

\bibitem{springel2005}
Springel, V. (2005). The cosmological simulation code GADGET-2.
\emph{Monthly Notices of the Royal Astronomical Society} \textbf{364},
1105--1134.
\href{https://doi.org/10.1111/j.1365-2966.2005.09655.x}{doi:10.1111/j.1365-2966.2005.09655.x}

\bibitem{benettin1994}
Benettin, G. and Giorgilli, A. (1994). On the Hamiltonian interpolation of
near-to-the-identity symplectic mappings with application to symplectic
integration algorithms. \emph{Journal of Statistical Physics} \textbf{74},
1117--1143.
\href{https://doi.org/10.1007/BF02188219}{doi:10.1007/BF02188219}

\bibitem{hairer2006}
Hairer, E., Lubich, C. and Wanner, G. (2006). \emph{Geometric Numerical
Integration: Structure-Preserving Algorithms for Ordinary Differential
Equations}, 2nd edition. Springer-Verlag, Berlin. ISBN 978-3-540-30666-5.

\bibitem{danby1987}
Danby, J. M. A. (1987). The solution of Kepler's equation, III.
\emph{Celestial Mechanics} \textbf{40}, 303--312.
\href{https://doi.org/10.1007/BF01235847}{doi:10.1007/BF01235847}

\bibitem{rein2012}
Rein, H. and Liu, S.-F. (2012). REBOUND: an open-source multi-purpose
$N$-body code for collisional dynamics. \emph{Astronomy \& Astrophysics}
\textbf{537}, A128.
\href{https://doi.org/10.1051/0004-6361/201118085}{doi:10.1051/0004-6361/201118085}

\bibitem{kokubo1998}
Kokubo, E., Yoshinaga, K. and Makino, J. (1998). On a time-symmetric Hermite
integrator for planetary $N$-body simulation. \emph{Monthly Notices of the
Royal Astronomical Society} \textbf{297}, 1067--1072.
\href{https://doi.org/10.1046/j.1365-8711.1998.01581.x}{doi:10.1046/j.1365-8711.1998.01581.x}

\bibitem{laskar1989}
Laskar, J. (1989). A numerical experiment on the chaotic behaviour of the
Solar System. \emph{Nature} \textbf{338}, 237--238.
\href{https://doi.org/10.1038/338237a0}{doi:10.1038/338237a0}

\bibitem{butcher2016}
Butcher, J. C. (2016). \emph{Numerical Methods for Ordinary Differential
Equations}, 3rd edition. Wiley, Chichester.
\href{https://doi.org/10.1002/9781119121534}{doi:10.1002/9781119121534}

\end{thebibliography}

\end{document}