为什么我的程序输出“nan”而不是双精度数（速度和位置？）

Question

我希望这个程序输出速度和位置的值，但它只输出NaN。我检查了我所有的功能，我没有在任何地方除以零，还能是什么？位置和速度应该模拟一个轨道。

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

double kx1 (double vx);
double ky1 (double vy);

double kvx1(double G, double M, double x, double y);
double kvy1(double G, double M, double x, double y);

double kx2(double G, double M, double h, double x, double y, double vx);
double ky2(double G, double M, double h, double x, double y, double vy);

double kvx2(double G, double M, double h, double x, double y, double vx);
double kvy2(double G, double M, double h, double x, double y, double vy);

double kx3(double G, double M, double h, double x, double y, double vx);
double ky3(double G, double M, double h, double x, double y, double vy);

double kvx3(double G, double M, double h, double x, double y, double vx);
double kvy3(double G, double M, double h, double x, double y, double vy);

double kx4(double G, double M, double h, double x, double y, double vx);
double ky4(double G, double M, double h, double x, double y, double vy);

double kvx4(double G, double M, double h, double x, double y, double vx);
double kvy4(double G, double M, double h, double x, double y, double vy);

int main()
{

    double M,G,y0,vx,vy,x,y,h,t, positionx,positiony,velocityx,velocityy;
    G=6.67*pow(10,-11);
    h=10;
    t=0;
    M=5.97*pow(10,24);
    y=0;
    //y0=0;
    x=1738*pow(10,3);
    vx=0;
    vy=1.023*pow(10,3);

    FILE *fp, *fopen();
    fp=fopen("Orbit", "w");
    for (t=0; t<=1000; t++)
    {

        printf("%lf,%lf,%lf,%lf,%lf\n", y,x,vx,vy,t);
        fprintf(fp,"%lf,%lf,%lf,%lf,%lf\n", y,x,vx,vy,t);

        positionx=x+(h/(double)6)*(kx1(vx)+2*kx2(G,M,h,x,y,vx)+2*kx3(G,M,h,x,y,vx)+kx4(G,M,h,x,y,vx));
        positiony=y+(h/(double)6)*(ky1(vy)+2*ky2(G,M,h,x,y,vy)+2*ky3(G,M,h,x,y,vy)+ky4(G,M,h,x,y,vy));
        velocityx=vx+(h/(double)6)*(kvx1(G,M,x,y)+2*kvx2(G,M,h,x,y,vx)+2*kvx3(G,M,h,x,y,vx)+kvx4(G,M,h,x,y,vx));
        velocityy=vy+(h/(double)6)*(kvy1(G,M,x,y)+2*kvy2(G,M,h,x,y,vy)+2*kvy3(G,M,h,x,y,vy)+kvy4(G,M,h,x,y,vy));

            /*if (positiony==y0)
            {
                break;
            }*/

        x=positionx;
        y=positiony;
        vx=velocityx;
        vy=velocityy;

    }

    fclose(fp);
    return 0;

}

double kx1 (double vx)
{
double kx1ans;
kx1ans=vx;
return kx1ans;
}

double ky1 (double vy)
{
double ky1ans;
ky1ans=vy;
return ky1ans;
}

double kvx1(double G, double M, double x, double y)
{
double kvx1ans;
kvx1ans=(-1*G*M*x)/(pow((pow(x,2)+pow(y,2)),1.5));
return kvx1ans;
}

double kvy1(double G, double M, double x, double y)
{
double kvy1ans;
kvy1ans=(-1*G*M*y)/(pow((pow(y,2)+pow(y,2)),1.5));
return kvy1ans;
}

double kx2(double G, double M, double h, double x, double y, double vx)
{
double kx2ans;
kx2ans=vx+(h*kvx1(G,M,x,y))/2;
return kx2ans;
}

double ky2(double G, double M, double h, double x, double y, double vy)
{
double ky2ans;
ky2ans=vy+(h*kvy1(G,M,x,y))/2;
return ky2ans;
}

double kvx2(double G, double M, double h, double x, double y, double vx)
{
double kvx2ans;
kvx2ans=(-1*G*M*(x+(h*kx1(vx))/2))/(pow((pow((x+(h*kx1(vx))),2)+pow(y,2)),1.5));
return kvx2ans;
}

double kvy2(double G, double M, double h, double x, double y, double vy)
{
double kvy2ans;
kvy2ans=(-G*M*(y+(h*ky1(vy))/2))/(pow((pow((y+(h*ky1(vy))),2)+pow(x,2)),1.5));
return kvy2ans;
}

