Miller  was the first to show that small changes in the initial conditions of an N-body system result in exponentially diverging solutions. Lecar  co-ordinated a study between many researchers, each of whom independently computed the solution to an N-body problem with identical initial conditions. They found that different algorithms and computers gave results in which some measures differed by as much as 100%. More recent work on the growth of errors includes Kandrup and Smith , who showed that under a large range of parameters, the time scale over which small perturbations grow by a factor of e (Euler's constant, ), called the e-folding time, is comparable to the crossing time (the average time it takes a particle to cross the system once). Goodman, Heggie, and Hut  developed a detailed theory of the growth of small perturbations, and verified it with simulation to show that the exponential growth of small errors results mostly from close encounters, which occur infrequently. This is an interesting result because it says that, even though the full phase-space solutions may experience exponential error growth, this growth is due mainly to a few particles that undergo close encounters. Perhaps the results of collisionless systems can be trusted for longer than collisional ones, since close encounters have a much smaller effect in the former. Kandrup and Smith  showed that as softening is increased (i.e., the collisionality is decreased), the e-folding time grows slightly faster than linearly. (i.e., the Lyapunov exponent decreases, so errors are magnified more slowly.) They agree with Goodman et al. that the error magnification is due more to the rare individual particles whose errors grow much more quickly than the average, although they claim the global potential also plays a role. Presumably the particles whose errors grow more quickly than average are ones that have suffered close encounters.
Since the time-scale for the growth of errors is so short (the errors can be magnified by a factor of each crossing time), the results of all N-body simulations may be suspect. If the relative error per crossing time for a simulation is , then after about p crossing times, a particle's position will have an error comparable to the size of the system -- in other words, all information will have been lost about the particle's position. Typical simulations today have a p between 4 and 8.
I am not aware of any convincing justification for the belief that statistical measures taken from these simulations are valid, although I tend to share the same ``warm fuzzy'' feelings that most astronomers have -- namely that, in some sense, large N-body simulations give valid statistical results. The problem is, in what sense are they valid, if any; how valid are they; and finally, how badly are we allowed to integrate before validity is lost?