Université de Rennes I

Campus de Beaulieu

35042 Rennes Cedex


Magistère première année

1996-1997


Automates cellulaires

stochastiques

Transitions de phases

dues à la propagation de défauts

Céline Roget et Mélissa Inglart

Travail dirigé par D. Petritis


Introduction

Ce projet informatique a pour but d’étudier numériquement le comportement asymptotique d’automates cellulaires stochastiques. On commencera bien sûr par définir ce qu’est un automate cellulaire stochastique, l’objet précis de notre étude puis nous vous présenterons le programme fortran inhérent à l’étude numérique de ces automates suivant deux axes : la densité asymptotique et la mémoire des conditions initiales traduite par la distance de Hamming. Par la suite, nous nous attacherons à l’interprétation des résultats obtenus, suivie d’une comparaison avec le ferromagnétisme.

I. Définition des Automates Cellulaires Stochastiques

a)Automates cellulaires déterministes

Un automate cellulaire déterministe est défini par un quadruplet (X, M, f, s(0)) où

X est l’ensemble des suites doublements infinies d’éléments dans {0,1}

M est un entier fixé

f est une fonction de {0,1}(2M+1) dans {0,1}

s(1) une condition initiale c’est à dire une suite de X initiale.

On définit l’évolution temporelle de cet automate par :

sn(t+1)=f (sn-M(t), ... , sn+M(t))

ainsi l’état du site n au temps t+1 dépend de l’état des sites n-M, ... , n+M au temps t. Dans notre étude, nous nous sommes limités à M=1.

Dans toute la suite, on indexe l’espace par i, i appartenant à  l’ensemble des entiers de 0 à L-1, et l’index d’espace est pris modulo L.

D’autre part, on considère dans toute la suite que la valeur au point n et au temps t+1 ne dépend que de la valeur des points n-1 et n+1 au temps t. Autrement dit, pour calculer par exemple la valeur en i=1 au temps t+1, on applique f à (sL(t), s2(t)), de même pour i=L, on aura sL(t+1)=f (sL-1(t), s1(t)).

Dans le cas d’un automate cellulaire déterministe, la valeur de f dans les quatre cas est fixée, par exemple :

f (0,0)=0

f (0,1)=1

f (1,0)=1

f (1,1)=1.

Finalement se donner un automate cellulaire, c’est se donner un réseau composé de L sites s(i,t) avec s(i,t) appartenant à {0,1} et regarder son évolution par pas de temps.

Rmq : on se limite ici à la dimension 1 ; mais toutes les conclusions pourront être étendues à une dimension n, n 1.

Explication de l’évolution sur un exemple.

Prenons  f (0,0)=0,  f (0,1)= f (1,0)=1,  f (1,1)=1.

Un point du réseau représentera une cellule, vivante si le point vaut 1, morte sinon.

s(1)   _V___J___J___V___J___V___J___V___J___V_

             0       1       1      0       1       0       1       0       1       0  

s(2)   _J___J___J___J___V___J___J___J___J___V_

             1       1       1       1       0       1       1       1        1      0  

s(3)   _J___J___J___J___J___J___J___J___J___J_

             1       1       1       1       1       1        1      1        1       1  

L’état à t+1 dépend  de l’état des voisins à t. On a ici un genre de cellule qui meurt vite de solitude ( f (0,0)=0), mais elles ont un fort instinct grégaire : elles aiment la vie de groupe ( f (0,1)= f (1,0)=1,  f (1,1)=1) ! Remarquons d’autre part que ce dernier état est stable vu la fonction  f, toutes ces cellules sont désormais immortelles !

b)Automates cellulaires stochastiques

Dans ce projet, nous étudions en fait les automates cellulaires stochastiques, en effet f n’est pas déterminée que par ses voisins les plus proches ; on définit des probabilités pour que f donne 1 en fonction de ses voisins :

P ( f(0,0)=1 ) = p0

P ( f(0,1)=1 ) = p1

P ( f(1,0)=1 ) = p2

P ( f(1,1)=1 ) = p3

