Cosmology is entering an era of percent level precision due to current large observational surveys. This precision in observation is now demanding more accuracy from numerical methods and cosmological simulations. In this paper, we study the accuracy
of $N$-body numerical simulations and their dependence on changes in the initial conditions and in the simulation algorithms. For this purpose, we use a series of cosmological $N$-body simulations with varying initial conditions. We test the influence of the initial conditions, namely the pre-initial configuration (preIC), the order of the Lagrangian perturbation theory (LPT), and the initial redshift, on the statistics associated with the large scale structures of the universe such as the halo mass function, the density power spectrum, and the maximal extent of the large scale structures. We find that glass or grid pre-initial conditions give similar results at $zlesssim 2$. However, the initial excess of power in the glass initial conditions yields a subtle difference in the power spectra and the mass function at high redshifts. The LPT order used to generate the ICs of the simulations is found to play a crucial role. First-order LPT (1LPT) simulations underestimate the number of massive haloes with respect to second-order (2LPT) ones, typically by 2% at $10^{14} h^{-1} M_odot$ for an initial redshift of 23, and the small-scale power with an underestimation of 6% near the Nyquist frequency for $z_mathrm{ini} = 23$. Moreover, at higher redshifts, the high-mass end of the mass function is significantly underestimated in 1LPT simulations. On the other hand, when the LPT order is fixed, the starting redshift has a systematic impact on the low-mass end of the halo mass function.
