-1.#IND0 错误在 C 中使用 Runge-Kutta 方法答案

【问题标题】：-1.#IND0 error using Range-Kutta method in C-1.#IND0 错误在 C 中使用 Runge-Kutta 方法
【发布时间】：2012-10-20 00:52:43
【问题描述】：

我试图使用 4 点 Range-Kutta 方法求解耦合的一阶微分方程。输出m 的值时，我得到-1.#IND0 错误。我知道这可能是 NaN，但对我来说没有意义，因为 m 的值应该会增加，而 -1IND0 在有效值之间。这是我的输出示例：

3110047776596300800000000000000000000.00000 35953700.00
-1.#IND0 35984000.00
-1.#IND0 36013700.00
3721056749337648900000000000000000000.00000 36042800.00
-1.#IND0 36071400.00
4132402773947312100000000000000000000.00000 36099500.00
-1.#IND0 36127200.00
4546861919240663800000000000000000000.00000 36154400.00

这是我的代码：

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

#define pi 3.141592654


double f(double p, double m, double r)
{
    return -0.000000000000000012812899255404507 * m * pow(p, 1.0/3) / (r * r);
}

double g(double p, double r)
{
    return 4 * pi * r * r * p;
}

int main()
{
    double  p_c,        //central density
            p,              //densities
            m,              //masses
            f_val[4],       //arrayed f
            g_val[4],       //arrayed g
            r = 1e-15,      //radius
        dr = 100,       //radius increment
        p_0 = 0.001;    //effective zero density
double p_min = 1e6;
double p_max = 1e14;
int i;                  //Loop counter

FILE *data=fopen("dwarf.txt", "w");//Output file

for(p_c = p_min; p_c <= p_max; p_c += (p_max - p_min) / 100)
{
    p = p_c;
    m = (4.0/3) * pi * r * r * r * p_c;

    while(p > p_0)
    {
        //fprintf(data, "%.5lf %.2lf %.2lf\n", p, m, r);

        f_val[0] = f(p, m, r) * dr;
        g_val[0] = g(p, r) * dr;

        f_val[1] = f(p + f_val[0]/2, m + g_val[0]/2, r + dr/2) * dr;
        g_val[1] = g(p + f_val[0]/2, r + dr/2) * dr;

        f_val[2] = f(p + f_val[1]/2, m + g_val[1]/2, r + dr/2) * dr;
        g_val[2] = g(p + f_val[1]/2, r + dr/2) * dr;

        f_val[3] = f(p + f_val[2], m + g_val[2], r + dr) * dr;
        g_val[3] = g(p + f_val[2], r + dr) * dr;

        m += (g_val[0] + 2 * g_val[1] + 2 * g_val[2] + g_val[3]) / 6; 
        p += (f_val[0] + 2 * f_val[1] + 2 * f_val[2] + f_val[3]) / 6;

        r += dr;
    }

    fprintf(data, "%.5lf %.2lf\n", m, r);
    printf("%.5lf %.2lf\n", m, r);
}
exit;
}

【问题讨论】：

-0.000000000000000012812899255404507 这完全是 16 位有效数字。标准 (double) 文字是否支持该精度的每一点并不明显。在long double 比double 长的平台上使用-0.000000000000000012812899255404507L 会得到更多。（您可能想要也可能不想要，这取决于...）
嗯...而且您只将 pi 定义为 10 位有效数字，因此 f 中的多余部分是徒劳的。考虑const double pi = 4.0 * atan(1.0)。
您能否尝试使用您的编译器使用建议的编辑来编译代码？我似乎无法让它工作。对不起，我是 c 的初学者。
你们有同样的错误吗？
顺便说一句---我的笔记不会影响您遇到的问题，我只是记下一些您可能想要了解的事情。但是你需要习惯调试......如果你坚持编程，你会从现在开始做很多事情。

标签： c physics numerical-methods

【解决方案1】：

我得到了 nan 的。在cygwin上编译运行：

3110047776596300799965078807132504064.00000 35953700.00
nan 35984000.00
nan 36013700.00
3721056749337648263817730951571570688.00000 36042800.00
nan 36071400.00
4132402773947312079489066295688691712.00000 36099500.00
nan 36127200.00
4546861919240663813565041399809703936.00000 36154400.00

自从我研究龙格-库塔已经有一段时间了...查看您的代码，我认为 r 是自变量，dr 是步长，m 是您要解决的因变量。我很困惑 p 是什么。你能给我们更多的细节吗？如果我能看到你试图求解的实际方程，那就更有意义了。

【讨论】：

德里克，在 StackOverflow 上，“答案”应该是真正回答问题的东西。寻求更多信息的回复属于原始问题的 cmets。欢迎使用 StackOverflow。
抱歉，我在该问题的任何地方都没有看到评论选项，但显然我可以对自己的答案发表评论......
不幸的是，您需要一定的声誉 (50) 才能评论其他人的问题或答案，IIRC。所以这样做很好。