  The effects of discreteness arising from the use of the N-body method on the
accuracy of simulations of cosmological structure formation are not currently
well understood. After a discussion of how the relevant discretisation
parameters introduced should be extrapolated to recover the Vlasov-Poisson
limit, we study numerically, and with analytical methods we have developed
recently, the central issue of how finite particle density affects the
precision of results. In particular we focus on the power spectrum at
wavenumbers around and above the Nyquist wavenumber, in simulations in which
the force resolution is taken smaller than the initial interparticle spacing.
Using simulations of identical theoretical initial conditions sampled on four
different "pre-initial" configurations (three different Bravais lattices, and a
glass) we obtain a {\it lower bound} on the real discreteness error. With the
guidance of our analytical results, we establish with confidence that the
measured dispersion is not contaminated either by finite box size effects or by
subtle numerical effects. Our results show notably that, at wavenumbers {\it
below} the Nyquist wavenumber, the dispersion increases monotonically in time
throughout the simulation, while the same is true above the Nyquist wavenumber
once non-linearity sets in. For normalizations typical of cosmological
simulations, we find lower bounds on errors at the Nyquist wavenumber of order
of a percent, and larger above this scale. The only way this error may be
reduced below these levels at these scales, and indeed convergence to the
physical limit firmly established, is by extrapolation, at fixed values of the
other relevant parameters, to the regime in which the mean comoving
interparticle distance becomes less than the force smoothing scale.
