0

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

1 回答 1

1

我得到了南的。在 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 是什么。你能给我们更多的细节吗？如果我能看到您要解决的实际方程式，那将更有意义。

于 2012-10-20T03:43:17.193 回答