et on en déduit la probabilité que f =0 par la loi de normalisation. Donc pour trouver la valeur du site i au temps t+1 on compare un nombre tiré aléatoirement entre 0 et 1 à la probabilité d’avoir pour résultat 1.

Ainsi, pour le site i au temps t+1, si on a s(i-1,t) = s(i+1,t) = 0 et r le nombre tiré, on obtient  f(0,0) = 1 si r < 0,  f(0,0) =0 sinon.

Dans le cadre de notre étude, on considère désormais que p0 = 0 et p1 = p2 (cas symétrique).

Exemple muni de probabilités.

Prenons p0 = 0, p1 = p2 = 1, p3 =0,2.

s(1)   _V___V___V___J___J___J___V___J___V___V_

             0      0       0       1       1       1       0       1       0       0 

s(2)   _V___V___J___J___V___J___V___V___J___V_

             0      0       1       1       0       1       0       0       1       0 

s(3)   _V___J___J___J___V___V___J___J___V___J_

             0      1       1        1      0       0       1       1        0      1 

s(4)   _V___J___V___J___J___J___J___J___V___V_

             0      1       0       1       1        1       1       1       0       0 

s(5)   _J___V___V___J___V___V___J___J___J___V_

             1       0      0       1       0       0       1       1       1       0 

c)Position du problème à résoudre

Le processus vu plus haut est utilisé pour la modélisation de plusieurs systèmes en recherche physique et biologique. Domany et Kinzel ont introduit cette schématisation. Dans la suite, nous étudierons le modèle de Domany et Kinzel utilisé en thermodynamique. Ce modèle montre une transition de phases d’une phase ordonnée à une phase désordonnée pour p0 = 0. La phase ordonnée est définie par :

  s(i)=0

On dit que l’état est adsorbé.

Le paramètre d’ordre est la densité asymptotique m = lim lim m(t,L) avec

Là+ tà+

m (t,L) = s(i,t).

On obtient une surface critique qui sépare les régions où m=0 et m>0. Remarquons que m=0 si et seulement si un nombre fini de sites du réseau sont à l’état 1.

Un autre axe d’étude de ce modèle est d’observer la propagation de défauts. Cette dernière consiste à considérer deux répliques s et h suivant le même modèle, ici celui des automates cellulaires stochastiques, mais ayant des conditions initiales différentes. On observe la différence entre les deux évolutions, qui est donnée au temps t et au site i par

h(i,t) = s(i,t) h(i,t)

 où  est la somme modulo 2 ; h est en fait le défaut étudié entre deux répliques (s et h). Dans le cas des transitions de phases dues à la propagation de défauts, le paramètre d’ordre est la distance de Hamming asymptotique 

H = lim lim H(t,L)

Là+ tà+

avec

H (t,L) = h(i,t),

cette dernière somme étant la somme usuelle.

De même il existe donc une surface critique qui sépare les régions où H=0 et celle où H>0.

Dans le cas où le défaut disparaît, c’est-à-dire H=0, toutes les répliques suivent asymptotiquement la même trajectoire, qui est indépendante des conditions initiales. Si on est dans ce cas, on peut dire que le système n’a pas mémoire de ses conditions initiales puisque, quelles qu’elles soient, on arrive toujours asymptotiquement à la même forme du réseau. Par contre si H est strictement positif, l’état asymptotique du système dépend de ses conditions initiales.

Nous allons étudier les relations entre la transition de phase due à la densité et celle due à un défaut, dans le modèle de Domany-Kinzel.

Tout d’abord nous allons exposer notre programme servant à l’obtention du diagramme de phases, en fonction de p3  et de p1 dans le cas où p0  =0 et p2 = p1. Puis nous ferons l’analyse de ce diagramme en expliquant les résultats obtenus.

II. Présentation du programme

a)But du programme, exploitation graphique

Il fallait faire une étude numérique du problème, retrouver expérimentalement ce fameux diagramme des phases, aussi le langage fortran s’est-il imposé à nous. Nous avons d’abord divisé l’étude en deux axes principaux : la densité asymptotique d’une part, la distance de Hamming (inhérente à la mémoire des conditions initiales) d’autre part.

