-1。＃IND0錯誤使用範圍庫塔方法在C

我試圖解決一個耦合一階微分方程使用4點距離 - 庫塔方法。當輸出m的值時，出現-1.#IND0錯誤。我知道這可能是NaN，但它對我來說沒有意義，因爲m的值應該增加，並且我得到-1IND0之間的有效值。這裏是我的輸出樣本：-1。＃IND0錯誤使用範圍庫塔方法在C

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; 
}

來源

2012-10-20 sk1ll3r

'-0.000000000000000012812899255404507'這完全是16位有效數字。這個標準（'double'）文字將支持這種精度的每一點並不明顯。使用'-0.000000000000000012812899255404507L'會在'long double'長度超過'double'的平臺上獲得更多信息。（你可能會也可能不想要，這取決於...） – dmckee

嗯......而你只定義了pi到10個有效數字，所以額外的f是零。考慮'const double pi = 4.0 * atan（1.0）'。 – dmckee

考慮使用long double。 – askmish

我男的。編譯並在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爲自變量，博士是步長，m是你正試圖解決的因變量。我很困惑p是什麼。你能給我們更多的細節嗎？如果我能看到你正在試圖解決的實際方程，那將會更有意義。

來源

2012-10-20 03:43:17 derekswanson08

Derek，在StackOverflow上，「答案」應該是真正回答問題的東西。尋求更多信息的回覆屬於對原始問題的評論。歡迎來到StackOverflow。 –

對不起，我沒有看到問題的任何地方的評論選項，雖然顯然我被允許評論我自己的答案... – derekswanson08

不幸的是，你需要一定的聲望（50）之前，你甚至可以評論其他人們的問題或答案，IIRC。所以這樣做很好。 – nneonneo

-1。＃IND0錯誤使用範圍庫塔方法在C

回答

相關問題