Untitled Document

The code

You will find below the fortran source that was used to process all the results presented in this report. It should be written in standard f77 but I used a f90 compiler to build the executable file.
I must apologise because the bassin of attraction section does not work... if anyone could fix the bug because I am affraid I am getting blind these days!

************************************************************************
************************************************************************
************************************************************************
         program main
      implicit none
************************************************************************
C     Declaration des variables
      REAL*8 x_old, y_old
      REAL*8 delta, epsx, epsy
      REAL*8 x_min, x_max, y_min, y_max
      REAL*8 nx, ny, l1, l2, truc

REAL*8 C_x, C_y
REAL*8 R

REAL*8 alpha, beta
INTEGER n_pts, n_iter

INTEGER option, defaut

10   WRITE(*,*) 'Choisir une option :'
      WRITE(*,*) '1-Attracteurs divers et varies'
      WRITE(*,*) '2-Points d''equilibre'
      WRITE(*,*) '3-Déformation d''un cercle'
      WRITE(*,*) '4-Sensibilite aux conditions initiales (SCI)'
      WRITE(*,*) '5-Tracer du bassin d''attraction'
      READ(*,*) option

************************************************************************
C     Attracteurs divers et varies
      if (option.eq.1) then
         WRITE(*,*) 'x0     y0     alpha   beta'
         WRITE(*,*) '0      0      1.4     0.3      Attracteur de Henon'
         WRITE(*,*) '0      0      -0.52   -0.998   Effondrement'
         WRITE(*,*) '0      0      -0.52   -1       Courbe fermée'
         WRITE(*,*) '0.132 -1.169 -0.651 -1       Comète...'
         WRITE(*,*) '0.132 -1.169 0.651   -1       Courbe fermée'
         WRITE(*,*) '-0.313 -0.166 -0.497 -0.999   Spirale à pôles'
         WRITE(*,*) '0.22   0.01   0.999   -0.999   Berlingot'
         WRITE(*,*) '0.44   0.114 1.801   0.064    Autre étrange...'
         pause
         WRITE(*,*) 'Entrer x0'
         READ(*,*) x_old
         WRITE(*,*) 'Entrer y0'
         READ(*,*) y_old
         WRITE(*,*) 'Entrer alpha'
         READ(*,*) alpha
         WRITE(*,*) 'Entrer beta'
         READ(*,*) beta
         WRITE(*,*) 'Enter le nombre d''iterations'
         READ(*,*) n_iter

call attract(x_old, y_old, alpha, beta, n_iter)
end if

************************************************************************
C     Points d'equilibre
      if (option.eq.2) then
         WRITE(*,*) 'Nous allons maintenant tester les
     &deux points d''equilibre'
         WRITE(*,*) 'Enter le nombre d''iterations'
         READ(*,*) n_iter
20      WRITE(*,*) 'Entrer les valeurs de alpha et beta'
         READ(*,*) alpha, beta

         delta = (1.0d0-beta) * (1.0d0-beta) + 4.0d0*alpha
         if (delta.le.0) then
            WRITE(*,*) 'Veuillez mieux choisir alpha et beta de sorte
     &qu''il existe des points d equilibre'
            WRITE(*,*) 'Notez que la condition alpha>0 est suffisante'
            GOTO 20
         end if

x_old = -(1.0d0 - beta + SQRT(delta)) / (2.0d0 * alpha)
y_old = beta * x_old

         truc = 0.5d0 * (1-beta+SQRT(delta))
         l1 = truc - (truc**2+beta)**(0.5d0)
         l2 = truc + (truc**2+beta)**(0.5d0)
         WRITE(*,*) 'Premier point'
         WRITE(*,*) 'premiere valeur propre: ', l1
         WRITE(*,*) 'premiere valeur propre: ', l2
         call attract(x_old, y_old, alpha, beta, n_iter)

x_old = -(1.0d0 - beta - SQRT(delta)) / (2.0d0 * alpha)
y_old = beta * x_old

         truc = 1-beta-SQRT(delta)
         l1 = truc - (truc**2+beta)**(0.5d0)
         l2 = truc + (truc**2+beta)**(0.5d0)
         WRITE(*,*) 'Deuxieme point'
         WRITE(*,*) 'premiere valeur propre: ', l1
         WRITE(*,*) 'premiere valeur propre: ', l2
         call attract(x_old, y_old, alpha, beta, n_iter)

end if