Les deux études ont été menées sur le même schéma. D’abord créer un damier de points, de coordonnées p1 et p3, ces deux quantités prenant « iter » valeurs différentes réparties uniformément sur [0 ;1] ; il s’agit donc de calculer les coefficients de la matrice correspondante, de dimensions iter*iter, coefficients représentant soit la densité asymptotique, soit la distance de Hamming pour la fonction  f alors considérée.

Ainsi à chaque point du damier correspond un couple (p1, p3). Notre étude se limitant à p0 =0 et p1 =p2 (cas symétrique), ce couple caractérise la fonction  f correspondante, qui sera modélisée dans le programme par le tableau (p0, p1, p2, p3 ).

 On fixe au préalable une condition initiale qui sera commune à chaque point du damier, puis pour chaque point du damier,  f étant ainsi définie, on calcule la distance asymptotique en faisant d’abord évoluer l’automate cellulaire du temps t=1 au temps T très grand puis en appliquant la formule m(t,L).

Quant à la mesure de la distance de Hamming, il suffit de fixer au préalable deux conditions initiales, et de faire évoluer les deux automates pour chaque point du damier et d’en calculer la distance de Hamming.

On obtient deux jolies matrices qu’il faudra exploiter graphiquement ; transformer la matrice en un tableau où les coefficients « non nuls » (supérieurs à un e donné) seront représentés par « * » (par exemple) et les coefficients « nuls » (inférieurs à e) par un « o ». Se dessine alors des zones, les phases, et des frontières, les transitions de phases. Notez que le traitement graphique n’est pas effectué dans le programme suivant, il a été effectué avec Matlab.

b)Listing fortran du programme

*----------------------------------------------------------------------

*

            program projet

*

*          Ce programme permet l'etude de la densite asymptotique et de la

*          distance de Hamming pour t, L, iter (precision du damier

*          de resultats) fixees par l’utilisateur ; les conditions intiales étant

*          générées aleatoirement et automatiquement.

*          ENTREE

*                 t : temps final "grand"

*                 L : longueur du tableau "grand"

*              iter : nombre de valeurs prises par p1 et p3

*               Yd1 : condition initiale 1

*               Yd2 : condition initiale 2

*          SORTIE

*              MatM : matrice iter*iter des resultats concernant la

*                        densite asymptotique

*              MatH : matrice iter*iter des resultats concernant la

*                        distance de Hamming

*

            integer t, L, iter, Yd1(5000), Yd2(5000)

            integer mresult, hresult, i, j, IX

            real MatM (21,21), MatH(21,21)

            real  U

*

            parameter (mresult=11, hresult=12)

*

*          Ouverture fichiers

            open (mresult, FORM = 'FORMATTED', STATUS = 'unknown',

     +        FILE='mresult')

            open (hresult, FORM = 'FORMATTED', STATUS = 'unknown',

     +        FILE='hresult')

*

*          Initialisation du generateur aleatoire

            write(*,*) 'Donner un nombre entier quelconque. '

            read (*,*) IX

*

*          Initialisation par l'utilisateur des differents parametres

            write(*,*) 'Donner t. '

            read (*,*) t

            write(*,*) 'Donner L (=<5000). '

            read (*,*) L

            write(*,*) 'Donner iter (=<21). '

            read (*,*) iter

*

*          Initialisation automatique aleatoire des conditions initiales

            do 10 i=1, L

               call unif(IX, U)

               if (U .LE. 0.5) then

                        Yd1(i)=0

               else

                        Yd1(i)=1

               endif

10       continue

            do 20 i=1, L

               call unif(IX, U)

               if (U .LE. 0.5) then

                        Yd2(i)=0

               else

                        Yd2(i)=1

               endif

20       continue

*

*          Lancement des etudes

            call EtudierM (t, L, Yd1, iter, MatM, IX)

            call EtudierH (t, L, Yd1, Yd2, iter, MatH, IX)

*

*          Recuperation des resultats

            do 30 i=1, iter

               do 40 j=1, iter

100           format (F5.3)     

                  write (mresult,100) MatM(i,j)

                  write (hresult,100) MatH(i,j)

40          continue        

