Exemples d'Applications

Plusieurs fonctions de la boite à outils « IsItChaos » peuvent être utilisées pour simuler des systèmes chaotiques. La boite à outils « modnum » présente également de telles fonctions Scilab en plus de plusieurs blocs Scicos. Nous allons présenter des exemples exploitant les deux boites.


Exemple IV.1: Le script suivant utilise la fonction « orb » pour simuler l'équation logistique: Xn=c..Xn-1.(1-Xn-1)

deff('[Xn]=f(Xn_1,c)','Xn=c*Xn_1*(1-Xn_1)');

x0=0.0001;

c=[3.5:.001:4] ;

orbits=[];

for i=c,

orbit=orb(20000,i,x0,f);

orbits=[orbits;orbit(1,$-999:$)];

end

Lest=[];

err=[];

save('Test_Logidt_K_orbits_c.dat',orbits,c);

for i=1:size(orbits,1),

[Nbr,L,Lesti,erri]=Lyap_K(orbits(i,:),,,,,1000,30, %T, 4);

Lesti=mean(Lesti);

erri=mean(erri);

if modulo(i,100)==0 then save('Test_Logidt_K_Lest_err.dat',Lest,err), end;

Lest=[Lest, Lesti];

err=[err, erri];

end;

plot2d(c,Lest);


20000 points de la trajectoire sont calculés pour chacune des valeurs de c allant de 3.5 à 4 avec une étape d'incrémentation de 0.001 et une condition initiale x0 égal à 0.0001. Seulement les 1000 derniers points de chaque trajectoire sont utilisés par la fonction « Lyap_K » pour estimer l'exposant de Lyapunov maximal par la méthode de Kantz.


À la fin un diagramme de la variation de l'exposant de Lyapunov en fonction c est généré. Le diagramme est présenté dans la figure IV-2 le long d'un diagramme des valeurs exactes de l'exposant de Lyapunov pris selon la référence [GOT02].


Exemple IV-2 : Le script suivant est identique au script de l'exemple 4-1 sauf qu'au lieu de calculer l'exposant de Lyapunov, la fonction « Chaos01 » est utilisée pour tester la nature chaotique de la série suivant la méthode 0-1 détaillée dans la section III.2.2.

deff('[Xn]=f(Xn_1,c)','Xn=c*Xn_1*(1-Xn_1)');

x0=0.0001;

c=[3.5:.001:4] ;

orbits=[];

for i=c,

orbit=orb(20000,i,x0,f);

orbits=[orbits;orbit(1,$-999:$)];

end

Kest=[];

C=[];

save('Test_Logidt_Chaos01_orbits_c.dat',orbits,c);

for i=1:size(orbits,1),

