Kinds of orbits in two-body problem I have a numerical simulation for two bodies moving under gravitational interaction, and some of the outputs are not cohering with my understanding of the theory.
My understanding of the theory:
In the standard derivation of the analytic solution for the two body problem, we replace $\vec{r}_1$, $\vec{r}_2$ with $\vec{r}=\vec{r}_1-\vec{r}_2$ and c.o.m. vector $\vec{R}$. Then after some steps, we end up trying to solve the central force problem for a particle of mass $\mu$ (reduced mass), with position $\vec{r}$. Looking at the effective potential, we end up with 5 different cases - hyperbolic, parabolic, elliptical, circular, orbits, and direct infall. 
In particular, we are interested in the bound orbits (elliptical/spherical). If we now map back to our original coordinates, we have $\vec{r}_1=\vec{R}+\frac{m_2}{m_1+m_2}\vec{r}$ and $\vec{r}_2=\vec{R}-\frac{m_1}{m_1+m_2}\vec{r}$. Setting $\vec{R}$ fixed, the $\vec{r}$ components sweep out 'opposite' ellipses, with sizes that are scaled by the mass ratios. So the elliptical case corresponds to the two bodies drawing out stable ellipses, each with a focus at $\vec{R}$ the barycenter.
However, for certain input parameters, my code is producing an ellipse that is precessing around.

The code is using a standard Euler method, where we update the velocities/positions over small timesteps via $x(t+\delta t)=x(t)+v(t)\delta t$, $v(t+\delta t)=v(t)+a(t)\delta t$.
Here the initial conditions of my two bodies are set such that the net momentum is zero. Why is this precession occurring - is it likely due to numerical errors, or is this actually allowed and my understanding of the two-body solution missing something?
 A: The precession is just one symptom of a larger problem - inaccuracy of the calculation. The inaccuracy is because of two things:
1) your algorithm and really all usual integration algorithms is only an approximation to the mathematical model that is the differential equation in time; this approximation can be made better by using smaller steps in time and using higher precision operations, but never exactly the same as differential equation (unless the equation is analytically solvable, which for $n>2$ particles it is not).
2) even when using very small steps, powerful computer and high precision, numerical errors are happening with floating point numbers in FPU of your computer because of finite precision of FPU (on x86 the double type uses around 15 significant digits for calculations, anything below that is lost). The compounding effect is to destroy the accuracy of the calculated physical quantities when the step gets too small. That is the nature of FPU and floats: they provide only an approximate results. Using floats, it is not easy to keep the error low or track bounds on the error of the results.
To make your program closer to the physical model, first improve the algorithm. Read about leap frog algorithm and adaptive step algorithms. You can get a decent (eye-satisfactory) behaviour with those.
Then, if you have apparently well behaving simulation, plot energy and angular momentum of your system. They should be constant; if they are not, you know something is wrong. If you are interested in how much error there is in the result of any quantity, you will need to use a more reliable arithmetical system than floating point arithmetics, such as interval arithmetics, or Gosper's continued fraction arithmetics, to keep tabs on the errors of the calculated quantities.
A: Euler method assumes constant force throughout the step, and when the planets are near the force is changing rapidly causes a) incorrect results and b) asymmetry in the loading between approach and separation.
You can visualize this with the following graph. A rapidly changing force (red) is sampled at fixed intervals and the value at the beginning of the interval is used for the area (change in velocity) of the bars.

A: Nowadays, using Euler's method for conservative forces is definitely a bad idea. There is no excuse for that bad idea. Velocity Verlet's method is definitely better (it is time-reversible, symplectic algorithm) and it is as simple to implement as Euler.
In the main loop over the integration times, one updates positions and velocities in the following way:
$$
\begin{eqnarray}
x(t+\delta t) &=& x(t) + v(t)\delta t + \frac{1}{2}a(t)\delta t^2 \\
v(t+\delta t/2)&=&v(t) + a(t)\frac{\delta  t}{2} \\
a(t+\delta t) &=& {\text {force}}(x(t+\delta t))/{\text {mass}} \\
v(t+\delta t)&=&v(t+\delta t/2) + a(t+\delta t)\frac{\delta  t}{2}
\end{eqnarray}
$$
with one evaluation of the force per cycle, exactly like with Euler's method. 
The discretization error in one time step id $O(\delta t^3)$ on position and velocity, which becomes $O(\delta t^2)$ over a finite time (global error). Much better than the $O(\delta t)$ global error of Euler formula.
Moreover, the angular momentum is conserved within machine precision, due to the symplectic nature of the algorithm. Notice that "double precision" i.e. standard 64 bits floating point numbers are more than adequate. Roundoff errors are usually negligible in typical calculations of this kind.
However even symplectic algorithm may have hard time to integrate a strongly eccentric orbit like the precessing orbit of your plot. In that case, one can easily expect large discretization errors in the part of the trajectory where the velocity increases significantly. In that case, variable step algorithms usually perform better and even  Velocity Verlet may be modified to work with a time step inversely proportional to the speed.
A final comment on the choice of integrating the equations of motion for the relative coordinate ${\bf r = r_2-r_1}$, introducing the reduced mass. Numerically it is almost trivial to go from the two-body to the general N-body problem. Therefore it is much ore general to solve the two-body problem directly as a coupled set of two (vector) differential equations, without need of introducing the reduced mass.