30       continue

*

*          Fermeture Fichiers

            close(11)

            close(12)

*

            end

*

*----------------------------------------------------------------------

*

            subroutine EtudierH (t, L, Yd1, Yd2, iter, MatRes, IX)

*

*          Cette procedure permet de creer la matrice resultat de l'etude

*          des mesures de Hamming, indicateur de memoire des conditions

*          initiales, pour un damier de p1, p3

*          ENTREE t : temps final "grand"

*                 L : longueur du tableau "grand"

*               Yd1 : condition initiale fixee pour la premiere replique

*               Yd2 : de meme pour la seconde replique

*              iter : dimension de la

*                IX : parametre du generateur aleatoire

*          SORTIE MatRes : matrice des resultats (iter)*(iter)

*

            integer t, L, Yd1(L), Yd2(L), iter, i, j, k, Y1(5000), Y2(5000)

            integer IX

            real f(4),res, H,  MatRes(21,iter)

*

            f(1)=0.

*          On fait varier p1

            do 10 i=1, iter

*             On fait varier p3

               do 20 j=1, iter

*                     On met a jour Y1 et Y2

                        do 30 k=1, L

                           Y1(k)=Yd1(k)

                           Y2(k)=Yd2(k)

30                   continue

*                     Mettre a jour f

                        res = i-1

                        f(2)=res/(iter-1)

                        f(3)=f(2)

                        res = j-1

                        f(4)=res/(iter-1)

                        call calH(t,L,Y1,Y2,f,IX, H)

                        MatRes (i,j) = H

20          continue

10       continue

*

            end

*

*----------------------------------------------------------------------

*

         subroutine EtudierM (t, L, Yd, iter, MatRes, IX)

*

*          Cette procedure permet de creer la matrice resultat de l'etude

*          des densites pour un damier de p1, p3

*          ENTREE t : temps final "grand"

*                 L : longueur du tableau "grand"

*                Yd : condition initiale fixee

*              iter : dimension de la matrice resultat

*                IX : parametre du generateur aleatoire

*          SORTIE MatRes : matrice des resultats (iter)*(iter)

            integer t, L, Yd(L), iter, i, j, k, Y(5000), IX

            real f(4), res, m, MatRes(21,iter)

*

            f(1)=0.

*          On fait varier p1

            do 10 i=1,iter

*             On fait varier p3

               do 20 j=1, iter

*                     Mettre a jour Y

                        do 30 k=1, L

                           Y(k)=Yd(k)

30                   continue

*                     Mettre a jour f

                        res = i-1

                        f(2)=res/(iter-1)

                        f(3)=f(2)

                        res = j-1

                        f(4)=res/(iter-1)

                        call calm(t,L,Y,f,IX,m)

                        MatRes (i,j)=m

20          continue

10       continue

*

            end

*                    

*----------------------------------------------------------------------

*

         subroutine calm(t, L, Y, f, IX, m)

*

*          Cette fonction permet de calculer la densite d'une replique

*          ENTREE t : temps final

*                 L : longueur du tableau

*                 Y : replique

*                 f : fonction de transfert caracterisee par 4 probas

*                IX : parametre du generateur aleatoire

*          SORTIE m : densite de Y

*          ! Y a change

*

            integer t, L, Y(L), i, IX

            real f(4), res, m

*

*          Faire evoluer la replique Y de temps=1 a temps=t

            call Evoluer (t, L, Y, f, IX)

*         

*          Calcul de la densite

            res = 0.          

            do 10 i=1, L

               res = res + Y(i)

10       continue

            m = res/L

*

            end

*

*----------------------------------------------------------------------

*

         subroutine calH(t, L, Y1, Y2, f, IX, H)

*

*          Cette fonction permet de calculer la distance de Hamming entre

*          deux repliques (meme fonction d'evolution, point de depart

*          different).

*          ENTREE t : temps final

*                 L : longueur du tableau

*                Y1 : premiere replique

*                Y2 : seconde replique

*                 f : fonction de transfert caracterisee par 4 probas

*                IX : parametre du generateur aleatoire

*          SORTIE H : distance de Hamming