************************************************************************

      if (option.eq.3) then
         WRITE(*,*) 'Deformation d''un cercle de conditions initiales'
         WRITE(*,*) ''
         WRITE(*,*) 'Si vous voulez des valeurs par defaut, taper 1'
         WRITE(*,*) 'Sinon, taper 0'
         READ(*,*) defaut
         if (defaut.eq.1) then
            alpha = 1.4d0
            beta   = 0.3d0
C           Coordonnées du centre
            C_x    = 0.25d0
            C_y    = 0.0d0
C           Rayon
            R      = 0.1d0
C           Nombre de points
            n_pts = 5000
C           Nb d'iterations avant la sortie de l'image
            n_iter = 1
         else
            WRITE(*,*) 'Entrer alpha et beta'
            READ(*,*) alpha, beta
            WRITE(*,*) 'Entrer les coordonnees du centre du cercle'
            READ(*,*) C_x, C_y
            WRITE(*,*) 'Entrer le rayon du cercle'
            READ(*,*) R
            WRITE(*,*) 'Entrer le nombre de points'
            READ(*,*) n_pts
            WRITE(*,*) 'Entrer le nombre d''itérations avant affichage'
            READ(*,*) n_iter
         end if
         call cercle(alpha, beta, C_x, C_y, R, n_pts, n_iter)
      end if

************************************************************************
      if (option.eq.4) then
         WRITE(*,*) 'Entrer x0 du premier point '
         READ(*,*) x_old
         WRITE(*,*) 'Entrer y0 du premier point '
         READ(*,*) y_old
         WRITE(*,*) 'Entrer l''ecart en x entre les points '
         READ(*,*) epsx
         WRITE(*,*) 'Entrer l''ecart en y entre les points '
         READ(*,*) epsy
         WRITE(*,*) 'Entrer alpha et beta '
         READ(*,*) alpha, beta
         WRITE(*,*) 'Entrer le nombre d''itérations '
         READ(*,*) n_iter
         call sci(alpha, beta, x_old, y_old, n_iter, epsx, epsy)
      end if

************************************************************************
      if (option.eq.5) then
         WRITE(*,*) 'Entrer x_min, x_max et le nombre de points en x '
         READ(*,*) x_min, x_max, nx
         WRITE(*,*) 'Entrer y_min, y_max et le nombre de points en y '
         READ(*,*) y_min, y_max, ny
         WRITE(*,*) 'Entrer alpha et beta '
         READ(*,*) alpha, beta

call bassin(alpha, beta, x_min, x_max, y_min, y_max, nx, ny)
end if

************************************************************************
C     Verification que le choix saisi est valide !!!
      if (option.le.0) then
         WRITE(*,*) 'Veuillez entrer un choix valide '
         WRITE(*,*) ''
         GOTO 10
      end if

      if (option.ge.6) then
         WRITE(*,*) 'Veuillez entrer un choix valide '
         WRITE(*,*) ''
         GOTO 10
      end if

end program
************************************************************************
************************************************************************
************************************************************************

************************************************************************
C     Cette routine fait les iterations du modele de Henon avec les valeurs
C     Passees pour alpha, beta, x0, y0 et le nombre d'iterations
C     On a fixe un critere d'arret du calcul a 1E100
      SUBROUTINE attract(x_old, y_old, alpha, beta, n_iter)
      implicit none

      REAL*8 x_old, y_old, x_new, y_new
      REAL*8 alpha, beta
      INTEGER n_iter, i
      CHARACTER*20 resultats

      WRITE(*,*) 'Entrer le nom du fichier de stockage des resultats '
      READ(*,*) resultats
      OPEN(1,FILE=resultats)
      WRITE(1,*) x_old, y_old, 'Points de depart '

      do i=1,n_iter
         x_new = y_old + 1.0d0 - alpha * x_old * x_old
         y_new = beta * x_old
         WRITE(1,*) x_new, y_new

x_old=x_new
y_old=y_new

         if (x_new.gt.1.0d100) then
            WRITE(*,*) 'La coordonnee "X" a depasse 1E100... le calcul
     &a surement diverge !!! '
            GOTO 10
         end if
         if (y_new.gt.1.0d100) then
            WRITE(*,*) 'La coordonnee "Y" a depasse 1E100... le calcul
     &a surement diverge !!! '
            GOTO 10
         end if
      end do

10 CLOSE(1)
end subroutine

