En septembre 1991, la revue Scientific American présentait à ses lecteurs une collection de très belles images obtenues à l'aide d'un algorithme à peine plus compliqué que celui utilisé pour générer les ensembles de Mandelbrot ou de Julia. Ces images furent créées en utilisant une extension de la formule logistique

Pi+1 = kPi(1 - Pi)    (*)

Pour obtenir de telles images, on a utilisé une quantité appelée exposant de Lyapunov à la mémoire du mathématicien russe Alexander Markus-Lyapunov (1857-1918).

Cette quantité sert à mesurer le degré de sensibilité d'un système comme celui du haut, suite à une variation infinitésimale de sa condition initiale Po. Un système sensible à de très petites variations de la condition initiale Po (une des caractéristiques des systèmes chaotiques) aura une quantité positive. La quantité sera négative si de petites variations de Po n'ont aucun effet à long terme sur le système. L'exposant de Lyapunov est en fait une mesure du degré de stabilité d'un système.

Prenons le système associé à la formule logistique du haut (*) et Po = 1/2. Le graphique du bas met en relation pour cette formule, l'exposant de Lyapunov (en bleu) et le diagramme de bifurcation (en vert) pour des valeurs de k variant de 0 à 4.

 

On voit très bien que les valeurs négatives de l'exposant de Lyapunov correspondent aux valeurs k pour lesquelles le système se stabilise avec le temps. Plus la valeur de l'exposant de Lyapunov est grande et négative, plus le système sera stable. Inversement pour les valeurs positives, le système passe dans une phase chaotique et devient sensible à de petites variations de sa condition initiale Po.

Comment arrive-t-on à créer de belles images en utilisant l'exposant de Lyapunov?

Pour obtenir ces images, on doit utiliser deux taux de croissance k1 et k2 que l'on fait varier selon une loi périodique appelée racine. Par exemple la racine k1, k2, k2, k1 est associée au système

Pi+1 = k1Pi(1 - Pi)
Pi+1 = k2Pi(1 - Pi)
Pi+1 = k2Pi(1 - Pi)
Pi+1 = k1Pi(1 - Pi)

En faisant varier k1 et k2 entre deux valeurs, on crée une grille rectangulaire. Pour chaque point ( k1, k2) de la grille, on calcule l'exposant de Lyapunov puis on colore ce point d'une couleur associée à la valeur obtenue. Par exemple, on peut choisir de colorier le point en bleu pour des valeurs positives et en jaune pour des valeurs négatives ou colorier le point en choisissant une multitude de teintes reliées à la valeur de . Évidemment, on obtient de meilleurs résultats lorsque le nombre n d'itérations est élevé mais on constate rapidement qu'une image demande énormément de calculs et par conséquent énormément de temps. Toutes les images qui apparaissent dans cette page ont été obtenues à l'aide du logiciel Maple. J'ai utilisé par la suite le logiciel de traitement d'image Photoshop pour améliorer les teintes et les contrastes.


racine : k2,k2,k1,k1


racine : k2,k2,k1,k1,k1,k1,k2


racine : k2,k2,k2,k1,k1


racine : k2,k2,k2,k2,k2,k2,k2,k2,k2,k2,k1,k1,k2,k2,k1,k1

Certaines parties des images ont des contrastes plus marqués. Ces parties d'image sont associées à des exposants ayant de très grandes valeurs négatives. L'erreur initiale tend à disparaître dans ces régions. L'arrière-plan de l'image correspond aux valeurs positives de l'exposant pour lesquelles le système est sensible à d'infimes variations de sa valeur initiale. Cette région peut engendrer des zones chaotiques.