*

            integer t, L, Y1(L), Y2(L), i, IX

            real f(1:4), res, H

*

*          Faire evoluer les repliques Y1 et Y2 de temps=1 a temps=t

            call Evoluer (t, L, Y1, f, IX)

            call Evoluer (t, L, Y2, f, IX)

*         

*          Calcul de la distance de Hamming

            res = 0.

            do 10 i=1, L

               res = res + abs( Y1(i)-Y2(i) )

10       continue

            H = res/L

*

            end     

*

*----------------------------------------------------------------------

*

         subroutine Evoluer (t, L, Y, f, IX)

*

*          Cette procedure permet de faire evoluer l'automate stochastique

*          de temps=0 a temps=t

*          ENTREE t : integer

*                 L : integer, longueur du tableau

*                 Y : automate represente dans un tableau de longueur L

*                 f : fonction d'evolution

*                IX : parametre du generateur aleatoire

*          SORTIE Y

*

            integer t, L, Y(L), temps, i, p, c, r, IX

            real f(4)

*

            do 10 temps=2, t

               p = Y(L)

               r = Y(1)

               do 20 i=1, L-1

                        c=Y(i)

                        Y(i)=calf(p,Y(i+1),f,IX)

                        p=c

20          continue

               Y(L)=calf(p,r,f,IX)

10       continue

*

            end

*

*----------------------------------------------------------------------

*

         function calf(p,s,f,IX)

*

*          Cette fonction permet de calculer f(p,s) grace aux quatre

*          probabilites de transition.

*          Le parametre du generateur aleatoire est deja initialise

*          dans le programme principal.

*

            integer p, s, IX

            real f(4), U

*

            call UNIF (IX, U)

            if (p .EQ. 0) then

               if (s .EQ. 0) then

                        if (U .LE. f(1)) then

                           calf=1

                        else

                           calf=0

                        endif

               else

                        if (U .LE. f(2)) then

                           calf=1

                        else

                           calf=0

                        endif

              endif

            else

               if (s .EQ. 0) then

                        if (U .LE. f(3)) then

                           calf=1

                        else

                           calf=0

                        endif

               else

                        if (U .LE. f(4)) then

                           calf=1

                        else

                           calf=0

                        endif

              endif

            endif

*

            end

*

*----------------------------------------------------------------------          *                    

            SUBROUTINE UNIF (IX, U)

*         

*          Cette procédure permet d’obtenir un nombre « aleatoire »

*          entre 0 et 1, U, de manière uniforme.

*          IX esl le parametre du generateur aleatoire

*

            K1 = IX / 127773

            IX = 16807*(IX-K1*12773)-K1*2836

            IF (IX .LT. 0) IX = IX + 2147483647

            U= IX * 4.656612875E-10

            END

*

*----------------------------------------------------------------------

c)Etude succinte de la complexité

Pour chaque étude (« EtudierH » et « EtudierM »), on a mis en place un damier de iter² points ; pour chacun de ces points, nous avons du faire évoluer un ou deux automates de longueur L sur (t-1) pas de temps, donc appeler L*(t-1) ou 2*L*(t-1) « calf », effectuer 2*L*(t-1) ou 4*L*(t-1) affectations, puis calculer m ou H, c’est-à-dire procéder à une somme plus ou moins complexe suivie d’une division, comptant tantôt L+1, tantôt 3*L+1 opérations arithmétiques élémentaires ; de plus « calf » nécessite 3 tests et 9 opérations arithmétique (on inclut les opérations nécessitées par « unif »).

Tout cela fait beaucoup, pour l’étude de H, on obtient iter²*(2*L*(t-1)*13+3*L+1). Finalement une complexité globale en O(iter²*t*L).

Nous prendrons dans l’ensemble des calculs qui suivent, iter=21. Dès lors prendre plusieurs milliers pour t et pour L occasionne une quantité de calculs faramineuse (pour t=5000 et L=1000, de l’ordre de 1e11 opérations arithmétiques et logiques élémentaires), d’où finalement un temps de calcul conséquent.

III . Densité asymptotique : résultats quantitatifs

