In the next decade, cosmological surveys will have the statistical power to detect the absolute neutrino mass scale. $N$-body simulations of large-scale structure formation play a central role in interpreting data from such surveys. Yet these simulations are Newtonian in nature. We provide a quantitative study of the limitations to treating neutrinos, implemented as $N$-body particles, in $N$-body codes, focusing on the error introduced by neglecting special relativistic effects. Special relativistic effects are potentially important due to the large thermal velocities of neutrino particles in the simulation box. We derive a self-consistent theory of linear perturbations in Newtonian and nonrelativistic neutrinos and use this to demonstrate that $N$-body simulations overestimate the neutrino free-streaming scale, and cause errors in the matter power spectrum that depend on the initial redshift of the simulations. For ${z}_{i}\ensuremath{\lesssim}100$, and neutrino masses within the currently allowed range, this error is $\ensuremath{\lesssim}0.5%$, though represents an up to $\ensuremath{\sim}10%$ correction to the shape of the neutrino-induced suppression to the cold dark matter power spectrum. We argue that the simulations accurately model nonlinear clustering of neutrinos so that the error is confined to linear scales.