TangenteX.com

Les systèmes critiques auto-organisés auto-organisés

La découverte des systèmes critiques

Depuis longtemps, les physiciens se sont aperçu que les systèmes pouvaient être classés en plusieurs familles selon leur dynamique. Les systèmes linéaires, dont la réponse est proportionnelle à l'excitation, sont les plus simples. Bien souvent, le caractère linéaire de la réponse est une approximation : vous connaissez tous le pendule simple, qui est linéaire par approximation. Les systèmes non linéaires, qui sont très sensibles aux conditions initiales. Leur réponse à une excitation n'est plus linéaire mais dépend de l'excitation et des conditions initiales. Parmi les systèmes non linéaires, citons les systèmes chaotiques. Nous en avons croisé quelques uns sur TangenteX.

Et enfin, une classe de systèmes très particulière, et pourtant très répandue dans notre univers. Ce sont des systèmes formés de trés nombreux composants, qui accumulent doucement de l'énergie suite à des faibles excitations, puis qui relâchent plus ou moins d'énergie de façon plus ou moins périodique, suite à une excitation tout aussi faible.

Voyons l'exemple des tremblements de terre. Expérimentalement, on constate que l'énergie s'accumule doucement dans le système, formé par deux ou plusieurs plaques continentales qui frottent entre elles. Souvent, le système relâche son énergie de friction par des petits tremblements de terre, de faible énergie. Mais parfois, la relaxation (le relâchement d'énergie) provoque la libération d'une très grande quantité d'énergie, sans que l'on ne sache pas très bien pourquoi, sauf que c'est plus ou moins périodique. Les petits tremblements sont bien plus fréquents que les grands, heureusement.

Les modèles non linéaires classiques, comme le modèle "stick and slip" que nous avons examiné ici, ne décrivent pas bien ce comportement. Les physiciens ont donc abordé le problème sous un autre angle.

Ils ont constaté que le spectre de fréquences de ces évènements suivait une loi de puissance de type 1/f (à peu près). Et ils ont également constaté qu'il existe de nombreux systèmes physiques qui répondent à la même dynamique : circulation automobile, débit des grands fleuves, émissions radio des quasar, tâches solaires, bruit de scintillement (flicker noise) dans une résistance électrique, etc.

En 1987, le physicien danois Per Bak et ses collaborateurs Chao Tang et Kurt Wiesenfeld, ont proposé dans Self-Organized Criticality : an explanation of the 1/f noise - Phy. Rew. vol 59 number 4 juillet 1987, une théorie des systèmes critiques auto-organisés ou SOC (Self-Organized Criticality), qu'ils illustrèrent par l'exemple de la dynamique d'un tas de sable que l'on constituait en ajoutant grain par grain le sable sur un réseau. Cette dynamique est particulière : le système "tas de sable" est globalement dans un état stable (on dira plutôt métastable) alors qu'il présente de nombreuses instabilités locales.

Per Bak a identifié une autre particularité des SOC : leurs variables dynamiques, par exemple la taille des avalanches sur un tas de sable, présentent des invariants d'échelle. En d'autres mots, elles présentent des structures fractales.

Caractérisation d'un système critique auto-organisé

Dans son ouvrage "Self-organized criticality - Emergent complex behavior in physical and biological systems" (Cambridge University Press), le physicien anglais Henrik Jensen propose quatre critères caractéristiques d'un SOC ou plus généralement de la criticalité auto-organisée :

Per Bak et ses collaborateurs ont démontrés que leur modèle de tas de sable (BTW) présentait toutes les caractéristiques d'un SOC.

Modélisation d'un tas de sable

Le modèle BTW

Le modèle BTW est basé sur le concept d'automate cellulaire. Il ne tient pas compte des caractéristiques physiques du milieu granulaire que constitue un tas de sable réel. Cela a pour inconvénient de ne pas modéliser idéalement le comportement d'un tas de sable réel, mais cela a l'immense avantage de proposer un modèle pour la grande famille des systèmes critiques auto-organisés.

