There are many different numerical techniques to follow the evolution of a system of many particles. For earlier reviews see Hockney & Eastwood (1981) and Sellwood (1987). The most frequently used methods for cosmological applications fall in three classes: Particle Mesh (PM) codes, Particle-Particle/Particle-Mesh (P 3M) codes, and TREE codes. All methods have their advantages and disadvantages.
PM codes use a mesh for the density and potential. As a result, their resolution is limited by the size of the mesh. The largest simulations done by the author have been on a 8003 mesh with 3 x 2563 = 1.5 x 108 particles. The SP2 parallel supercomputer at Cornell has been used to run simulations with a 16003 = 4.096 x 109 mesh (Gross, 1997). There are two advantages of the method: i) it is fast (the smallest number of operations per particle per time step of all the other methods), ii) it typically uses very large number of particles. The latter can be crucial for some applications. There are a few variants of PM codes. A ``plain-vanilla'' PM was described by Hockney & Eastwood (1981), and this includes a Cloud-In-Cell density assignment and a 7-point discrete analog of the laplacian operator. Higher order approximations improve the accuracy on large distances, but degrade the resolution (e.g. Gelb (1992)). In an effort to reduce the order of approximation and to increase the resolution, Melott (1986) introduced the staggered mesh. It gives a better resolution on cell-size distances, but particles get self-forces (an isolated particle experiences a force from itself), which might be not a welcome feature.
P 3M codes are described in detail in Hockney & Eastwood (1981) and in Efstathiou et al.(1985). They have two parts: a PM part, which takes care of large-scale forces, and a PP part, which adds a small-scale particle-particle contribution. The simulations usually have 643- 1003 particles. Because of strong clustering at late stages of evolution, the PP part becomes prohibitively expensive once large objects start to form in large numbers. Significant speed is achieved in a modified version of the code which introduces subgrids (next levels of PM) in areas with high density (AP 3M code of Couchman (1991)). With modification this code runs as fast as TREE code even for heavily clustered configurations (Couchman 1991).
TREE codes are the most flexible code in the sense of the choice of boundary conditions (Appel 1985, Barnes & Hut 1986, Hernquist 1987). They are also more expensive than PM: it takes 10-50 times more operations. Bouchet & Hernquist (1986) and Hernquist, Bouchet & Suto (1991) have extended the code for the periodical boundary conditions, which are important for simulating large-scale fluctuations.
Multigrid methods were introduced long ago, but only recently have they started to show a potential to produce real results (Anninos, Norman & Clarke 1994, Suisalu & Saar 1995, Kravtsov et al.1997). At present the most advanced and fastest multigrid code has been developed by Kravtsov et al.(1997).