Hydrogen Ground State Energy Using the non-relativistic Schrodinger equation, the energy levels for hydrogen are found to be $E_{n} = -\frac{1}{2n^2}mc^2\alpha^2$. Using the relativistic Dirac equation, the energy levels for hydrogen are found to be $$E_{n} = -m_{e} c^{2} \left(1 - \left(1 + \left(\frac{\alpha}{\left(n - j - 1/2 + \sqrt{\left(j + 1/2\right)^2 - \alpha^2}\right)}\right)^2\right)^{-\frac{1}{2}}\right)\ .$$ Often, for the non-relativistic energy levels, the mass of the electron is used for $m$. Though the mass used should actually be the reduced mass of a hydrogen atom, $\mu$, since it's a two body problem. Using $m_{e}$, the ground state energy of hydrogen is calculated to be $-13.60569 eV$. Using the reduced mass, the ground state energy is calculated to be $-13.5983 eV$. Since using the reduced mass is supposed to be more accurate than just using the electron's mass, I would expect it to be closer to the relativistic energy values. Yet using $n=1$ and $j=\frac{1}{2}$, the ground state energy is calculated to be $-13.60587 eV$, which is closer to the value of the energy calculated using just the electron's mass, which isn't what I expected. Why is this? Should the reduced mass be used for the Dirac energy levels as well?
 A: The relativistic generalization of the formula; and the replacement of the electron mass by the reduced mass are clearly two basically independent steps (at least in the leading approximation), and both of them have to be applied to agree with the experiment.
The most accurate experimental value of the energy is $-13.59844\,{\rm eV}$, see NIST:

http://webbook.nist.gov/cgi/cbook.cgi?ID=C12385136&Units=CAL&Mask=20

This takes the relativistic corrections into account, as well as the motion of the nucleus (i.e. the replacement of the electron mass by the reduced mass). Your value $-13.60587\,{\rm eV}$ was a theoretical value calculated with the electron mass, not reduced mass (although someone may have used slightly different values of the fine-structure constant etc. than you have, therefore the tiny discrepancy), so this value can't be expected to match the experiment too accurately.
If one goes to higher orders, these two effects (relativistic; proton's motion) cease to be truly independent of each other, and higher-order terms generally exist that can't be incorporated by a simple replacement of $m_e$ by $\mu$ or by the simple relativistic square roots.
The exact calculations at that precision actually require the full quantum electrodynamics, see e.g.

https://arxiv.org/abs/quant-ph/0511197

but only the most precise experiments may be sensitive to these corrections. The binding energy isn't measured too accurately. On the other hand, the transition energies are measured much more accurately due to spectroscopy and they're sufficiently good checks to verify whether one has incorporated all the corrections – and their mutual interactions – correctly.
Note that in the exact treatment, there are additional small terms that have to be taken into account to get "truly precision values" to match the experiment, such as the interaction of the electron's and nucleus' magnetic moment; and the nonzero size of the nucleus that truncates the Coulomb potential for small enough $r$ (distance of the electron and the nucleus).
