- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
I have two codes to reproduce the situation:
Code1
PROGRAM test REAL(KIND=8) dt, t, tmax INTEGER i, imax PARAMETER (tmax = 10.2D0, dt = 0.2D0, imax = 10) t = 0 i = 0 DO WHILE (t .LT. tmax) PRINT *, "before", t, i t = t + dt i = i + 1 PRINT *, "after", t, i END DO END PROGRAM test
Code 2
PROGRAM test REAL(KIND=8) dt, t, tmax INTEGER i, imax PARAMETER (tmax = 10.2D0, dt = 0.2D0, imax = 10) t = 0 i = 0 DO WHILE (i .LT. imax) PRINT *, "before", t, i t = t + dt i = i + 1 PRINT *, "after", t, i END DO END PROGRAM test
If I look at the output, the last 3 terms code 1 prints:
after 10.2000000000000 51 before 10.2000000000000 51 after 10.4000000000000 52
and the code 2
after 1.80000000000000 9 before 1.80000000000000 9 after 2.00000000000000 10
The code 2 is doing what it is supposed to, code 1 is doing one more loop it is not supposed to do. I believe the reason is the machine representation of numbers. But what do I do to get around it, to get a real "less than" condition?
Link Copied
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
You could determine the loop count before you start and use an integer counter.
E.G. nloops = NINT((tmax-tmin)/dt)
assuming tmin is set.
Testing reals for equality is always going to lead to grief unless you allow for the possibility of round off error.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Thank you for that. But then, there is an interesting (and for me embarassing). How do you determine the loop cound (I think you should have INT( (tmax-tmin)/dt ) +1 in your post, but anyway, this is valid only for the loop of the type where I wanted DO WHILE (t .LT. tmax) but this no longer works for the loop type DO WHILE (t .LE. tmax) and there is a difference from the above number of loops by 1 and this depends on dt by I have not been able to figure out how. It is always one less or the same as INT( (tmax-tmin)/dt ) +1
To clarify, I want the same number of loops as would have happened in the case t .LE. tmax in the case we were able to do precise arithmentics (no round off errors etc).
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
The phenomenon you are seeing is precisely the reason why real numbers are tricky to use as a loop variable. One improvement:
Use t = tmin + dt * (i-1)
instead of accumulating the variable t. Accumulation means that also the rounding error is accumulated.
If you seek a lower-than loop:
n = int( (tmax - tmin - 0.5d0*dt) / dt)
ought to do it
For a lower-equal:
n = int( (tmax - tmin + 0.5d0*dt) / dt)
Caveat: I have not tested these formulae in any formal sense, so I may be wrong :).
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
To avoid the accumulation of roundoff errors...
you avoid accumulation altogether.
PROGRAM test REAL(KIND=8) dt, t, tmax INTEGER i, imax PARAMETER (tmax = 10.2D0, dt = 0.2D0, imax = 10) t = 0 i = 0 DO WHILE (i*dt .LT. tmax) PRINT *, "before", t, i t = t + dt i = i + 1 PRINT *, "after", t, i END DO END PROGRAM test
Jim Dempsey
- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page