Accurately calculating gravity for large systems

Question

I'm in the middle of a program to simulate gravity between particles for > 50000 particles. It works so far, but the way I calculate gravity seems to be a tad askew. Generally, it works. It fails when particles get too close (distance between them approaching 0, gravitational force approaching infinity). I believe this is because of how the implementation chunks up time. It can only calculate the gravity between frames, and I can shorten the time between frames but the problem still exists.

Here's the implementation:

//The time between frames
static final double planc = 0.002;

//Calculates distance^2 between two particles
public double d2(Particle p) {

    return (x - p.x) * (x - p.x) + (y - p.y) * (y - p.y);

}

//Gets Vector from on particle pointing to the next
public Vector getVec(Particle p) {

    return new Vector(p.x - x, p.y - y);

}

//The function that actually (independently) moves the particles according 
//to the forces acting on them
public void move() {

    //v is a Vector that represents the force (not acceleration) 
    //acting on the particle
    x += v.dx * planc / m;
    y += v.dy * planc / m;

}

//Actually calculates gravity
public boolean grav(Particle p) {

    double d2 = d2(p);

    //Happens once because of how the function is used (unimportant)
    if (d2 == 0)
        return false;

    //-1000 is the G value
    double f = -1000 * m * p.m / d2;

    //Creates acceleration-vectors
    Vector v = getVec(p);
    v.normTo(f * planc, Math.sqrt(d2));//sets v's magnitude to f * planc. 2nd parameter saves a calculation
    //Basically formats a vector which represents gravitational acceleration

    //"adds" force-vector to both particles
    this.v.sub(v);
    p.v.add(v);

    return false;

}

This is fast but sacrifices accuracy. Is there a better way to calculate gravity for large n-body systems (that doesn't screw over performance)?

Answer 1

Gravity for small distances

I've seen two simple ways that are used to fix this problem:

1. Add a small amount to each distance:

   double d2 = d2(p);
   
   d2 += EPSILON_DIST;
   

2. Use a minimum distance:

   double d2 = d2(p);
   
   if (d2 < MIN_DIST)
       d2 = MIN_DIST;
   

Either way will prevent small distances from causing huge accelerations. Since d2 is in units of distance squared, your EPSILON_DIST and MIN_DIST should also be squared distances.

More efficient gravity computation

Your gravity calculation could be improved a little. I answered another almost identical question here, so you can take a look at the solution I came up with there. The main point is that you can reduce the number of divisions to only one in grav() and zero in move(). Although it doesn't exactly match your Particle class, I think you can extrapolate from it. Here is the code, pasted from that other question:

public void interact(Body other) {
    double dx = other.getX() - x;
    double dy = other.getY() - y;
    double r  = calculateDistance(dx, dy);
    double inv_r3 = 1.0 / (r * r * r);

    /* Precalculate force component (1/r^2) times direction (dx/r) = dx / r^3 */
    dx *= inv_r3;
    dy *= inv_r3;

    /* calculate accelerations for both bodies */
    vx += other.getMass() * dx;
    vy += other.getMass() * dy;
    other.addToVx(mass * -dx);
    other.addToVy(mass * -dy);
}

public void update() {
    x += vx;
    y += vy;
}