Nous avons lancé le programme pour t=5000, L=1000 et iter=21.

La matrice obtenue se traduit par : chaque ligne correspond à une probabilité p1 différente, p1 croissant par pas de 0.05 de 0 à 1 et chaque colonne correspond à une probabilité p3 différente, p3 croissant par pas de 0.05 de 0 à 1.

Colonne de 1 à 11.

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

On a ainsi le coefficient 0.000 en gras qui correspond à la densité m(5000,1000) pour f caractérisée par p1=0.05 et p3=0. De même le  coefficient 0.197, en gras, correspond à m(5000,1000) pour f caractérisée par p1=0.75 et p3=0.5.

 
0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.009 0.000 0.000 0.000 0.033 0.019 0.017 0.024 0.051 0.197

0.015 0.077 0.189 0.251 0.293 0.231 0.316 0.372 0.380 0.421 0.466

0.367 0.361 0.386 0.424 0.471 0.456 0.456 0.478 0.532 0.544 0.556

0.453 0.440 0.417 0.460 0.497 0.476 0.513 0.516 0.555 0.590 0.605

0.465 0.479 0.496 0.517 0.542 0.545 0.567 0.589 0.593 0.627 0.662

0.506 0.509 0.528 0.553 0.525 0.544 0.618 0.615 0.621 0.641 0.675

Colonne de 12 à 21.

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.538

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.017 1.000

0.000 0.000 0.000 0.000 0.000 0.000 0.010 0.059 0.404 1.000

0.000 0.000 0.000 0.000 0.025 0.152 0.288 0.543 0.792 1.000

0.018 0.049 0.079 0.170 0.387 0.482 0.631 0.755 0.866 1.000

0.262 0.398 0.441 0.518 0.570 0.684 0.718 0.806 0.896 1.000

0.492 0.527 0.580 0.639 0.662 0.716 0.775 0.867 0.929 1.000

0.601 0.615 0.634 0.660 0.720 0.746 0.795 0.849 0.929 1.000

0.614 0.654 0.681 0.720 0.754 0.792 0.861 0.898 0.959 1.000

0.678 0.683 0.697 0.735 0.767 0.821 0.833 0.899 0.946 1.000

0.699 0.722 0.726 0.781 0.811 0.848 0.883 0.938 0.954 1.000

IV . Distance de Hamming : resultats quantitatifs

Nous avons effectué plusieurs essais, le premier avec t=5000 et L=1000, le second avec t=30000 et L=1000 enfin le troisième avec t=15000 et L=5000. Notons que chacune de ces applications comportent énormément de calculs, donc le temps CPU de calcul peut devenir lui aussi long (quelques heures dans notre cas).

Ce qui nous a donné trois matrices différentes, nous avons constaté une légère diminution des coefficients non nuls mais trop légère pour sentir l’état d’équilibre tout proche. Nous présumons donc que cet état d’équilibre, approximation du comportement asymptotique de H, fonction de p1 et p3, ne peut être atteint que pour t et L beaucoup plus grand.

Voici donc la matrice correspondante à l’application t=5000 et L=1000, très similaire aux autres.

Colonne 1 à 11.

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.005 0.000 0.000 0.024 0.017 0.048 0.069 0.067 0.046 0.313

0.193 0.169 0.324 0.378 0.422 0.350 0.424 0.461 0.502 0.492 0.500

0.441 0.486 0.450 0.490 0.497 0.493 0.513 0.512 0.510 0.529 0.464

0.463 0.489 0.506 0.480 0.487 0.485 0.510 0.500 0.496 0.487 0.463

0.491 0.475 0.513 0.493 0.490 0.467 0.486 0.484 0.504 0.483 0.466

0.470 0.514 0.483 0.489 0.501 0.495 0.461 0.506 0.466 0.445 0.423

Colonne 12 à 21.

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.786

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.058 0.000

0.000 0.000 0.000 0.000 0.000 0.000 0.012 0.124 0.516 0.000

0.000 0.000 0.000 0.000 0.062 0.236 0.417 0.520 0.346 0.000

0.086 0.073 0.125 0.314 0.453 0.511 0.481 0.407 0.273 0.000