Voici les trois programmes que j'ai utilisés pour écrire cette page. Si vous avez le logiciel Maple, copiez les programmes (en rouge) sur une feuille de travail du logiciel et changez au besoin les paramètres indiqués. Le dernier programme est celui que j'ai construit pour créer les fractales de Lyapunov. Soyez patient si vous utilisez ce programme, il demande énormément de calculs.

    Programme I  (L'exposant de Lyapunov et le diagramme de bifurcation en fonction de k)

Lyapunov:=proc(Table,xo,a,b,n)
local x,i,j,k,t,init,iter,exposant;
init:=200;
iter:=400;
j:=0;
for t to n do;
k:=a+t*(b-a)/n;
x:=xo;
for i from 1 to init do;
x:=k*x*(1-x);
od;
exposant:=1;
for i from 1 to iter do;
x:=k*x*(1-x);
exposant:=exposant*(k-2*k*x);
od;
if exposant<>0 then exposant:=log(abs(exposant))/iter else exposant:=-infinity fi;
j:=j+1;Table[j,1]:=k; Table[j,2]:=exposant;
od;

end:

Bifurcation:= proc(Table,xo,a,b,n)
local x,i,j,k,t;
j:=0;
for t to n do;
k:=a+t*(b-a)/n;
x:=xo;
for i to 120 do;
x:=k*x*(1-x);
if i>100 then j:=j+1; Table[j,1]:=k; Table[j,2]:=x fi;
od;
od;
end:

xo:=0.5:
a:=0:
b:=4:
n:=1000:


Table:= array(1..n,1..2):
evalhf(Lyapunov(var(Table),xo,a,b,n)):
g1:=plots[pointplot](Table, symbol=POINT, connect=true, color=blue):

Table:= array(1..20*n,1..2):
evalhf(Bifurcation(var(Table),xo,a,b,n)):
g2:=plots[pointplot](Table, symbol=POINT, color=green,labels=["k",""]):

plots[display]({g1,g2});

Programme II  (Les fractales de Lyapunov)

couleur:=proc(a,b)
local x,i,init,iter,exposant;
x:=.5;
init:=200;
iter:=400;
for i from 1 to init do;
x:=a*x*(1-x);
x:=b*x*(1-x);
x:=b*x*(1-x);
x:=a*x*(1-x);

od;
exposant:=0;
for i from 1 to iter do;
x:=a*x*(1-x);
exposant:=exposant+log(abs(a-2*a*x));
x:=b*x*(1-x);
exposant:=exposant+log(abs(b-2*b*x));
x:=b*x*(1-x);
exposant:=exposant+log(abs(b-2*b*x));
x:=a*x*(1-x);
exposant:=exposant+log(abs(a-2*a*x));

od;
exposant:=exposant/iter:
if exposant<0 then exposant:=exposant*(-1);
else exposant:=0;
fi;
end:

plot3d(0,1..4,1..4,orientation=[-90,0],style=patchnogrid, axes=framed, numpoints=30000,color=couleur);


Fractale de Lyapunov obtenue à l'aide de Maple.

Fractale de Lyapunov après traitement dans Photoshop.
  • Le dernier programme utilise la racine k1, k2, k2, k1 (abba). On peut étudier d'autres racines en changeant les lignes vertes du programme.
  • Pour effectuer un ZOOM, modifiez l'étendue 1..4,1..4, mais faites attention aux décimales. Maple 7 gère très mal les étendues contenant des décimales dans certaines macro-commandes. Écrivez 21/10..34/10 plutôt que 2.1..3.4.
  • Pour améliorer la résolution de l'image, changez la valeur de numpoints. Mais attention si on augmente la résolution, le programme prendra beaucoup plus de temps. Pour obtenir chacune des images de cette page, il a fallu environ une heure à chaque fois en travaillant sur un Mac avec un processeur PowerPC G4 de 800 MHz.
  • Notez que la première boucle que l'on répète 200 fois (init:=200) est ignorée dans le calcul de . Elle est présente dans ce programme uniquement pour permettre au système de s'adapter.
  • Pour calculer , le programme effectue 400 itérations (iter:=400) .On peut améliorer la précision de l'image en augmentant la valeur donnée à iter. En revanche, il faudra attendre plus longtemps pour obtenir cette image.