重力计算导致 NaN。没有明确的理由

Question

我的 C 代码有问题。 它是二维 N 体问题的强制解决方案。 有时我将 NaN 作为结构实例参数值。

我的猜测是除法有问题。我分析了很多案例，但仍然找不到导致值为 NaN 的发生模式。

这是我的代码：

#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

在updateBodyForces 中测试两个浮点值是否相等。它们的差异可能只有最后一个 bit，大约 1/10,000,000。

在这之后，你取它们的差平方的平方根，所以结果可能是0（真的是零，0.0000000...），这不是问题，但是那么你除这个数字。那是 NaN 的来源。

替换这部分

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

基于FLT_EPSILON 进行更明确的测试。有关详细说明，请参阅Floating point equality and tolerances。

经过一些测试：epsilon 值很难猜到。既然您可以在极端情况下使用 dist = 1，请将其添加到测试下方，force 行上方，以确保：

if (dist < 1)
    dist = 1;

所以你肯定不会得到任何 NaN。这导致了这个更简单的功能：

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;
}

您也可以通过将1 替换为更小的值来使摇摆不定的时空破解变得不那么明显。

Answer 2

-NaN 结果的常见（也是迄今为止最可能的解释）产生是sqrt 的否定参数，原因是 (i) @987654324 的 base 参数@ 为负数，或 (ii) 在浮点变量中累积笑话数字：bodyInstance->vx += &c。会累积舍入误差。

在调用sqrt之前检查该情况。

您还将获得带有类似0.0 / 0.0 表达式的NaN，但我从未见过平台在这种情况下产生否定的NaN。

所以这里的道德是更喜欢dx * dx 而不是pow(dx, 2)：前者更准确，不易受到带有负面dx 的意外结果的影响，当然也不会更慢。更好的是，使用 C 标准库中的 hypot。

参考：http://en.cppreference.com/w/c/numeric/math/hypot