Direct N-body simulation of star clusters or other stellar systems has proven itself an extremely powerful tool to study the structure and evolution of stellar systems, since the pioneering work by von Hörner ([1960,1963]). The only known way to do experiments on stellar systems is to construct their models in computers, because we cannot do laboratory experiments on stellar systems,
Of course, N-body simulation is not the only way to construct computer models of star clusters. One could use Monte-Carlo (Henon [1971], for recent development see Giersz [1998]) or direct integration of Fokker-Planck (FP) equations (Cohn [1980], Takahashi [1996], Drukier et al. [1999]). In particular, recent advance in the treatment of the two-dimensional [f(E,J)] FP equation made it possible to use FP code to the study of the evaporation of clusters in the tidal field of the parent galaxy.
The main advantage of these methods over direct N-body simulation is
the calculation cost. While it is still out of reach to perform the
direct simulation of, say, a typical globular cluster with 
,
the FP equation is valid in the limit of 
. Therefore, when the assumption of large-N limit is good,
methods based on the FP approximation can deliver reliable results for
the calculation cost orders of magnitudes smaller than that of direct
N-body calculation. In principle, one can do the N-body simulation
with small N and then try to extrapolate that result to larger
N. However, this problem, which is sometimes called ``the scaling
problem'' has proven itself a complex and difficult
problem(Heggie et al. [1998]). The main reason of the difficulty is
that there are many physical processes of which timescales depend on
the number of particles in different ways, and that there is no way to
adjust all relevant timescales consistently.
On the other hand, methods based on FP approximation do have their own
limitations. For example, it's pretty difficult to extend the
formulation beyond the spherical symmetry. Moreover, they also suffer
the problem of the time scaling, only difference is that they are
affected from the opposite direction. For example, when we want to
include the dynamical effect of the slowly varying potential, FP
approximation goes into trouble. Consider the tidal stripping of
globular clusters with non-circular orbit. The orbital timescale is of
the order of 
years. The orbital timescale of the stars around
the tidal boundary of the cluster is also 
years. On the other
hand, the half-mass relaxation time is of the order of 
years. Therefore, the timestep to integrate the FP equation is at
the largest 
years, and in practice much smaller than that. Clearly,
the assumption that the dynamical timescale is smaller than the
Fokker-Plank timestep is broken, and FP calculations tends to grossly
overestimate the escaper rate.
As beautifully demonstrated by Takahashi and Portegies Zwart ([1998]), one can incorporate the escaping timescale into FP calculation. However, in order to do so the model parameter must be adjusted so that the agreement with reference N-body calculations is achieved. Moreover, this agreement has been so far achieved only with N-body simulations with simple static ``tidal cutoff'', where no tidal field is imposed and stars outside the tidal radius are just removed. Whether or not FP calculations can reproduce the result of N-body calculations with tidal fields is an open question.
The ultimate solution for the scaling problem is to perform N-body simulations with sufficiently large number of particles. In the first half of this paper, we describe our effort in that direction, the GRAPE project to develop and use special-purpose computers for N-body simulations. In the last half of this paper, we'd like to discuss briefly our recent work on N-body simulations of galactic centers with massive black holes, where again the scaling problem shows up.