我的C代码有问题。它是2D中N体问题的强制解决方案。有时我将NaN作为struct instance params值。
我的猜测是分裂出了问题。我已经分析了很多案例,但仍然找不到导致值为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);
}
}
}
由于(i)-NaN
的基本参数为负,或者(ii)浮动中的笑话数字累积,sqrt
结果的一个常见(并且到目前为止最可能的解释)生成pow
是负面的论据。点变量:bodyInstance->vx +=
&c。将累积舍入误差。
在调用sqrt
之前检查一下这个案例。
你也会得到一个像NaN
这样的表达式的0.0 / 0.0
,但我从未见过平台在这种情况下产生负面的NaN
。
所以这里的道德是更喜欢dx * dx
和pow(dx, 2)
:前者更准确,不容易受到负面dx
的意外结果,当然也不会慢。更好的是,使用C标准库中的hypot
。
在updateBodyForces
中,您测试两个浮点值是否相等。它们的差异可能与最后一点差不多,约为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。
经过一些测试:ε值难以猜测。由于你可以使用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
,你可以使wibbly-wobbly的时空黑客变得不那么明显。