floating point exception in C code !

Question

#include<stdio.h>
#include<math.h>

int main ()
{
    FILE *fp;
    fp=fopen("output","w");
    float t,y=0,x=0,e=5,f=1,w=1;
    for (t=0;t<10;t=t+0.01)
    {
        if( y==inf && y== nan) 
            break;
        fprintf(fp,"%lf\t%lf\n",y,x);
        y = y + ((e*(1 - x*x)*y) - x + f*cos(w*t))*t;
        x = x + y*t;
    }
    return 0;
}

Why is the ouput giving infinite and NAN values?

Answer 1

Perhaps the conditional statement should be if (y == inf || y == nan)? y cannot be both inf and NaN at the same time.

Answer 2

Because the break will never happen, as toft pointed out, and you never close the file, no data is written and then it all breaks down when the exception actually happens. Try to flush the data to the file in the for loop.

Answer 3

Comparing to inf or nan is done via isnan() and isinf() functions, not this way.
%lf is for double, not float.

And, for god's sake, fclose() your file! (first X lines are some meaningful numbers.)

Answer 4

Your calculation is blowing up. Just look at the values printed out for x and y and you will see they start to get very large and then turn info inf. Because your conditional is wrong, you wind up using inf in the calculation which turns into nan.

Answer 5

Be aware also that every comparison involving NaN returns false regardless of what operator you use (<, <=, etc) and what you compare it against.

Answer 6

I plugged your equations into Excel and when t gets to 0.39 x and y are 2E+270 and 5.2E+270 respetively. After that, they're too big for Excel to handle. Here is a trace of the other values if you're interested (lousy formatting, I know):

t   x       y
0.00    0       0
0.01    9.9995E-05  0.0099995
0.02    0.000719864 0.03099345
0.03    0.002688085 0.065607372
0.04    0.007431658 0.118589309
0.05    0.017321769 0.197802239
0.06    0.036281294 0.315992075
0.07    0.070850729 0.493849065
0.08    0.132072044 0.765266446
0.09    0.238828626 1.186184238
0.10    0.423641428 1.848128022
0.11    0.741634646 2.89084744
0.12    1.277397835 4.464693237
0.13    2.107123319 6.382503724
0.14    3.048840735 6.726552977
0.15    3.361369798 2.083527082
0.16    3.297988472 -0.396133288
0.17    3.231095608 -0.393487431
0.18    3.156811159 -0.412691383
0.19    3.07386543  -0.436556472
0.20    2.980485636 -0.466898968
0.21    2.874086586 -0.506662146
0.22    2.750700903 -0.56084401
0.23    2.603841953 -0.638517177
0.24    2.42203123  -0.757544679
0.25    2.18283838  -0.9567714
0.26    1.836632623 -1.331560601
0.27    1.255566799 -2.152095647
0.28    0.052253914 -4.297546016
0.29    -2.923971917    -10.2628477
0.30    -2.375067218    1.82968233
0.31    -1.600801299    2.497631996
0.32    0.082957072 5.261744912
0.33    4.774399642 14.21649263
0.34    -20.07952886    -73.09978971
0.35    3522.569432 10121.85417
0.36    -16277355468    -45214886085
0.37    1.64003E+30 4.43252E+30
0.38    -1.72156E+90    -4.53043E+90
0.39    2.0423E+270 5.2366E+270
0.40    #NUM!       #NUM!

So I guess the question is what is it this code is supposed to calculate?

Answer 7

Among other issues that others have addressed, the conditional break will never be triggered. y == nan is false for every value y, including NaN, so the condition is just (y == inf && false), which is, of course, false.

If you want to break when y is inf or nan, you should use:

if (isinf(y) || isnan(y)) break;

Or, if you want to use comparisons instead:

if (fabs(y) == inf || y != y) break;

(The fabs is included because you -- presumably -- also want to break if y is -inf.)

Answer 8

Others have pointed out that you are having nan/inf problems which is true, but here is how to fix your code to give you the results that I believe you are looking for.

Since no one else has really pointed it out (that I've noticed), you are trying to solve a system of differential equations using the Euler method. The coupled differential equations that you are solving are:

dy/dt = e*(1 - x * x) * y - x + f * cos(w * t)

dx/dt = y

However, your solution is faulty which gives huge numerical instability (and the wrong answer). These two lines:

    y = y + ((e*(1 - x*x)*y) - x + f*cos(w*t))*t;
    x = x + y*t;

should be:

    y = y + ((e*(1 - x*x)*y) - x + f*cos(w*t))*.01;
    x = x + y*.01;

where I have changed t to your delta t (time step) because that's what Euler's method calls for. I would make a new variable called delt or something like that so that you can easily change the time step. The solution is beautifully stable now and plotting x vs. t and y vs. t gives some very nice plots. I'd post them, but I have a feeling that this might be homework.

Also, if with different equations you need more stability, you can use smaller time steps or some better numerical ODE methods like Runge-Kutta or implicit methods.