Pour adapter le modèle BTW, je me suis trés largement inspiré de l'ouvrage "How nature works - the science of self-organized criticality" de Per Bak, aux éditions Springer Science, dont je vous recommande trés vivement la lecture.

Son principe

Soit un tas de sable modélisé sur un réseau carré de coté N, sur lequel je définie une grandeur z(x,y,t) où x et y représentent les coordonnées d'un noeud (ou site) du réseau et t la variables discrète de temps, pour reprendre les notations de Peter Bak. Dans BTW, la grandeur z(x,y,t) représente le nombre de grains de sable empilés sur un site du réseau. Comme chaque grain est considéré comme un cube parfait d'arête 1, z(x,y,t) mesure la hauteur de l'accumulation de grains sur un site à l'instant t.

Au temps t=0, mon tas de sable ne contient aucun grain et donc pour tout site (x,y), z(x,y,t=0) = 0.

La dynamique de mon tas de sable est gouvernée par quatres règles très simples :

La règle 3 caractérise le comportement d'un système critique auto-organisé. Elle met en jeu un seuil qui provoque une réaction non-linéaire du système et elle permet une réponse du système qui peut être complétement disproportionnée au forçage du système. L'ajout d'un seul grain sur un site choisi aléatoirement peut provoquer l'avalanche de plusieurs dizaines ou centaines de grains de sable.

Quelle sera l'évolution de mon tas de sable, si je pars d'un réseau vide ? Au début de la simulation, le tas de sable va grossir progressivement et nous pourrons observer localement quelques avalanches mettant en jeu assez peu de grains de sable. Puis le tas va gagner en hauteur et en largeur et la fréquence et la taille des avalanches vont croître, jusqu'à ce que le système ait atteint sa valeur d'équilibre dynamique.

A ce stade, on s'apercevra que la distribution temporelle des avalanches présente un spectre particulier, caractéristique des SOC, un spectre en 1/f, qu'on désigne dans la littérature sous le nom de "bruit en 1/f".

Le tas de sable comme modèle de système critique auto-organisé

Plus généralement, notre modèle de tas de sable présente les caractéristiques suivantes :

D'après la définition de Per Bak, il s'agit donc d'un représentant de la famille des SOC. Le modèle BTW caractérise un système critique auto-organisé comme notre tas de sable, par deux propriétés:

Les observables du modèle

J'ai choisi de m'interesser à trois observables dans mon adaptation du modèle BTW (ils n'ont rien d'original !) :

A l'aide de ces trois observables, je vais essayer de caractériser le système constitué par mon tas de sable simulé. Il est bien sur possible de s'intéresser à d'autres observables, par exemple l'évolution de la pente moyenne du tas de sable ou encore la répartition entre sites critiques et non critiques en fonction du temps.

La simulation

Je vous propose maintenant d'implémenter mon adaptation du modèle BTW dans un script Python, dont je vais détailler les principales composantes ci-dessous. Le code complet est inclus dans le script TasSableSimu.py. Accessoirement, la simulation étant assez longue vu son nombre d'itérations (105 à 106 pour avoir des résultats probants), il serait opportun et assez facile d'écrire ce bout de code en C.

Boucle d'évolution du tas de sable

A l'état initial, le tas de sable ne contient aucun grain de sable. Il est créé par les instructions :

N = 40                  # taille du réseau carré
Tas = zeros([N,N],int)  # réseau pour la simulation

où N est la taille du réseau carré sur lequel je vais bâtir notre tas de sable. C'est une hypothèse parmi d'autres, dans certaines adaptations de BTW, le tas de sable initial est rempli aléatoirement de grains, sous son niveau critique.

Pour l'édifier, j'utilise une boucle avec un nombre d'itérations NbIter très élevé (105 dans ma simulation). Dans cette boucle (règle 4), je procède aux actions suivantes, en appliquant les règles que j'ai défini ci-dessus :

Ce qui me donne le code :

