Nested iterators produce only zeros as data (gfortran)

[agsiar]

Hello. I'm supposed to plot a bifurcation diagram for the logistic map. What I'm trying to do is calculate 1128 iteration of y=a*y*(1-y), scrap the first 1000 and write the last 128 points of data in the LogisticMapDiagram.dat file and then plot that with gnuplot, but the value of y stays in zeroes all the time.
The code:

program BifurcationDiagram
implicit none
real y, a
integer i
open (unit=2, file="LogisticMapDiagram.dat")
write (2,*) "# a"," y"
y=0.4
a=0.0
print *, "beginning do"
do i=1,1128,1
do a=0.,4.,0.01
y=a*y*(1-y)
if (i .ge. 1001) then
write (2,*) a, y
end if
end do
end do
print *, "i'm done"
end program​

I've noticed that running a simplified version of the above program yields part of the results that I'm looking for:

program Test
real y, a
y = 0.4
a=3.7
open (unit=2, file="Test.dat")
do i=1,1128,1
y=a*y*(1-y)
if (i .ge. 1001) then
write (2,*) a, y
end if
end do
print*, "i'm done"
end program​

I've also noticed that if I don't specify a value of a and change do i=1,1128,1 for do a=0.,4.,0.01 in the immediately above code, I go back to getting nothing but zeroes for y. Thing is, I really do need to sweep the value of a, preferably in steps of 0.001 instead of steps of 0.01.

I've been told that fortran doesn't like using real numbers in its do statement, but I did another test tabulating sin(x) using x as the counter for the iterator and it worked. The test in question:

program TabulatingSine
real x
open (unit=2, file="SenoTabulado.dat")
write (2,10) "#x","seno(x)"
10 format (a,4x,a)
do x=0.,6.28318530718,0.001
write (2,20) x, sin(x)
20 format (f10.5,4x,f10.5)
end do
print *, "the program did something"
print *, "bye."
end program​

I've thought of not using a as the counter, but I can't come up with any way to make it work and I honestly don't know what to do here anymore...

P.S: The indentation got kinda screwed after CP-ing the code...

EDIT:
Nope, changing the counter from a real number to an integer and make the value of a vary within that do makes no difference... Code:

program DiagramaBifurcacion
implicit none
real y, a
integer i, j
!double precision y, a
open (unit=2, file="MapaLogisticoDiagrama.dat")
write (2,*) "# a"," y"
y=0.4
a=-0.01
print *, "beginning do"
do j=0,400,1
a=a+0.01
do i=1,1128,1
y=a*y*(1-y)
if (i .ge. 1001) then
write (2,*) i, a, y
end if
end do
end do
print *, "i'm done"
end program​

Now I'm truly at a complete loss of what to do...

 

[gsal]

First, please use code and /code tags enclosed in square brackets '[' and ']' surrounding source code...makes it easier for everybody to read.

It is best to only use integer number in loop variables. You have the right idea in your last example:

do j=0,400
  a=a+0.01
end do

or, simply

do j=0,400
  a=j/100.0
end do

may be more accurate in certain instances, because adding small numbers to large ones sometimes gets lost...you are o.k., now, though.

Noticed in the second example above that in the division I do use a real "100.0" and not an integer "100" The thing is that if I had used an integer then it becomes integral division which always yields zero when the numerator is smaller than the denominator.

Anyway...one thing you need to observe is that if you start with a=0.0 and evaluate "y = a*y*(1-y)", then, no matter what the initial value of y is, it will end up zero right away; then, because y is a multiplier on the right-hand-side, it will continue to be zero no matter what the value of 'a' is...

So, in your very first piece of code, the very first value for 'a' in the inner loop is zero...that's it, you are done...now, you will never be able to get out of there...y will be zero forever! Follow? So, your program is doing exactly what you have programmed...that's a problem with most computers, they do exactly what you tell them to do.

So, first fix is: do not use a=0...you know what the result is...zero.

Second, I don't know why, but I would switch the order of the loops...

do j=1,400
    a = j/100.0
    do i = 1,1128
        y = a*y*(1-y)
    end do
end do