[Kesti,Ci]=Chaos01(orbits(i,:)');

C=[C,Ci];

disp(i)

disp(string(round(i/size(orbits,1)*100))+'%')

Kest=[Kest, median(Kesti)];

save('Test_Logidt_Chaos01_Kest_C.dat',Kest,c,C)

end;

plot2d(c,Kest);


À la fin un diagramme de la variation du paramètre K en fonction de c est généré. Le diagramme est également présenté dans la figure IV-2 le long des valeurs exactes de l'exposant de Lyapunov prises selon la référence [GOT02] et le résultat de l'exemple 4-1.





                                                                                                                     




                                                                     (a)


                                                                                                                 






                                                                    (b)





                                                                                                                 



                                                                     (c)



             Figure IV-2 Comparaison entre (b) l'estimation de l'exposant de Lyapunov par la méthode de KANTZ, (c) L'évaluation du paramètre K par la méthode 0-1

et (a) les valeurs exactes de l'exposant de Lyapunov

On constate que la fonction « Lyap_K » donne des résultats meilleurs que celles de « Chaos01 », mais « Chaos01 » garde l'avantage de ne pas avoir besoin d'une reconstitution de l'espace de phase par les variables à retard. Par contre, nous n'avons pas senti l'impact de cet avantage sur le temps de calcul.

En effet la première implémentation de la fonction « Chaos01 » demandait un temps de calcul excessivement long (de l'ordre de dizaines de jours) par rapport au temps demandé par la méthode « Lyap_K » (plusieurs heurs).

En premier lieu, nous avons attribué ceci au fait que « Lyap_K » interface une commande TISEAN et bénéfice ainsi d'une exécution machine directe, alors que « Chaos01 » passe par l'interpréteur intégré de Scilab. Ainsi, une deuxième implémentation de la fonction « Chaos01 » a été réalisée en s'assurant l'utilisation d'exécution machine directe dans la partie essentielle de la fonction. Malheureusement, l'écart de la durée d'exécution est resté considérable.

Par la suite, nous avons réalisé une inspection de l'algorithme de la méthode 0-1 et nous sommes venu à la conclusion que l'algorithme a une complexité exponentielle. Ce qui veut dire qu'il faut faire très attention au choix des paramètres de calcul de la méthode, et qu'une petite augmentation dans les paramètres peut induire un temps de calcul supplémentaire démesuré.

Ainsi, nous avons choisi des valeurs pour les paramètres de calcul différents des valeurs suggérées par les auteurs de la méthode, mais assez proches pour avoir des résultats qualitativement similaires.


Exemple IV-3: L'instruction suivante démontre l'efficacité gagnée par le mélange de Scilab et TISEAN. Dans cette instruction la fonction « Bifur » est utilisée pour évaluer la bifurcation du système de Henon simulé par la fonction « Henon_orbit » (voir Annexe A).

orbits=Bifur(Henon_orbit,,,0,2,.001);

La bifurcation est étudiée pour les valeurs du paramètre A allant de 0 à 2 et incrémentées par une étape de 0.001. Ce que nous ne pouvons pas montrer dans le manuscrit, c'est qu'en plus de la génération de la Figure IV-3 représentant la bifurcation, la fonction affiche à la fin de chaque itération la trajectoire finale correspondante à la valeur itérée de A, réalisant ainsi une animation du processus de bifurcation.














                                Figure IV-3 La bifurcation du système de Henon calculée dans l'exemple 4-3


Exemple IV-4: Le diagramme de simulation « Lorentz » de « modnum » v3.0 peut être utilisé pour simuler numériquement le fameux système (sys.I.2):

Dans cet exemple les paramètres sont choisis comme suite:

a=10 b=28 c=8/3 ci=[x0;y0;z0]=[5.5;5;20]

Période échantillonnage Tsampl=3e-3 s

Temps de la fin de simulation Tfin=20 s

Le diagramme Scicos est montré dans la figure suivante:












    



                F
igure IV-4: Bloc de simulation « Lorentz » de la boite « modnum » v3.0


On ouvre le bloc sous Scicos et on choisit la commande « run » du menu « simulate ». À la fin de la simulation, deux figures sont crées, la première (Figure IV-5 (a)) représente le plan de phase x contre z, et la deuxième (Figure IV-5 (b)) représente les variations temporelles de x,y et z.











                                                                                                                                                                                                                                             Temps (itérations)

                                           (a)                                                                                                    (b)

Figure IV-5 Résultats de la simulation du bloc « Lorentz »,
(a) plan de phase de x contre z, (b) évolution temporelle de x,y et z

En ajoutant un bloc « Mux » et un bloc « write ouput to file » on peut récupérer les valeurs de x,y et z dans un fichier pour pouvoir les exploiter par la suite. Nous allons supposer que le fichier s'appelle « Sim_Lorentz.dat » et nous allons l'utiliser dans les exemples suivants.

Exemple IV-5: Le script suivant lit le fichier « Sim_Lorentz.dat », récupère les valeurs de x et y, et dessine le plan de phase de x contre y. ensuite, le script dessine la reconstitution de l'espace de phase par les variables à retard. La reconstitution utilise la fonction « Delay_var » de « IsItChaos ». L'opération est répétée pour plusieurs valeurs de la dimension M et le délai L.

tmp_x_y_z=read('Sim_Lorentz.dat',-1,4);

tmp=tmp_x_y_z(:,1);

x=tmp_x_y_z(:,2);

y=tmp_x_y_z(:,3);


plot(x(20000:22500),y(20000:22500));

scf();


M=[2 3 4]

L=[5 25 60];

i=0;

for Mi=M;

for Li=L;

x1_x2_x3=Delay_var(x(20000:22500),2500,Mi,Li);

x1=x1_x2_x3(:,1);

x2=x1_x2_x3(:,2);

i=i+1;

subplot(3,3,i);

plot(x1,x2);

xtitle('M='+string(Mi)+' L='+string(Li),' x1','x2');

end;

end;

Les diagrammes résultants sont présentés dans les figures IV-6 et IV-7. On remarque qu'une valeur basse de L dans l'espace de phase, produit une représentation presque linéaire entre x et y. C'est à cause de la forte corrélation entre les points proches dans la série x. sur l'autre coté, une grande valeur de L produit une grande distorsion qui est due à la corrélation entre les points éloignés qui se croisent dans l'attracteur original.

La dimension M par contre n'a pas trop d'influence sur les résultats qui suggèrent qu'une dimension M=2 est suffisante.








        Figure IV-6: Plan de phase de x contre y de la simulation « Lorentz »

















Figure IV-7: Différentes représentations dans l'espace de phase reconstruit en utilisant plusieurs dimensions M et plusieurs retards L.


L'exemple suivant est similaire au précédent, mais au lieu d'étudier l'espace de phase en deux dimensions nous allons utiliser 3 dimensions.

Exemple IV-5: Le script suivant idem au précédent, récupère les valeurs de x, y et z, et dessine le plan de phase de x contre y et z. Ensuite, le script dessine une reconstitution de l'espace de phase par les variables à retard pour plusieurs valeurs de la dimension M et le délai L.

tmp_x_y_z=read('Sim_Lorentz.dat',-1,4);

tmp=tmp_x_y_z(:,1);

x=tmp_x_y_z(:,2);

y=tmp_x_y_z(:,3);

z=tmp_x_y_z(:,4);

plot(x(20000:22500),y(20000:22500));

scf();


M=[3 4 5 6]

L=[25];

i=0;

for Mi=M;

for Li=L;

x1_x2_x3=Delay_var(x(20000:30000),10000,Mi,Li);

x1=x1_x2_x3(:,1);

x2=x1_x2_x3(:,2);

x3=x1_x2_x3(:,3);

i=i+1;

subplot(2,2,i);

plot3d(x1,x2,x3);

xtitle('M='+string(Mi)+' L='+string(Li),'x1','x2','x3');

end;

end;

Les diagrammes résultants sont présentés dans les figures IV-8 et IV-9. On remarque que peu de différence existe dans les reconstructions de l'espace de phase à partir de la dimension 4. Ce qui est en accord avec les théorèmes présentés dans la section III-1.1.

Il faut noter que l'attracteur reconstruit n'a pas toujours la même forme visible que l'attracteur original, mais il a les mêmes caractéristiques topologiques estimées par les mesures invariantes.












                        Figure IV-8: Plan de phase du système de Lorentz simulé.













                Figure IV-9: Reconstruction de l'espace de phase pour différentes valeurs de M

Il serait intéressant de pouvoir estimer rapidement, et sans inspection visuelle, la dimension M suffisante et le délai L optimal. L'exemple suivant utilise la fonction « false_nearest » de la boite « IsItChas » pour résoudre le problème du choix de M.

Exemple VI-6: Le script suivant lit les valeurs de la variable x de la simulation « Lorentz » et les utilise dans la fonction « false_nearest » pour calculer le pourcentage des faux voisins en fonction de la dimension M (voir Section III.1.6). Le calcul se fait pour différentes valeurs de L. A la fin les diagrammes représentant les résultats sont dessinés.

tmp_x_y_z=read('Sim_Lorentz.dat',-1,4);

tmp=tmp_x_y_z(:,1);

x=tmp_x_y_z(:,2);

i=0;

for Li=[1 10 25 50];

[dim,FalseNberFrac,NberSize,NberSize2]=false_nearest(x(30000:35000),5000,2,Li,,8);

i=i+1;

subplot(2,2,i);

plot(dim,FalseNberFrac);

xtitle('L='+string(Li),'Dimension M', 'Pourcentage de faux voisins');

end;












                        Figure IV-10: Pourcentage de faux voisins en fonction de la dimension M.

On remarque que la pente de décroissance commence à se stabiliser à partir de M=4, ce qui est en accord avec les résultats de l'exemple 4-5. L'inspection visuelle ici n'est pas nécessaire et peut facilement être remplacée par un algorithme d'estimation assez simple.

Une méthode complémentaire pour trouver le délai L à utiliser sans inspection visuelle est suggérée dans [KAN02]. La méthode est originalement proposée par Fraser et Swinney et basée sur le calcul de l'information mutuelle retardée et prend en compte les corrélations non linéaires. Elle est implémentée dans la fonction « mutual » utilisée dans l'exemple suivant.


Exemple IV-7: Le script suivant utilise les valeurs de x résultantes de la simulation « Lorentz » pour estimer un délai L approprié pour la reconstruction de l'espace de phase.

tmp_x_y_z=read('Sim_Lorentz.dat',-1,4);

tmp=tmp_x_y_z(:,1);

x=tmp_x_y_z(:,2);


[Del,MutulaInf,OccupBox,ShanEnt]=mutual(x,10000,200,50);

plot(Del,MutulaInf);

xtitle('','Délais L', 'Information mutuelle ');

La figure IV-11 montre le changement de l'information mutuelle en fonction du délai. Les auteurs de la méthode proposent d'utiliser la valeur du délai qui correspond au premier minimum local comme valeur de L.

On peut voir que le premier minimum local se situe autour de 50, ce qui est relativement pas loin de la valeur que nous avons choisi (25), mais qui est quand même une valeur non optimale. Toutefois, cette valeur donne des résultats acceptables dans l'estimation des mesures invariantes.













 F
igure IV-11: Changement de l'information mutuelle en fonction du délai pour la variable x de la simulation « Lorentz ».


Dans tous les exemples précédents, le nombre d'échantillons de x utilisés pour le calcul a été donné directement. Mais comment choisir ce nombre? Et est-ce qu'il existe une limite optimale à ne pas dépasser? La méthode de calcul de la séparation spatio-temporelle proposée par Provenzale et al. peut être utilisée pour estimer une limite optimale de la taille de série. Cette méthode est utilisée avec les « graphes de récurrence » pour tester la non-stationnarité d'une série temporelle.

La méthode de « séparation spatio-temporelle » et celle des « récurrences » sont implémentées respectivement par les fonctions « stp » et « nstat_z » explorées dans l'exemple suivant.

Exemple IV-8: Le script suivant applique les fonctions « stp » et « nstat_z » aux valeurs de x extraites de la simulation « Lorentz », et dessine le résultat. Le dessin de la fonction « stp » est présenté dans la figure IV-12 et montre les variations des lignes de probabilité constante pour trouver un certain nombre de voisins dans l'espace de phase à une distance donnée (l'axe y) en fonction du délai qui sépare les échantillons. La figure IV-13 montre les graphes de récurrence traduits par l'erreur normalisée (axe y) résultante de l'utilisation d'un segment de la série d donnée pour estimer un autre segment.

mp_x_y_z=read('Sim_Lorentz.dat',-1,4);

tmp=tmp_x_y_z(:,1);

x=tmp_x_y_z(:,2);


strR=stp(x(30000:35000),5000,4,25,,500);

plot2d(stpR(:,1),strR(:,2),0*ones(strR(:,2)));

[forecast,forcasted,err]=nstat_z(x,40000,4,25,100);

scf();

plot2d(forecast,err,0*ones(err));

plot(ones(err)*mean(err),'r');





                                                                           

 





                                                                                            

Figure IV-12: La variation des lignes de probabilité constante du nombre de voisins à une distance donnée en fonction du délai entre les points.










            Figure IV-13: L'erreur normalisée de l'estimation croisée des segments dans la série x.


On peut interpréter la figure IV-12 comme une mesure de la corrélation spatio-temporelle, et cela montre la distance temporelle minimale entre les points pour qu'on puisse les considérer indépendants du point de vue des mesures invariantes. Ainsi il suffit de prendre une fenêtre de Theiler égale à 100 dans la recherche des voisinages. Le genre de minima observé pour autour de 240 est dû aux croisements résultants de repliement de l'attracteur chaotique sur lui même. Ce minima nous renseigne qu'il faut au moins 250 points pour une représentation minimale de l'attracteur, mais de préférence il faut prendre plusieurs multiples de ce nombre pour améliorer les estimations.

La figure IV-13 montre que la taille choisie des segments (40000/100=400) donne des erreurs normalisées de l'estimation-croisée d'une moyenne proche de 1, c.-à-d. une erreur proche de la variance des données elles-mêmes. Nous concluons qu'un nombre d'échantillons égal à 400 est représentatif de la série et que la série complète peut être considéré stationnaire.

L'exemple suivant explore l'utilisation de la fonction « poincare » pour réaliser des sections de Poincaré dans l'espace de phase reconstruit.

Exemple IV-9: Le script suivant réalise des sections de Poincaré dans l'espace de phase reconstruit par la série x avec une dimension M égale à 4 et un délai L égale à 25. Les sections sont réalisées suivant l'hyperplan perpendiculaire au troisième axe de l'espace de phase et pour différentes valeurs de l'intersection du plan avec l'axe. La section qui donne le plus d'intersections est sélectionnée pour le dessin. La figure IV-14 montre les lois de variations des points d'intersection avec la section de Poincaré selon les trois axes de l'hyperplan.

mp_x_y_z=read('Sim_Lorentz.dat',-1,4);

tmp=tmp_x_y_z(:,1);

x=tmp_x_y_z(:,2);


Maxsize=0;

Maxz0=-10;

for z0=[-10:0.5:10];

poinc=poincare(x,40000,4,25,0,3,z0);

if size(poinc)>Maxsize then Maxzo=z0; end;

end;

poinc=poincare(x,40000,4,25,0,3,Maxz0);

subplot(3,1,1);

plot2d(poinc(1:$-1,1),poinc(2:$,1),-1*ones(poinc(2:$,1)));

xtitle('Evolution du MAP de la section de Poincaré suivant l''axe x1',' x1(n-1)','x1(n)');

subplot(3,1,2);

plot2d(poinc(1:$-1,2),poinc(2:$,2),-1*ones(poinc(2:$,2)));

xtitle('Evolution du MAP de la section de Poincaré suivant l''axe x2',' x2(n-1)','x2(n)');

subplot(3,1,3);

plot2d(poinc(1:$-1,3),poinc(2:$,3),-1*ones(poinc(2:$,3)));

xtitle('Evolution du MAP de la section de Poincaré suivant l''axe x3',' x3(n-1)','x3(n)');

On peut constater sur les dessins des traces du caractère déterministe de la série en inspectant chaque figure et également en comparant les figures entre elles.
























                            Figure IV-14 Projection de la section de Poincaré de la série x suivant trois axes


Exemple IV-10 Le script suivant calcule l'exposant de Lyapunov maximal pour la série x par la méthode de Kantz. La fonction Lyap_K est utilisée pour différentes dimensions M et en utilisant une fenêtre glissante de la taille de 1000 échantillons. L'histogramme des valeurs estimées de l'exposant maximal de Lyapunov pour chaque fenêtre est dessiné et la valeur moyenne estimée.


tmp_x_y_z=read('Sim_Lorentz.dat',-1,4);

tmp=tmp_x_y_z(:,1);

x=tmp_x_y_z(:,2);

retard=1000;

fenetre=1000;

DimMin=4;

DimMax=6;

Delai=25;

Lest=[];

err=[];

NbrWin=floor(size(x,1)/retard);

for ii=0:NbrWin-1,

if ii*retard+fenetre<=size(x,1)...

then

orbit=x(ii*retard+1:ii*retard+fenetre);

else

break;

end;

[Nbr,L,Lesti,erri]=Lyap_K(orbit,,DimMin,Delai,,,30, %T,DimMax,500);

Lesti=mean(Lesti);

erri=mean(erri);

if modulo(ii,100)==0 then save('Test_Lyap_K__chaosimage_Lest_err.dat',Lest,err), end;

Lest=[Lest, Lesti];

err=[err, erri];

end;

Sur la figure IV-15, si on néglige les valeurs isolées, on peut voir que les valeurs de l'exposant tournent autour de l'intervalle [0.001,00.4] avec une moyenne de 0.0034433 ce qui correspond après ajustement par la division sur la période échantillonnage (0.003) du système « Lorentz » à 1.14.
















                          Valeurs de l'exposant de Lyapunov    
         
                                    

    Figure IV-15: Histogramme des valeurs de l'exposant de Lyapunov pour la série x.


L'exemple suivant montre l'utilisation de la fonction « dvv » pour l'application de la méthode DVV introduite dans la section III.3.2.


Exemple VI-11: Le script suivant explore la série x en utilisant la méthode DVV. En premier lieu, la fonction « dvv » est appliquée à la série elle-même, ensuite elle est appliquée à une série de substitution « xsur » qu'est générée en appelant la fonction « surrogates ».

tmp_x_y_z=read('Sim_Lorentz.dat',-1,4);

tmp=tmp_x_y_z(:,1);

x=tmp_x_y_z(:,2);


[DVV,rd]=dvv(x,1000,3,25,3,3);

xsur=surrogates(x(20000:30000),2,1);

[DVVs,rds]=dvv(xsur,1000,3,25,3,3);


plot(1:size(rd,2),DVV);

xtitle('M='+string(M)+' L='+string(L),'Distance','variance du cible');

scf();

plot(1:size(rd,2),DVV);

plot(1:size(rds,2),DVVs,'r-');

xtitle('M='+string(M)+' L='+string(L),'Distance','variance du cible Origine=Bleu Surro=Rouge');

Les figures suivantes traduisent les résultats obtenus.










F
igure IV-15: La variation de l'indice DVV de la série x en fonction de la distance (standardisée en commençant à μd - nd.σd).









Figure IV-16 Résultats du test DVV pour la série x (en bleu) comparés à ceux de la série régénérée (en rouge)

On peut constater de la figure IV-15 que la variance DVV suit un parcourt similaire au parcourt associé aux systèmes chaotiques traités dans la section III.3.2. On peut aussi déduire de la figure IV-16 que la série x est loin d'être linéaire.

Le dernier exemple illustre l'utilisation de la fonction « d2 » pour estimer la dimension de corrélation D2


Exemple IV-12: Le script suivant applique la fonction « d2 » à la série x et dessine les courbes résultantes.

tmp_x_y_z=read('Sim_Lorentz.dat',-1,4);

tmp=tmp_x_y_z(:,1);

x=tmp_x_y_z(:,2);


D2=d2(x(20000:30000),10000,12,50,,500);

plot2d(D2(:,1),D2(:,2),-1*ones(D2(:,2)))

xtitle('','Epsilon','Dimension');


La figure VI-18 montre le résultat de l'estimation. Malgré une déformation non attendue pour les petites distances (epsilon) la moyenne des résultats (2.176476) est proche de la dimension D2 actuelle de l'attracteur de Lorentz (2.06). Ces résultats restent discutables et nous classons la fonction « d2 » avec les fonctions sensibles exigeant l'interprétation d'un expert.













Figure VI-18 Résultats de l'application de la fonction « d2 » à la série x.


Comments