La suite logistique

Histoire

La suite logistique est un modèle mathématique qui est utilisé pour simuler l'évolution d'une population animale.
C'est une équation simple qui ressemble à ça:

		x(t+1) = r * x(t) * (1 - x(t))

Où:

x(t) est la population courante qui est comprise entre 0 et 1.
x(t+1) est la population à l'intervalle de temps suivant. Disons l'année suivante.
r est un "taux de fertilité". Sa valeur doit rester entre 0 et 4 pour garder la population entre 0 et 1.

La première partie de cette équation "r * x(t)" est facile à comprendre: La population courante multipliée par le
taux de fertilité nous donne la nouvelle population.
La deuxième partie "(1 - x(t))" est un peu plus difficile a expliquer. Elle représente un "taux de décès" à cause
de la limitation de la nourriture et des interactions entre les animaux.

Même si cette équation parait simple, son évolution est assez complexe quand r varie.
Et elle présente un comportement chaotique pour les valeurs de r les plus élevées.

Elle ne sert pas seulement à modéliser la population, elle est aussi utilisée en économie ou dans les réseaux de
neurones.
Stanislaw Ulam et John von Neumann voulaient l'utiliser comme générateur de nombres aléatoires.
Et cette équation est aussi liée à l'ensemble de Mandelbrot comme vous pouvez le voir dans une des vidéos de la
section "Liens".

Comportements varés

On va écrire un simple programme pour afficher l'évolution d'une population initiale de 0.5 avec une valeur donnée
de r.
Chaque valeur que l'on calculera avec l'équation sera séparée de 10 pixels le long de l'axe x.

		#define WIDTH   10

On définit notre valeur de "r":

		float r = 0.95f;

