Brandt Kurowski December 16, 1996 Astrophysics Final Project

best with Netscape 3

N-BODY SIMULATION
direct summation + fast fourier transform

The physics of stellar systems, especially galaxies, may seem easy at first. After all, it's just straight Newtonian mechanics applied to systems of self gravitating particles. However, a typical spiral galaxy contains about 10¹² stars, making any sort of analytic approach beyond impossible. The alternative is to do numerical simulations. With the computer, brute force approaches to solving otherwise impossible physical systems become a reality.

As my final project for this course, I originally set out to simulate the collision of two spiral galaxies. My first attempt was with a direct summation program. My primary concern was the number of particles, n, to be used. I figured that after I wrote the program, I could adjust n to be as high as possible while still allowing for a reasonable running time. Of course, it would obviously be much smaller than the number of stars (~10¹²).

Softened gravitational potential (hi-res)

This meant that the particles that I would be simulating wouldn't really be stars, but instead could be thought of as clusters of stars (of equal mass). As such, I couldn't really consider them to be point masses, so I used the standard practice of substituting a softened potential for Newton's original one. The equation for gravitational force became:

where e is the "softening parameter," and can be thought of as the radial size of the clusters. This kept the force felt between two particles from becoming arbitrarily large when they were near one another, so the size of the discrete time step wouldn't be so critical. The size of the timestep, dt was the other important parameter. I measured it in units of ¹/_T, where T represents the period of rotation for a star located halfway from the center of the disk. Thus dt equals the fraction of a rotation that the disk goes through in each time step. It seemed to be as logical a choice as any. I would keep dt at the largest possible value that gave results similar to those with much smaller dt.

I soon discovered the importance of the order of complexity for algorithms, and became familiar with big-O notation. Direct summation runs in O(n²) time, which meant that as I added particles to make the simulation more realistic, it quickly became an impossible problem. In two-dimensions with n set at only 200, and dt as high as ¹/₁₀, one galactic rotation took 10 seconds. Of course, these numbers are ridiculous for any realistic simulation. Only slightly less useless would be n=500 and dt=¹/₅₀. This took 5 minutes per revolution, though. And with n=1000 and dt=¹/₁₀₀, one revolution required 45 minutes of calculation. These times are for a 75 MHz Pentium, which may not be the epitome of computational number crunching power, but is at least half as fast as any other computer on campus. It became clear to me that with a whole weekend of time on a fast computer, I wouldn't be able to do a simulation with n much more than 1000, even with a pretty coarse time step. This wasn't going to work for simulating the collision of a pair of galaxies, so I investigated other algorithms.

The next obvious choice was to compute the potential at every point in space on some sort of Cartesian grid, and use the gradient of the potential to find the total force acting on a given particle. Ordinarily this would take even longer than direct summation, unless one uses a fast fourier transform (FFT) to perform a convolution of a generic potential at the position of each mass. The alogrithm I used computed the mass inside each cell of the grid, transformed it, multiplied it by an already transformed, pre-computed gravitational potential, then transformed the resulting potential back. What mattered most this time was the resolution, N of the three dimensional grid that I used. The FFT runs in O(4N³ log₂(2N³)), which seemed much more promising than direct summation. This time, the number of particles n really didn't matter, because it just changed a constant in front of the order of complexity. I have summed up the running times in the following tables:

N=8 (512 cells) dt=0.01 n particles time (seconds) 500 20 1000 25 2000 30

N=16 (4096 cells) dt=0.01 n particles time (seconds) 500 149 1000 152 2000 160

Screen image, after 1 revolution

At first, I though that these times were great. I was limited to N=16 because of DOS memory constraints, but I planned on running the finished version on akbar and moving up to N=32 (for a total of 32768 cells). However, I noticed strange things happening in the disk. Four "spiral arms" formed early in the simulation. At first I thought "Wow, I've got spiral arms forming," but then realized that there shouldn't be any spiral arms in my simulation. A bar instability was a possibility, but I had included nothing to generate spirals. Not to mention the near impossibility of four armed spirals. After trying to remove any sort of periodicity that could account for a "crystal lattice" of galaxies being simulated, I realized that the discrete potential I was generating had an "X" of zero force. Because the potential was symmetric about the x and y axes, there were troughs that particles would settle into, and because of the discrete nature of the grid the trough was wide and flat. The only way to try to compensate for this was to, for each particle, do a direct sum of the gravitational force generated by the other particles in its grid cell.

Discrete softened potential (hi-res)

After all, a particle's neighbors weren't affecting it at all because of the coarseness of the grid. However, adding a direct sum over local neighbors in dramatically increased the running time. The effect was inversely related to the grid size, because as N increased, each cell was smaller and therefore each particle has fewer neighbors. The increase in time was still drastic at N=16, though, as is summed up in the tables below:

N=8 (512 cells) dt=0.01 n particles time (seconds) 500 303 1000 1232

N=16 (4096 cells) dt=0.01 n particles time (seconds) 500 381 1000 1118

I finally decided that this was the wrong type of project to undertake with the computing resources available at Marlboro at this time. Instead, it would've been better to simply use potentials that resembled those of spiral galaxies, and use massless test particles. Then I could observe the effects of the colliding potentials of two galaxies on the stars within the galaxies, assuming that the overall structure of each wasn't severely altered by the interaction. I do, however, think that this project was a valuable exercise in numerical programming, and it may become viable when the college acquires more computing power.

SOURCE CODE

1. NBODY.CPP Two-dimensional direct summation
2. REAL.CPP Three-dimensional FFT with optional local direct sums
   Each of these can be run in real-time graphics mode, or in text mode with a timer displayed. Some system calls may need to be altered when compiling this code on operating systems other than DOS.

REFERENCES

1. Binney, J. and Tremaine, S., Galactic Dynamics, Princeton University Press, 1987
2. Cormen, Leiserson, and Rivest, Introduction to Algorithms, The MIT Press, 1990
3. James, J.F., A student's guide to Fourier Transforms, Cambridge University Press, 1995
4. Shu, Frank, The Physical Universe, University Science Books, 1982

RELATED LINKS

• Grand Challenge Cosmology Consortium
• UW HPCC/ESS N-Body Shop, especially Galaxy Harassment
• Large Scale Structure and Galaxy Formation