0.376 0.477 0.485 0.498 0.453 0.403 0.392 0.336 0.210 0.000

0.490 0.492 0.453 0.438 0.428 0.384 0.351 0.300 0.162 0.000

0.477 0.493 0.434 0.435 0.365 0.340 0.323 0.244 0.143 0.000

0.464 0.446 0.409 0.398 0.356 0.315 0.316 0.216 0.137 0.000

0.455 0.426 0.410 0.377 0.344 0.312 0.261 0.195 0.123 0.000

0.436 0.377 0.392 0.366 0.309 0.314 0.243 0.177 0.090 0.000

V . GRAPHIQUES

Nous avons effectué le traitement graphique avec Matlab.

Il s’agit de prendre chaque coefficient de la matrice résultat (représentée par un vecteur dans les fichiers hresult et mresult), de le comparer à un epsilon « epsil », nous avons d’ailleurs pris à chaque fois 0.05, et suivant ce test placer un « o » rouge ou une étoile bleu dans le graphe dont la taille est déterminée par iter : graphe [0 ; 1]*[0 ; 1], p1 en abscisse et p3 en ordonnée.

Pour cela, après avoir effectué l’application Fortran, nous chargeons les fichiers obtenus sous Matlab qui les considérera dès lors comme des vecteurs.

Le traitement avec la fonction graphe est alors possible.

function graphe(mat, iter, epsil)

hold on ;

for i=1:iter

   for j=1:iter

      if mat ( (i-1)*iter+j ) <= epsil, plot ( (i-1)/(iter-1) , (j-1)/(iter-1) , 'ro') ; end ;

      if mat ( (i-1)*iter+j ) > epsil, plot ( (i-1)/(iter-1) , (j-1)/(iter-1) , 'c*') ; end ;

   end ;

end ;

hold off

Voici d’abord le graphe des transitions de phases théorique, recueilli dans l’article de monsieur Bagnoli intitulé On Damage-Spreading Transitions que nous avions à étudier.

 

Voici  maintenant les deux graphes obtenus expérimentalement, le premier sur la densité asymptotique m et le second sur la distance de Hamming H.

Vous remarquerez que la zone bleue du diagramme de H contient bien la zone définie par la transition de phases de H sur le graphe théorique.



VI . Explication du diagramme de phases

de la densité et de la propagation de défauts ,

comparaison avec le Ferromagnétisme

Rappels :

p0 = P ( f(0,0) = 1) donc f(0,0)=0 constamment.

p1 = P ( f(0,1) = 1) = p2 = P ( f(1,0) = 1), cas symétrique.

P3 = P ( f(1,1) = 1).

On appelle réseau ordonné un réseau dont tous les sites ont pour valeur 0.

a)Explication du diagramme de  phase pour la densité

Lorsque m=0, cela signifie qu’il existe au plus un nombre fini de sites non nuls ; cela équivaut donc à dire que le réseau est presque ordonné après un temps très long.

On remarque sur ce diagramme que p1 joue le rôle le plus important dans la détermination de m. Par exemple si p1 est inférieur ou égal à 0.5, c’est-à-dire que la probabilité d’avoir  f (0,1)=0 est grande, alors la valeur de p3 n’intervient pas sur le résultat, on aura toujours m=0.

Ceci s’explique intuitivement par le fait qu’on a plus de chance de tirer un nombre aléatoire appartenant à [p1 ; 1] qu’à [0 ; p1]. On aura donc plus souvent  f (0,1)=0. Finalement  f (1,1) ne crée plus assez de 1. Pas après pas, « les 0 gagnent du terrain » et la densité asymptotique est donc nulle, ce qui correspond au cas ordonné. On obtient le même phénomène si p1 est strictement supérieur à 0.8. En effet quelle que soit la valeur de p3 on aura plus de chance de tirer un nombre aléatoire sur [0 ; p1] qu’entre [p1 ; 1], donc  f (0,1) va souvent donner 1. Finalement 1>m>0 dans cette région et on obtient un réseau désordonné.

