  We present a new method for generating initial conditions for numerical
cosmological simulations in which massive neutrinos are treated as an extra set
of N-body (collisionless) particles. It allows us to accurately follow the
density field for both Cold Dark Matter (CDM) and neutrinos at both high and
low redshifts. At high redshifts, the new method is able to reduce the shot
noise in the neutrino power spectrum by a factor of more than $10^7$ compared
to previous methods, where the power spectrum was dominated by shot noise at
all scales. We find that our new approach also helps to reduce the noise on the
total matter power spectrum on large scales, whereas on small scales the
results agree with previous simulations. Our new method also allows for a
systematic study of clustering of the low velocity tail of the distribution
function of neutrinos. This method also allows for the study of the evolution
of the overall velocity distribution as a function of the environment
determined by the CDM field.
