0

I was trying to solve a coupled first order differential equations using the 4 points Range-Kutta method. When outputting the values of m, I get the -1.#IND0 error. I am aware that this can be NaN, but it doesn't make sense to me, because the value of m should be increasing and I get -1IND0 in between valid values. Here is sample of my output:

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

and here is my code:

#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; }
Improve this question
7
  • -0.000000000000000012812899255404507 That's fully 16 significant figures. It's not obvious that that standard (double) literal will support every bit of that precision. Using -0.000000000000000012812899255404507L will get a bit more on platforms where long double is longer than double. (Which you might or might not want, depending...) Commented Oct 20, 2012 at 1:02
  • Hmmm...and you only defined pi to 10 significant figures, so the extra is in f is for naught. Consider const double pi = 4.0 * atan(1.0). Commented Oct 20, 2012 at 1:04
  • could you please try to compile the code with your compiler using the suggested edits? i can't seem to get it working. sorry i'm a beginner in c. Commented Oct 20, 2012 at 1:17
  • do you guys get the same errors? Commented Oct 20, 2012 at 1:58
  • BTW---My notes will not affect the problem you are having, I was just taking note of some things that you might want to be aware of going forward. But you need to get used to debugging...if you stick with programming you'll be doing it a lot from here on out. Commented Oct 20, 2012 at 3:19

1 Answer 1

Reset to default
1

I got nan's. compiled and run on 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

It's been a while since I studied Runge-Kutta... Looking at your code, I think that r is the independent variable, dr is the step size, and m is the dependent variable you are trying to solve for. I am confused what p is. Can you give us more details? It would make more sense if I could see the actual equations you are trying to solve.

Improve this answer
Sign up to request clarification or add additional context in comments.

3 Comments

Derek, on StackOverflow, an "Answer" is supposed to be something that actually answers the question. Replies that seek more information belong in the comments to the original question. Welcome to StackOverflow.
Sorry, I don't see a comment option anywhere on the question, although apparently I am allowed to comment on my own answer...
Unfortunately, you need a certain amount of reputation (50) before you can even comment on other people's questions or answers, IIRC. So doing it this way is fine.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.