double kx3(double G, double M, double h, double x, double y, double vx)
{
double kx3ans;
kx3ans=vx+(h*kvx2(G,M,h,x,y,vx))/2;
return kx3ans;
}

double ky3(double G, double M, double h, double x, double y, double vy)
{
double ky3ans;
ky3ans=vy+(h*kvy2(G,M,h,x,y,vy))/2;
return ky3ans;
}

double kvx3(double G, double M, double h, double x, double y, double vx)
{
double kvx3ans;
kvx3ans=(-G*M*(x+(h*kx2(G,M,h,x,y,vx)/2)))/(pow((pow((x+(h*kx2(G,M,h,x,y,vx))),2)+pow(y,2)),1.5));
return kvx3ans;
}

double kvy3(double G, double M, double x, double y, double h, double vy)
{
double kvy3ans;
kvy3ans=(-G*M*(y+(h*ky2(G,M,h,x,y,vy))/2))/(pow((pow((y+(h*ky2(G,M,h,x,y,vy))),2)+pow(x,2)),1.5));
return kvy3ans;
}

double kx4(double G, double M, double h, double x, double y, double vx)
{
double kx4ans;
kx4ans=vx+(h*kvx3(G,M,h,x,y,vx));
return kx4ans;
}

double ky4(double G, double M, double h, double x, double y, double vy)
{
double ky4ans;
ky4ans=vy+(h*kvy3(G,M,h,x,y,vy));
return ky4ans;
}

double kvx4(double G, double M, double h, double x, double y, double vx)
{
double kvx4ans;
kvx4ans=(-G*M*(x+(h*kx3(G,M,h,x,y,vx))/2))/(pow((pow((x+(h*kx3(G,M,h,x,y,vx))),2)+pow(y,2)),1.5));
return kvx4ans;
}

double kvy4(double G, double M, double h, double x, double y, double vy)
{
double kvy4ans;
kvy4ans=(-G*M*(y+(h*ky3(G,M,h,x,y,vy))/2))/(pow((pow((y+(h*ky3(G,M,h,x,y,vy))),2)+pow(x,2)),1.5));
return kvy4ans;
}

Answer 1

在您的代码y=0 和vx=0 中。您可能会将 0 除以 0，而 0 定义为 NaN。你需要看看函数，比如kvx2(G,M,h,x,y,vx)。

更多info on what can cause Non a Number

编辑：

你在打电话

double kvy1(double G, double M, double x, double y)
{
  double kvy1ans;
  kvy1ans=(-1*G*M*y)/(pow((pow(y,2)+pow(y,2)),1.5));
  return kvy1ans;
}

其中 y 为 0。(pow((pow(y,2)+pow(y,2)),1.5)) = 0，显然是 G*M*y= 0。这导致了上述情况。我没有检查进一步的错误，因为这个就足够了。

只要你使用这种神秘而令人困惑的变量，你就会遇到这种错误。

Answer 2

我把你的任务改成了positiony 到

    double tmp1 = y;
    double tmp2 = h/(double)6;
    double tmp3tmp1 = ky1(vy);
    double tmp3tmp2 = 2*ky2(G,M,h,x,y,vy);
    double tmp3tmp3 = 2*ky3(G,M,h,x,y,vy);
    double tmp3tmp4 = ky4(G,M,h,x,y,vy);
    double tmp3 = tmp3tmp1 + tmp3tmp2 + tmp3tmp3 + tmp3tmp4;
    positiony = tmp1 + tmp2 * tmp3;
    if (t < 2) {
        printf("tmp1 is %f\n", tmp1);
        printf("tmp2 is %f\n", tmp2);
        printf("tmp3tmp1 is %f\n", tmp3tmp1);
        printf("tmp3tmp2 is %f\n", tmp3tmp2);
        printf("tmp3tmp3 is %f\n", tmp3tmp3);
        printf("tmp3tmp4 is %f\n", tmp3tmp4);
        printf("tmp3 is %f\n", tmp3);
        printf("positiony is %f\n", positiony);
    }

并且，当我运行程序时，我得到了

tmp1 为 0.000000 tmp2 是 1.666667 tmp3tmp1 是 1023.000000 tmp3tmp2 是 nan tmp3tmp3 是 2042.120516 tmp3tmp4 是 -295.258623 tmp3 是南 定位是南 南,1731387.732835,-1324.550838,南,1.000000

所以我得出结论，ky2() 中的某些内容是错误的......