for k in xrange(NbIter):  
    # désignation d'un site aléatoire
    i,j = TirageSiteUniforme(N)
    # ajout de grain de sable
    Tas[i][j] += PasGrain
    # application des règles de simulation d'avalanche 
    Tas, NbAval,NbGr = DetectionAvalanche(Tas,N)
    NbAvalanches.append(NbAval)
    NbGrains.append(NbGr)
    # calcul et stockage du remplissage
    Remplissage.append(Tas.sum())

Traitement des avalanches

Le traitement des avalanches est assurée par la fonction DetectionAvalanche(), qui implémente les règles 2 et 3  :

def DetectionAvalanche(tas,n):
    nbav = 0       # nombre d'avalanches
    nbgrains = 0   # nombre de grains impliqués dans les avalanches
    indicav = True # indicateur présence avalanche
   
    # boucle de traitement des avalanches   
    while indicav :       
        # traitement des sites intérieurs
        for i in xrange(1,n-1):
            for j in xrange(1, n-1):                     
                # détection de pente critique
                if tas[i][j] - tas[i][j+1] > Zc or \
                   tas[i][j] - tas[i][j-1] > Zc or \
                   tas[i][j] - tas[i+1][j] > Zc or \
                   tas[i][j] - tas[i-1][j] > Zc :                     
                       # éboulement des grains sur le voisinage                     
                       tas[i][j] = tas[i][j] - 4   
                       tas[i+1][j] += 1
                       tas[i-1][j] += 1
                       tas[i][j+1] += 1
                       tas[i][j-1] += 1
                       nbgrains += 4
                       nbav += 1
                else:
                     indicav = False                
        # traitement des bords (conditions aux limites = 0)
        for i in xrange(n):
            tas[i][0] = 0
            tas[i][-1] = 0
        for j in xrange(n):
            tas[0][j] = 0
            tas[-1][j] = 0
    return tas,nbav,nbgrains

