c - -1。在C中使用Runge-Kutta方法出现＃IND0错误

我试图使用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;
}

最佳答案

我有南的。编译并在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

自从我学习Runge-Kutta以来已经有一段时间了……看一下代码，我认为r是自变量，dr是步长，m是您要解决的因变量。我很困惑p是什么。您能给我们更多细节吗？如果我能看到您要求解的实际方程式，那就更有意义了。