|The Challenge||The Approach||Publications||Investigators||Research Sponsors|
|Trajectories of 10 bodies interacting through a gravitational force over 5000 time steps. Initial positions and velocities follow a 3D Plummer distribution.|
A central problem in many physical simulations is the evolution of a system of many elements or bodies interacting over large distances. This problem arises in cosmological simulations, where the long-range force is gravity, in molecular dynamics simulations, where the long range force is primarily electrostatics, and in many other areas including vortex formation in computational fluid dynamics and radiosity computation in computer graphics.
The simplest approach to determine the total force experienced by each body due to all other bodies in the system requires O(N2) work per integration timestep and quickly becomes intractable for large systems that must be evolved over many timesteps. The celebrated fast multipole methods reduce the cost to O(N log N) or even O(N) by the hierarchical approximation of the force due to a group of bodies at a distance. These methods have rigorous error bounds but require a spatial decomposition of the bodies to be created and maintained. Furthermore, the efficient implementation of the fast multipole methods for 3D long-range forces is extraordinarily complex, particularly when the bodies are non-uniformly distributed. In general, several thousands of bodies are required before the methods become competitive with alternatives with larger asymptotic complexity.
In our case, we are particularly interested in molecular dynamics simulations, where N may be relatively small (on the order of a few thousand atoms), the number of integration time steps is very large (on the order of 106 or more), and the accuracy required is very high. For this class of problems in particular, as well as N-body problems in general, we are investigating a new class of algorithms that are adaptive in time, and have work complexity that grows as O(N4/3) or even O(N) per integration step when amortized over (N) steps.
Instead of using a fixed time step for all interactions (which is too conservative for most interactions), an approach that equalizes impulse is followed. Mathematically, the constant time step t is traded for constant impulse I defined as Fij ×Tij, where Fij is the force between particles i and j, and Tij becomes the time step used to re-evaluate Fij. Using a constant impulse instead of a constant time step leads to an amortized execution complexity of O(N4/3) per simulation step. Algorithmic improvements that rely on the first and second derivatives of force further reduce the amortized work per step to O(N log N) and O(N), respectively.
|Distribution of the 45 pairwise distances between bodies after 100 time steps in the system shown previously.|
|Schedule for the next force evaluation according to the equal impulse method. Over half of all interactions need not be re-evaluated within the next 10 timesteps, and 20% of the interactions need not be re-evaluated within the next 100 timesteps.|