La boucle while tourne tant que le tas de sable n'a pas atteint son état d'équilibre, c'est à dire qu'il n'y ait plus d'avalanche sur le tas (l'indicateur indicav devient alors False).

Pour détecter ce retour à l'équilibre, une double boucle for parcourt l'intérieur du réseau (hors des bords) et vérifie pour chaque site, si la pente avec ses voisins ne dépasse pas la pente critique Zc Si c'est le cas, une avalanche est déclenchée. Les conditions aux limites sont appliquées aux bords gauche et droit puis bas et haut (règle 2). Sur les bords, les grains sont éjectés du tas (les sites des bords sont remis à 0).

la fonction DetectionAvalanche() retourne :

Ces données seront stockées respectivement dans la matrice Tas et dans les listes :

Stockage des données pour exploitation ultérieure

Cette partie du code, qui permet le stockage des données de simulation pour une exploitation ultérieure, met en oeuvre la fonction pickle() de Python. Elle permet la sauvegarde sous format binaire de n'importe quel objet Python, et donc de nos paramètres de manip et de nos listes de résultats :

fic = open('Data/TasSable.dat','w')
pickle.dump(N,fic)
pickle.dump(Zc,fic)
pickle.dump(PasGrain,fic)
pickle.dump(Tas,fic)
pickle.dump(NbAvalanches,fic)
pickle.dump(NbGrains,fic)
pickle.dump(Remplissage,fic)
fic.close()

Les paramètres et les données de simulation sont stockées dans le fichier TasSable.dat, que je sauvegarde sur le sous-répertoire Data de mon répertoire courant.

Note : si vous voulez écrire ce code en C, il vous faudra stocker les données dans un fichier texte, par exemple un fichier csv, qui pourra être exploité soit par le script TasSableAnalyse.py (après adaptation) soit par GnuPlot ou autre.

Récupérer les données de simulation pour les analyser

Dans le script TasSableSimu.py, j'ai sauvegardé les paramètres et les données de simulation pour traitement ultérieur. Pour analyser ces données, je vais donc devoir les récupérer dans le fichier TasSable.dat. Le code Python suivant procède à cette opération :

fic = open('Data/TasSable.dat','r')
N = pickle.load(fic)
Zc = pickle.load(fic)
PasGrain = pickle.load(fic)
Tas = pickle.load(fic)
NbAvalanches = pickle.load(fic)
NbGrains = pickle.load(fic)
Remplissage = pickle.load(fic)
fic.close()

Attention, l'ordre de récupération par la fonction pickle() doit être identique à celui de la sauvegarde.

Le script TasSableAnalyse.py contient la totalité du code utilisé dans la suite pour l'analyse de cette simulation.

Le tas de sable en 2D et 3D

Pour une grille de dimension N = 40 et les paramètres suivants :

NbIter = 100000
PasGrain = 1
Zc = 3

J'obtiens un tas de sable dont voici l'aspect en 2D :

Tas de sable - projection 2D

et en 3D :

Tas de sable - Vue 3D

Notre tas de sable théorique présente une forme pyramidale, de hauteur moyenne d'environ 25 grains.

Remplissage du tas de sable

Le tas de sable grossit de PasGrain (PasGrain = 1 dans la simulation) grains par itération. L'accroissement de la taille du tas de sable est approximativement linéaire jusqu'à une limite d'itérations (environ 31 000 itérations dans la simulation), limite au délà de laquelle l'accroissement du tas cesse et sa taille devient pratiquement constante. Cette taille est d'environ 14 600 grains dans la simulation. Le tas de sable est dans un état quasi-stationnaire : il évolue très peu et de manière qui semble aléatoire autour d'une valeur moyenne de remplissage.

Nous constatons ce phénomène sur la courbe de remplissage ci-dessous :

Tas de sable - Courbe de remplissage

Lors de la phase stationnaire de l'évolution du tas, disons à partir de l'itération 40000, le nombre de grains ajoutés au tas est compensé par le nombre de grains qui s'échappent du tas par les bords, suite aux avalanches.

Essayons maintenant de zoomer sur une petite portion de cette courbe, dans l'état quasi-stationnaire. Nous obtenons ceci :

Tas de sable - zoom Courbe de remplissage

Nous constatons que le nombre de grains croît régulièrement et de façon approximativement linéaire, puis que périodiquement (à peu près), il décroît brutalement. Cette décroissance brutale est due à une ou plusieurs avalanches qui ont provoqué l'éjection de grains du tas.

Autre chose, le tas de sable a une forme approximativement pyramidale. Comme vous le savez, le volume d'une pyramide de base N2 et de hauteur h est de \( V = \dfrac{N^2 h}{3}\), soit avec les valeurs de simulation N= 40 et h = 25, un volume d'environ 13 300 unités. C'est très proche du volume du tas à l'état quasi-stationnaire (14 600). La hauteur h est obtenue en sortant le max de la variable Tas.

Le nombre d'avalanches

Etude de la distribution du nombre d'avalanches

Pendant une itération, c'est à dire après l'ajout d'un grain de sable sur un site aléatoire du tas, nous assistons à un ajustement de l'organisation du système, une relaxation. Cet ajustement doit contribuer à diminuer l'énergie totale du tas. Cela se traduit par un certain nombre d'avalanches. Le nombre d'avalanches mesure l'énergie dissipée lors de chaque itération. Le nombre d'avalanches pendant une itération mesure aussi le temps de relaxation du système "tas de sable".

La liste NbAvalanches contient, pour chaque itération, le nombre d'avalanches survenues sur l'ensemble du tas de sable lors d'une itération. La courbe ci-dessous illustre l'évolution du nombre d'avalanches en fonction de l'itération :

Tas de sable - Nombre d'avalanches

Cette courbe montre un signal aléatoire, qui n'est pas facilement compréhensible en l'état. 

Pendant la phase de croissance du tas, disons entre les itérations 0 et 25 000, le nombre d'avalanches est assez faible. Il croît avec le volume du tas de sable, phénomène que l'on aperçoit bien entre les itérations 25 000 et 30 000.

Nous observons également que le nombre d'avalanches par itération peut varier de façon considérable pendant l'état quasi-stationnaire. Dans cet état, le volume du tas reste quasi-constant, mais il maintient son organisation au prix d'un nombre plus ou moins grand de restructurations que sont les avalanches.

La taille des avalanches

Etude de la taille des avalanches

La liste NbGrains nous donne, pour chaque itération, le nombre de grains impliqués dans les avalanches. La courbe ci-dessous illustre l'évolution de la taille des avalanches en fonction de l'itération :

Tas de sable - Taille des avalanches

Distribution de probabilité du taille des avalanches

Pour mieux discerner les motifs dans ce signal, j'ai calculé la distribution de probabilité de la taille des avalanches à l'état quasi-stationnaire avec ce code :

Classes, DistribAv = DistributionProba(NbGrains[T:])
fig7 = figure(7,figsize=(10,8))
grid(True)
xlabel(u"Taille de l'avalanche (en nombre de grains)")
ylabel(u'Probabilité')
title(u"Distribution de probabilité d'une avalanche en fonction de sa taille")
loglog(Classes,DistribAv)

Le code de la fonction DistribAvalanches() est classique lorsqu'il s'agit de calculer une distribution de probabilité :

def DistributionProba(Liste):
    densiteproba = []
    Max = max(Liste)+1
    n = float(len(Liste))
    classes = range(4,Max,4)
    for i in classes:
        nbelem = 0
        for j in Liste:
            if j == i:
                nbelem += 1
        # normalisation à 1 de l'histogramme -> densité de proba        
        densiteproba.append(nbelem/n)       
    return classes,densiteproba

J'ai découpé l'histogramme de la taille des avalanches par classes de largeur = 4 , qui correspond à la taille la plus petite d'une avalanche. Cette largeur est assez petite pour me permettre d'assimiler l'histogramme normalisé comme une approximation de la fonction de densité de probabilité.

Nous obtenons la distribution suivante :

Tas de sable - Distribution taille avalanches

La taille minimum d'une avalanche est de 4 grains, conformément à notre modèle. Nous constatons que les avalanches les plus petites sont aussi les plus probables. Les avalanches les plus grandes, qui traduisent le plus grand relâchement d'énergie lors d'une itération, sont les moins probables.

Le spectre des tailles d'avalanches

Pour tracer ce spectre, j'utilise le code suivant :

# calcul et affichage du spectre de la taille des avalanches en régime stationnaire
T = 50000                # itération de début du régime stationnaire
NbEch = 2048             # taille de l'échantillon pour FFT (sur NbEch itérations)
Ech = NbGrains[T:T+NbEch]
X = abs(fft(Ech))/NbEch
freq = fftfreq(len(X))
freq = freq[range(NbEch/2)]
X = X[range(NbEch/2)]
fig5 = figure(5,figsize=(10,8))
grid(True)
xlabel(u"Fréquence (Hz)")
ylabel(u'Amplitude')
title(u"Taille des avalanches - Spectre normalisé")
semilogy(freq,X)

Rien de bien nouveau : je calcule la transformée de Fourier sur un échantillon de 2048 points, prélevé dans l'observable NbGrains, en me limitant à la portion correspondante à l'état stationnaire de mon système.

Voici la courbe obtenue :

Tas de sable - Spectre taille avalanches

La densité spectrale de puissance

La taille des avalanches est une variable pseudo-aléatoire en fonction des itérations (du temps). Le calcul de la densité spectrale d'un tel signal est un peu particulière : il faut découper le signal en petites tranches temporelles, puis calculer le spectre de chacune de ces tranches. Puis enfin, on calcule la moyenne de ces spectres. C'est assez casse-pied, mais heureusement, comme c'est aussi assez commun, Python a tout prévu. La fonction matplotlib.psd nous aidera à faire ce calcul. Sa documentation de cette fonction est disponible ici.

Echantillon = NbAvalanches[T:]
fig9 = figure(9,figsize=(10,8))
title(u"Energie totale - Densité spectrale")
spectre, frequence = psd(Echantillon,NFFT=512,Fs=1.0,window=mlab.window_none)

Cette fonction mérite quelques explications. Le signal contenu dans Echantillon est découpé en segments de longueur NFFT. NFFT représente le nombre de points de la transformée de Fourier qui sera appliquée au segment. Eventuellement, il est possible d'appliquer une fenêtre sur le signal pour calculer son spectre (je ne l'ai pas fait). La fonction psd() calcule donc le spectre de chaque segment et fait la moyenne des spectres. Elle retourne un vecteur contenant la densité spectrale de puissance du signal et un vecteur contenant les fréquences en abscisse. Nous les utiliserons dans le paragraphe suivant. Et enfin la fonction psd() trace la courbe de densité spectrale :

