Gravity computations resulting in NaN. No clear reason

Question

I have problem with my C code. It's force solution on N-bodies problem in 2D. Sometimes I get NaN as struct instance params value.

My guess is that something is wrong with division. I have analyzed a lot of cases, but still could not find an occuring pattern that results in values being NaN.

Here is my code:

#include <stdio.h>
#include <math.h>
#include <time.h>
#include <stdlib.h>

double G;  
int N;
int T;
int COUNT;

typedef struct 
{
    double rx, ry;
    double vx, vy;
    double fx, fy;
    double mass;
} Body;

void updateBody (Body* bodyInstance, int timestamp) {
    bodyInstance->vx += timestamp * bodyInstance->fx / bodyInstance->mass;
    bodyInstance->vy += timestamp * bodyInstance->fy / bodyInstance->mass;
    bodyInstance->rx += timestamp * bodyInstance->vx;
    bodyInstance->ry += timestamp * bodyInstance->vy;
};

void updateBodyForces (Body* bodyA, Body* bodyB) {
    double dx, dy, dist, force;
    dx = bodyB->rx - bodyA->rx;
    dy = bodyB->rx - bodyA->rx;
    // collision/same place in spacetime hack
    if (bodyB->rx == bodyA->rx && bodyB->ry == bodyA->ry) {
        dist = 1;
    } else {
        dist = sqrt(pow(dx, 2) + pow(dy, 2));
    }
    force = (G * bodyA->mass * bodyB->mass) / (pow(dist, 2) + 100);
    bodyA->fx += force * dx / dist;
    bodyA->fy += force * dy / dist;
}

void resetBodyForces (Body* bodyInstance) {
    bodyInstance->fx = 0;
    bodyInstance->fy = 0;
}

void getRandomBody (Body* bI) {
    bI->rx = rand() % 10;
    bI->ry = rand() % 10;
    bI->vx = rand() % 10;
    bI->vy = rand() % 10;
    bI->fx = 0;
    bI->fy = 0;
    bI->mass = 20;
}

int main( int argc, char *argv[] )  {
    G = argc >= 2 ? atof(argv[1]) : 0.01;  
    N = argc >= 3 ? atoi(argv[2]) : 3;
    T = argc >= 4 ? atoi(argv[3]) : 1;
    COUNT = argc >= 5 ? atoi(argv[4]) : 10;
    srand(time(NULL));
    Body bodies[N];
    for (int i=0; i<N; i++) {
        getRandomBody(&bodies[i]);
    }
    for (int i = 0; i < COUNT; i++) {
        for (int j = 0; j < N; j++) {
            resetBodyForces(&bodies[j]);
            for (int k = 0; k < N; k++) {
                if (j != k) {
                    updateBodyForces(&bodies[j], &bodies[k]);
                }
            }
        }
        for (int j = 0; j < N; j++) {
            updateBody(&bodies[j], T);
        }
    }
}

Answer 1

In updateBodyForces you test two floating point values for equality. They may differ by as little as the very last bit, about 1/10,000,000.

Right after this you take the square root of their difference squared, and so the result may be 0 (really really zero, 0.0000000...), which is not a problem, but then you divide by that number. That is the source of the NaN.

Replace this part

// collision/same place in spacetime hack
if (bodyB->rx == bodyA->rx && bodyB->ry == bodyA->ry) {
    dist = 1;
}

with a more explicit test based on FLT_EPSILON. See Floating point equality and tolerances for a longer explanation.

After some testing: the epsilon value is difficult to guess. Since you are okay with a dist = 1 for corner cases, add this below the test, above the force line, to be sure:

if (dist < 1)
    dist = 1;

so you won't get any NaNs for sure. That leads to this simpler function:

void updateBodyForces (Body* bodyA, Body* bodyB) {
    double dx, dy, dist, force;
    dx = bodyB->rx - bodyA->rx;
    dy = bodyB->ry - bodyA->ry;
    dist = sqrt(dx*dx + dy*dy);
    // collision/same place in spacetime hack
    if (dist < 1)
        dist = 1;
    force = (G * bodyA->mass * bodyB->mass) / (pow(dist, 2) + 100);
    bodyA->fx += force * dx / dist;
    bodyA->fy += force * dy / dist;
}

You can make the wibbly-wobbly space-time hack a bit less obvious by replacing the 1 with a smaller value as well.

Answer 2

A common (and by far the most likely explanation here) production of the -NaN result is a negative argument to sqrt, due to either (i) the base parameter to pow being negative, or (ii) an accumulation of joke digits in your floating point variables: bodyInstance->vx += &c. will accumulate rounding errors.

Check that case prior to calling sqrt.

You will also get a NaN with an expression like 0.0 / 0.0, but I've never seen a platform yield a negative NaN in that instance.

So the moral here is to prefer dx * dx to pow(dx, 2): the former is more accurate, not vulnerable to unexpected results with negative dx, and certainly not slower. Better still, use hypot from the C standard library.

Reference: http://en.cppreference.com/w/c/numeric/math/hypot