Par contre lorsque p1 appartient à [0.5 ; 0.8], p3 détermine aussi l’état final. D’ailleurs si p3 reste relativement grand, les 0 ne l’emportent pas et donc 1>m>0, on est dans le cas désordonné. Lorsque p3 =1 et p1 =0.5, on est en M, le point critique du diagramme de phase. C’est en fait à partir de ce point que m ne vaut plus constamment 0 et que sa valeur dépend aussi de p3. C’est le premier point de la surface de transition rencontré lorsqu’on augmente p1.

Exemple : point où m=1 et point où m=0.

On se place dans des cas deterministes

·      Explication du point (1,1) :dans ce cas on a   f(0,0)=0

                                                                                   f(0,1)=1

                                                                                              f(1,1)=1

quand t tend vers l’infini s(i)=1 pour tout i donc m=1 quelle que soit la condition initiale (car on a en plus H=0).

·      Explication du point (0,0) : dans ce cas on a     f(0,1)=0

                                                                                     f(1,1)=0

                                                                                     f(0,0)=0

Donc quand t tend vers l’infini, s(i)=0 pour tout i donc m=0 quelle que soit la condition initiale.

b) Explication du diagramme de phase du au defaut

            Ce diagramme met en évidence deux régions du plan séparées suivant la valeur de la distance de Hamming H. Cette distance donne une idée de la différence après un temps « infini » entre deux réseaux suivant le même modèle mais ayant des conditions initiales différentes.

Lorsque H=0, cela signifie qu’il existe au plus un nombre fini de distances non nulles entre les deux répliques étudiées s et h (cf §I). Ceci est donc équivalent à dire que le système n’a pas la mémoire des conditions initiales. En effet nous avions deux réseaux avec des conditions initiales différentes et finalement après un temps infini, les deux réseaux sont égaux partout sauf en un nombre fini de points.

Par exemple, lorsque p1£0.8, quel que soit le couple (p1 ,p3), on obtient H nul.

Il existe seulement une toute petite région pour laquelle H est strictement positive, c’est-à-dire pour laquelle l’état final du réseau dépend de l’état initial du réseau. En effet si p1<0.3, la probabilité pour que f(1,1)=0 est élevée et si p1>0.8 la probabilité pour que f(0,1)=1 est élevée. Ceci explique que dans cette région  il n’y a pas convergence (le réseau a du mal à se stabiliser et sa configuration dépend des conditions initiales). Il faut donc un temps plus grand pour atteindre l’équilibre des deux phases et obtenir la transition de phases pour H.

c) Récapitulation des trois zones du diagramme de transition de phases

d) Extension et interprétation en termes de ferromagnetisme

            Les matériaux ferromagnétiques sont des substances qui, comme le fer, acquièrent ou possèdent des aimantations importantes. Le ferromagnétisme correspond à un état ordonné où tous les moments magnétiques élémentaires s’orientent parallèlement les uns avec les autres. En fait un matériau ferromagnétique change d’état selon la température. Si T>Tf, le matériau devient paramagnétique : il est dans une phase desordonnée (Tf : température de Curie). Par contre si T<Tf, on est dans l’état ferromagnétique ordonné. On peut interpréter les zones du diagramme des phases pour la densité comme étant des zones où le matériau est ferromagnétique ou paramagnétique.

On modélise les moments magnétiques par des sites valant 0 ou 1 selon la règle suivante :

Polarisation ­

Polarisation ¯

Site ayant la valeur 0

Site ayant la valeur 1

          Dans ces deux cas, on a un système stable, donc le matériau est dans une phase ferromagnétique.


Conclusion

Nous avons donc étudier expérimentalement le comportement asymptotique des automates cellulaires stochastiques, en se limitant à un cas très précis de fonction d’évolution, mais déjà une application pratique dans le domaine du ferromagnétisme a été possible.

Finalement, ce « jeu de la vie », dont le principe à l’air si simple, où certains ne voient que le moyen de faire de jolis économiseurs d’écran, a de multiples applications, les unes plus euristiques que les autres,  dans la modélisation de phénomènes physique, thermodynamique ou en biologie cellulaire.

Bibliographie

On Damage-Spreading Transition, Franco Bagnoli.