Tas de sable - Densité spectrale de puissance

Cette courbe ressemble diablement à celle d'une fonction en 1/x, ce qui signifierait que la densité spectrale du signal diminue de façon inversement proportionnelle à sa fréquence. Autrement dit, plus une avalanche globale lors d'une itération est importante, la quantité d'énergie relaxée grande, et moins sa fréquence est grande.

Le bruit en 1/f

Nous avons supposé que notre modèle BTW du tas de sable décrivait un système de la famille des SOC. Sa dynamique devrait donc laisser apparaitre une loi de puissance, en l'occurence un bruit en 1/f. Essayons de vérifier cette hypothèse.

Caractérisation du bruit en 1/f

Pour caractériser le bruit en 1/f en faisant apparaitre la droite en coordonnées log-log, j'utilise la fonction matplotlib.loglog() en l'appliquant sur le vecteur de densité spectrale calculée ci-dessus :

fig10 = figure(10,figsize=(10,8))
grid(True)
xlabel(u"Fréquence (Hz)")
ylabel(u'Amplitude spectrale')
title(u"Caractérisation du bruit en 1/f")
loglog(frequence,spectre)

J'obtiens la courbe suivante :

Tas de sable - Bruit en 1/f

Nous observons un plateau pour les basses fréquences puis une droite de pente environ égale à -1. Ce schéma est typique d'un bruit de type \( f^{-\alpha} \). BTW a un peu rapidement fait l'hypothèse que \( \alpha = 1 \), ce qui a été contredit en partie depuis.