That way, your consecutive 1128 iterations kind of make more sense as the seed (a) was set at the beginning of this set of iterations and that's it; then, re-seeding for the next set and so on.

 

[agsiar]

Thanks, I didn't know of the code wrapping.
I've tried what you mentioned and sadly it doesn't help. Even before the change, the problem wasn't that a always had the same value (zero), it was that y changed from it's initial value (0.4) to zero and stayed there. Here's the modified code, which is pretty much the same:

    program DiagramaBifurcacion
        implicit none
        real y, a
        integer i, j
        !double precision y, a
        print *, "abro archivo ''MapaLogisticoDiagrama.dat'' en unidad 2"
        open (unit=2, file="MapaLogisticoDiagrama.dat")
        print *, "escribo ''# a'' y ''    y'' en archivo"
        write (2,*) "# a","    y"
        y=0.4
        print *, y
        print *, "empiezo el do"
        do j=1,400
            a = j/100.0
            do i = 1,1128
                y = a*y*(1-y)
                if (i .ge. 1001) then
                    write (2,*) a, y
                end if
            end do
        end do
        print *, "ya terminé"
    end program

This is a chunk of the data it produces:
3.4900000 0.0000000
3.4900000 0.0000000
3.4900000 0.0000000
3.4900000 0.0000000
3.4900000 0.0000000
3.4900000 0.0000000
3.4900000 0.0000000
3.4900000 0.0000000
3.5000000 0.0000000
3.5000000 0.0000000
3.5000000 0.0000000
3.5000000 0.0000000
3.5000000 0.0000000
3.5000000 0.0000000
3.5000000 0.0000000
3.5000000 0.0000000
3.5000000 0.0000000
3.5000000 0.0000000
3.5000000 0.0000000
3.5000000 0.0000000
3.5000000 0.0000000
3.5000000 0.0000000
As you see, the value of a does change, but y stays the same throughout the whole execution.

This is a chunk of data produced by the code previously shown with a constant value of a:
3.5000000 0.82694072
3.5000000 0.50088418
3.5000000 0.87499726
3.5000000 0.38281968
3.5000000 0.82694072
3.5000000 0.50088418
3.5000000 0.87499726
3.5000000 0.38281968
3.5000000 0.82694072
3.5000000 0.50088418
3.5000000 0.87499726
3.5000000 0.38281968
3.5000000 0.82694072
3.5000000 0.50088418
3.5000000 0.87499726
3.5000000 0.38281968
3.5000000 0.82694072
3.5000000 0.50088418
3.5000000 0.87499726
3.5000000 0.38281968
3.5000000 0.82694072
3.5000000 0.50088418
3.5000000 0.87499726
3.5000000 0.38281968
As you see a doesn't change, but it produces data that can be plotted with gnuplot. Making a program and manually change the value of a to make different files doesn't sound like the best idea, that's why I'd rather make all data fall into one file without me having to change from 1.01 to 1.02 all the way to 4. But I don't know how to make gfortran do that or gfortran is being derp...

Edit:
Tried silverfrost's ftn95 and had the same issue. So I guess it has something to do with the program and not the compiler.

Edit2:
It's fixed! and now when plotted gives me the bifurcation diagram I wanted! A friend online helped me see the error that I was oblivious of for almost a whole day! Here's the final code:

    program DiagramaBifurcacion
        implicit none
        real y, a
        integer i, j
        !double precision y, a
        print *, "abro archivo ''MapaLogisticoDiagrama.dat'' en unidad 2"
        open (unit=2, file="MapaLogisticoDiagrama.dat")
        print *, "escribo ''# a'' y ''    y'' en archivo"
        write (2,*) "# a","    y"
        print *, y
        print *, "empiezo el do"
        do j=1,400
            y=0.4
            a = j/100.0
            do i = 1,1128
                y = a*y*(1-y)
                if (i .ge. 1001) then
                    write (2,20) a, y
20                  format (f10.5,4x,f15.9)
                end if
            end do
        end do
        print *, "ya terminé"
    end program

The change was to move the y=0.4 inside the first loop, and then it worked! :)
yeah I know, I'm a noob.