************************************************************************
C     Cette routine calcule la deformation d'un cercle de conditions
C     initiales.
      subroutine cercle(alpha, beta,C_x, C_y, R, n_pts, n_iter)
      implicit none

      REAL*8 C_x, C_y
      REAL*8 R
      REAL*8 x_old, y_old
      REAL*8 x_new, y_new
      REAL*8 beta, alpha
      INTEGER n_pts, n_iter
      INTEGER i, j

OPEN(1,FILE='cercle.xls')

      do i=1,n_pts
         x_old = C_x + R * COS(2.0*i*3.1416/n_pts)
         y_old = C_y + R * SIN(2.0*i*3.1416/n_pts)
         WRITE(1,*) x_old, y_old
         do j=1,n_iter
            x_new = y_old + 1 - alpha * x_old * x_old
            y_new = beta * x_old
            x_old = x_new
            y_old = y_new
         end do
         WRITE(1,*) x_new, y_new
      end do

CLOSE(1)
end subroutine

************************************************************************
C     Cette routine calcule les trajectoires de deux points initialement
C     voisins de epsilon.
      subroutine sci(alpha, beta, x_old1, y_old1, n_iter, epsx, epsy)
      implicit none

      REAL*8 x_old1, y_old1, x_new1, y_new1
      REAL*8 x_old2, y_old2, x_new2, y_new2
      REAL*8 epsx, epsy, ecart
      REAL*8 alpha, beta
      INTEGER n_iter, i
      CHARACTER*20 resultats

      x_old2 = x_old1 + epsx
      y_old2 = y_old1 + epsy
      WRITE(*,*) 'Entrer le nom du fichier de stockage des resultats'
      READ(*,*) resultats
      OPEN(1,FILE=resultats)
      ecart = SQRT((x_old1 - x_old2)**2 + (y_old1 - y_old2)**2)
      WRITE(1,*) i, ecart

      do i=1,n_iter
         x_new1 = y_old1 + 1.0d0 - alpha * x_old1 * x_old1
         y_new1 = beta * x_old1
         x_new2 = y_old2 + 1.0d0 - alpha * x_old2 * x_old2
         y_new2 = beta * x_old2
         ecart = SQRT((x_new1 - x_new2)**2 + (y_new1 - y_new2)**2)
         WRITE(1,*) i, ecart

         x_old1=x_new1
         y_old1=y_new1
         x_old2=x_new2
         y_old2=y_new2
      end do

CLOSE(1)
end subroutine

************************************************************************
C     Cette routine sert a tracer le bassin d'attraction. La methode
C     retenue consiste a dire qu'un point est en dehors du bassin
C     d'attraction si un des 100 premiers iteres de ce point a une
C     coordonnee qui depasse 1e100. Sinon il est considere comme
C     interieur au bassin.
      subroutine bassin(alpha, beta, x_min, x_max, y_min, y_max, nx, ny)
      implicit none

      REAL*8 x_min, x_max, y_min, y_max
      REAL*8 x_old, x_new, y_old, y_new
      REAL*8 delta_x, delta_y, x, y, nx, ny
      REAL*8 alpha, beta
      INTEGER n_iter, i, j, k
      CHARACTER*20 resultats

      n_iter = 100
      WRITE(*,*) 'Entrer le nom du fichier de stockage des resultats'
      READ(*,*) resultats
      OPEN(1,FILE=resultats)

delta_x = (x_max - x_min) / nx
delta_y = (y_max - y_min) / ny

      do i=1, nx+1
         x_old = x_min + (i-1) * delta_x
         do j=1, ny+1
           y_old = y_min + (j-1) * delta_y
           do k=1,n_iter
              x_new = y_old + 1.0d0 - alpha * x_old * x_old
              y_new = beta * x_old

              if (x_new.gt.1.0d10) then
                 GOTO 30
              end if
              if (y_new.gt.1.0d10) then
                 GOTO 30
              end if
              if (x_new.LT.-1.0d10) then
                 GOTO 30
              end if
              if (y_new.lt.-1.0d10) then
                 GOTO 30
              end if

              x_old = x_new
              y_old = y_new
            end do
           x = x_min + (i-1) * delta_x
           y = y_min + (j-1) * delta_y
           WRITE(1,*) x, y
30        x_old = x_min + (i-1) * delta_x
         end do
      end do

CLOSE(1)
end subroutine

"That's all folks" so please go back to the main page if there is anything you would like to check out... or have a look at what the other students did.