0

4点Range-Kutta法を使用して、結合された一次微分方程式を解こうとしていました。の値を出力するとm、-1.#IND0エラーが発生します。mこれが NaN になる可能性があることは承知していますが、 の値が増加する必要があり-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 に答える