De la simulation à la réalité

Les expériences autour d'un tas de sable sont légions. Il est possible de les faire assez simplement et elles font d'ailleurs l'objet de bon nombre de TIPE.

Le modèle BTW et ceux qui l'ont suivis sont très instructifs et ont donné lieu au développement d'une théorie des systèmes critiques auto-organisés très approfondie.

Cependant, il faut bien constater que les tas de sable de la réalité ne se comportent pas vraiment selon la loi de puissance que nous venons d'identifier. Sans doute parce que cet algorithme ne tient pas compte de la physique complexe (!) des milieux granulaires.

Dès 1989, des physiciens ont essayé de monter des expériences qui permettraient de vérifier le modèle BTW. Plusieurs difficultés apparurent : procéder exactement au dépôt sur le site aléatoire choisi et ne déposer qu'un seul grais de sable à la fois ! Ces expériences ne furent pas vraiment un succés : ils ne constatèrent pas l'apparition de la loi de puissance dans la distribution de la taille des avalanches.

Le modèle BTW ne tient pas compte de phénomènes tels que la friction entre les grains qui entraine un effet d'hystérésis. Contrairement au modèle, on n'observe pas systématiquement une avalanche lorsque la pente critique est dépassée, ce qui contrarie notre modèle... D'autre part, les expérimentateurs constatèrent que les grains de sable avaient tendance à rouler, en dissipant ainsi de l'énergie (cinétique) autrement qu'en relaxant comme le prévoit le modèle BTW.

Finalement, en 1996, des physiciens norvégiens (Frette, Christensen et all.) ont réussi à s'approcher du modèle BTW en utilisant des grains de riz en lieu et place de grains de sable. Et pas n'importe quel type de grains de riz : du riz long. Ils ne dissipent pas leur énergie cinétique en roulant ! Avis à ceux qui voudraient vérifier expérimentalement le modèle BTW...

Il n'en reste pas moins vrai qu'il est stupéfiant de voir les études que l'on peut tirer d'une simulation presque triviale, qui tient en quelques lignes de code.

Contenu et design par Dominique Lefebvre - www.tangenteX.com mars 2019 Licence Creative Commons Contact :