Puis dans la boucle principale on efface l'écran et on initialise x à 0.5

		while (sys.isQuitRequested() == false)
		{
			gfx.clearScreen(Color(0, 0, 0, SDL_ALPHA_OPAQUE));
			float x = 0.5f;

On boucle sur la largeur de l'écran 10 pixels par 10 pixels:

			for (int t = 0; t < SCREEN_WIDTH - 1; t += WIDTH)
			{

On stocke la valeur courante de x et on calcule la nouvelle valeur avec l'équation

				float oldX = x;
				x = r * x * (1 - x);

On multiplie ces valeurs par la hauteur de l'écran pour avoir une coordonnée y de pixel.

				int oldY = (1 - oldX) * (SCREEN_HEIGHT - 1);
				int y    = (1 - x)    * (SCREEN_HEIGHT - 1);

Et finalement on dessine une ligne entre ces 2 points:

				gfx.line(t, oldY, t + WIDTH, y, Color(255, 255, 255));
			}

		Télécharger le code source
		Télécharger l'exécutable pour Windows

Maintenant voyons ce qui se passe avec différentes valeurs du taux de fertilité.
Quand r est entre 0 et 1, la population meurt rapidement.
Voici le résultat pour r = 0.95:

Entre 1 et 2 la population se stabilise à une valeur en-dessous de celle de départ.
Ici pour r =1.3:

Entre 2 et 3.1 la population oscille un peu puis se stabilise à une valeur supérieure à celle de départ.
Pour r = 2.7:

Entre 3 et 3.44 la population oscille entre 2 valeurs.
Pour r = 3.2:

Entre 3.44 et 3.54 la population oscille entre 4 valeurs.
Pour r = 3.5:

Au-dessus de 3.54 la courbe semble de plus en plus aléatoire
Pour r = 3.8:

On peut apporter une petite modification au code pour animer la valeur de r, et comprendre plus facilement ce qui
se passe.
Au départ, on initialise r à 0:

		float r = 0.0f;

Puis à la fin de la boucle principale on va ajouter quelques lignes pour la faire augmenter lentement et la remettre
à 0 si on dépasse 4:

		sys.wait(16);
		r += 0.002f;

		if (r > 4.0f)
			r = 0.0f;

		Télécharger le code source
		Télécharger l'exécutable pour Windows

Vous pouvez voir le résultat sur ma chaine youtube. Voir le lien à la fin de cet article.

Le diagramme de bifurcation

Une autre façon de visualiser les points stables de ce système est de dessiner un diagramme où r varie le long de
l'axe x.
Pour chacune de ces valeurs de r on va dessiner un point là où se trouve la population après N itérations.
Où N est une grande valeur comme 10000.
Mais comme le système peut osciller on va aussi dessiner un point après N+1 itération, après N+2 itérations, etc.
Disons qu'on va afficher 150 points par valeur de r.

Alors au début du code on va écrire une fonction qui itère la formule n fois:

		float process(float r, int n)
		{
			float x = 0.5f;

			for (int i = 0; i < n; i++)
				x = r * x * (1 - x);

			return x;
		}

Puis dans la fonction principale, après avoir ouvert la fenêtre, on initialise les variables.
"r" va aller de 2.3 à 4.0.

		int t = 0;
		float startR = 2.3f;
		float endR   = 4.0f;

Dans la boucle principale on va dessiner 1 colonne par tour de boucle. "t" sera la coordonnée x de cette colonne
Et on commence notre boucle de 150 points.

		while (sys.isQuitRequested() == false)
		{
			if (t < SCREEN_WIDTH)
			{
				for (int i = 0; i < 150; i++)
				{

On convertit la coordonnée x courante en une valeur de r.
Et on appelle notre fonction process pour 10000+i itérations.

					float r = startR + (endR - startR) * t / (SCREEN_WIDTH - 1);
					float x = process(r, 10000 + i);

Comme avant on convertir le résultat de la formule en une coordonnée y et on affiche un point à cet endroit.

				float y = (1 - x) * (SCREEN_HEIGHT - 1);
				gfx.setPixel(t, y, Color(255, 255, 255));
			}

			t++;
		}

		Télécharger le code source
		Télécharger l'exécutable pour Windows

Et voilà le résultat

Si vous suivez ce diagramme de gauche à droite vous pouvez retrouver les observations qu'on a fait avant.

Pour les petites valeurs de r il n'y a qu'un point stable.
Puis la ligne se sépare en 2 branches, quand la population oscille entre 2 valeurs
Puis les 2 branches se divisent en 4, etc. Et très rapidement ça devient chaotique.

Une chose intéressante c'est qu'on peut voir que la partie droite n'est pas complètement chaotique. Il y a quelques
îlots de stabilité qui peuvent être expliqués par la vidéo sur l'ensemble de Mandelbrot.

Cette image ressemble à celle que l'on voit en bas de la page wikipedia anglaise même si ce n'est pas le même
modèle.
Mais au-dessus il y a d'autres images qui ont l'air plus belles où on peut voir plus de détails dans la zone
chaotique.

Le beau diagramme

Quand on observe ces images on voit que les pixels ne sont pas totalement blancs ou noirs, il y a différentes
teintes de gris.
En fait, ces images montrent plusieurs courbes superposées. On commence avec un fond vide et à chaque fois qu'on
dessine un point, on augmente un peu les valeurs R, V et B du pixel.

Alors voilà ce qu'on va changer dans le code.
Dans la boucle principale maintenant on va dessiner 2000 points par colonne.

		for (int i = 0; i < 2000; i++)
		{

Puis, comme avant, on calcule les coordonnées. Mais on doit bidouiller un peu la valeur de r pour éviter des trous
dans les lignes.

			float r = startR + (endR - startR) * (t + (i % 10) / 10.0f) / (SCREEN_WIDTH - 1);
			float x = process(r, 10000 + i);
			float y = (1 - x) * (SCREEN_HEIGHT - 1);

Maintenant, au lieu de dessiner un point, on va lire la valeur du pixel à ces coordonnées:

			Color col = gfx.getPixel(t, y);

Puis on augmente les valeurs des composants R, V et B du pixel, et on vérifie qu'elles ne dépassent pas 255:

			int c;

			c = col.r + 6;
			if (c > 255)
				c = 255;
			col.r = c;

			c = col.g + 6;
			if (c > 255)
				c = 255;
			col.g = c;

			c = col.b + 6;
			if (c > 255)
				c = 255;
			col.b = c;

Et finalement on réécrit le pixel à l'écran.

			gfx.setPixel(t, y, col);
		}

		Télécharger le code source
		Télécharger l'exécutable pour Windows

Et voilà l'image finale: