1

I am writing a code for a Monte Carlo simulation in Fortran, but I am having a lot of problems because of the small numbers involved.

The biggest problem is that in my code particle positions are not updated; the incriminated code looks like this

x=x+step*cos(p)*sin(t)

with step=0.001. With this, the code won't update the position and I get a infinite loop because the particle never exits the region. If I modify my code with something like this:

 x=x+step

or 

x=x+step*cos(t)

there is no problem. So it seems that the product step*cos(t)*cos(p)(of the order 10**-4) is too small and is treated as zero.

x is of the order 10**4.

How do I solve this problem in portable way?

My compiler is the latest f95.

francescalus
  • 30,576
  • 16
  • 61
  • 96
mattiav27
  • 655
  • 2
  • 9
  • 27
  • 3
    Please give details of the values of the variables and their declarations. You don't care so much the magnitude of `y` in `x=x+y` but its size in comparison with the spacing around `x`. Without further details we can't be more precise. – francescalus Mar 26 '19 at 18:13
  • @francescalus `x` is of the order 10**4. I have edited the post. – mattiav27 Mar 26 '19 at 18:21
  • 1
    You are probably adding a tiny number of less than 1e-4 to a single precision real of size 1e4. This is similar to a human able to remember maximum 3 digits of a number and he adds 10 together with 0.000001. since he only can remember 3 digits, the result stays 10. – kvantour Mar 26 '19 at 18:24
  • Here it is not important, but next time it might be. f95 does not specify any particular compiler or its specific version. We can only guess that you are in fact using gfortran on Ubuntu Linux or something like that. – Vladimir F Героям слава Mar 26 '19 at 19:32

1 Answers1

2

Your problem is essentially the one of this other question. However, it's useful to add some Fortran-specific comments.

As in that other question, the discrete nature of floating point numbers mean that there is a point where one number is too small to make a difference when added to another. In the case of this question:

if (1e4+1e-4==1e4) print *, "Oh?"
if (1d4+1d-4==1d4) print *, "Really?"
end

That is, you may be able to use double precision reals and you'll see the problem go away.

What is the smallest number you can add to 1e4 to get something different from 1e4 (or to 1d4)?

print *, 1e4 + SPACING(1e4), 1e4+SPACING(1e4)/2
print *, 1d4 + SPACING(1d4), 1d4+SPACING(1d4)/2
end

This spacing varies with the size of the number. For large numbers it is large and around 1 it is small.

print*, EPSILON(1e0), SPACING([(1e2**i,i=0,5)])
print*, EPSILON(1d0), SPACING([(1d2**i,i=0,5)])
end
francescalus
  • 30,576
  • 16
  • 61
  • 96