Discussion :

[lopppez]

Bonsoir à toutes et tous
e cherche un code qui peut ploter l'histogramme ( la distribution) d'une suite
et aussi de comparer entre deux distribution de deux suites
le fichier en attachement
merci

[lopppez]

salut;
Lors de ma recherche je trouve ce code il plot la mesure invariante d'une suite récurente appelé logistique [x(n+1)=ax(n)(1-x(n))] une suite de dimension 1

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
%%%%%%% INVARIANT MEASURE FOR LOGISTIC MAP
clear all;
 
%%%%%%%%%%%% PARAMETER OF THE MAP 
a=4;
 
%%%%%%%%%%%% NUMBER OF BOXES; BOX SIZE
nx=128;        dx=1/nx; 
 
c=zeros(nx,1);
 
%%%%%%%%%%% INITIAL CONDITION
x=0.4;
 
%%%%%%%%%%% NUMBER OF ITERATIONS
N=240000;
 
%%%%%%%%%%% MAIN CYCLE
tic
for ic=1:N
 x = a*x*(1-x);
 ix=floor(x(1)/dx)+1; 
 c(ix)=c(ix)+1;
end
toc
 
%%%%%%%%%% NORMALIZATION 
c=c/N/dx;
 
 
figure(3); hold off;
plot((dx/2:dx:1-dx/2),c); axis('square');
 
%%%%%%%%%%%%% CHECK NORMALIZATION
sum(c*dx)

at aussi un autre code pour la timeseriie de cette suite car
une autre méthode je peut dresser l'estogramme (ou la mesure) de trouver la serie temporelle ensuite utiliser EXEL pour ploter l'histogramme

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
a=4          % Parameter value
x0=0.8         % Initial condition
N=100;        % Number of iterations
 
x(1) = x0;       % Matrices in matlab cannot have zero index
 
% compute the orbit and print out results
fprintf(1,'Year        Density\n');
for ic=1:N
 x(ic+1) = a*x(ic)*(1-x(ic));
 fprintf(1,'%d     %f\n',ic,x(ic+1));
end
 
% graph the orbit
figure(1); hold off;
plot(x,'r*'); 
hold on; 
plot(x,'y');
axis([1 N 0 1]);

MOI je veut faire la meme chose mais en dimension deux avec par exemple
cette suite

H: (x, y) -> (y, 1-ay2+bx)

mais j'arive pas
quelqu'un peut m'aider

[lopppez]

et pour la suite de 2 dimension je trouve que le code qui plot l'attracteur et la serie temporelle

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
% Iterating the Henon map
 
clear
 
henon_map = inline('[a-x.^2+b*y x]','a','b','x','y');
 
a=input('a [1.4] = ');
b=input('b [0.3] = ');
 
if length(a)==0
 a=1.4
end
 
if length(b)==0
 b=0.3
end
 
%x0=0.1;y0=-0.5;
x0=-0.001;y0=0.05;
 
Ntransient=1000;
Nshow=10;
Npts=5000;
Ntotal=Ntransient+Npts;
fs=15;
 
X=x0;
Y=y0;
xo=x0;
yo=y0;
 
index=[0];
 
for i=1:Ntotal
 index=[index i];
 H=henon_map(a,b,xo,yo);
 xn=H(1);
 yn=H(2);
 X=[X xn];
 Y=[Y yn];
 xo=xn;
 yo=yn;
end
 
figure(1)
clf
plot(X(1),Y(1),'b.','LineWidth',1,'MarkerSize',22)
hold on
plot(X(Ntransient+1:Ntotal),Y(Ntransient+1:Ntotal),'r.','MarkerSize',5)
plot(X(1:Nshow),Y(1:Nshow),'b--*','LineWidth',1)
%axis([-2.5 2.5 -2.5 2.5]); 
axis tight;
set(gca,'FontSize',[fs]);
xlabel('x_n'); ylabel('y_n');
title(['Henon map: a=',num2str(a),', b=',num2str(b),', (x_0,y_0)=(',num2str(x0),',',num2str(y0),')']);
 
fprintf('Press any key to continue...\n');
pause
 
figure(2)
clf
 
disp=150;
subplot(2,1,1)
plot(index(1:disp),X(1:disp),'r-o')
axis tight
set(gca,'FontSize',[fs]);
xlabel('n');ylabel('x_n');
title(['Henon map: a=',num2str(a),', b=',num2str(b),', (x_0,y_0)=(',num2str(x0),',',num2str(y0),')']);
subplot(2,1,2)
plot(index(Ntotal-disp:Ntotal),X(Ntotal-disp:Ntotal),'r-o')
axis tight
set(gca,'FontSize',[fs]);
xlabel('n');ylabel('x_n');

[lopppez]

salut;
en fin je trouve et dans ce forum la réponse à ma question qui donne l'histogramme d'une suite de dimension 2. et le voila:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
N = 10000 ;% nombre d'itération
% pre-allocation
x = zeros(N,1);
y = zeros(N,1);
% conditions initiales
x(1) = 0.0004;
y(1) = 0.01;
% Récurrence
for i = 1:N-1;
    x(i+1) = 1+y(i)-1.4*x(i)^2 ;
    y(i+1) = 0.3*x(i);
end
y
hist(y)

il donne bien l'histogramme pour des n petit MAIS lorsque j'augmente le nombre de par exemple 10.000.000.000 il donne l'erreur suivante

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
??? Error using ==> zeros
Maximum variable size allowed by the program is
exceeded.
 
Error in ==> updfed3 at 3
x = zeros(N,1);
 
??? Error using ==> zeros
Maximum variable size allowed by the program is
exceeded.

j'ai essaie de cherche dans le help je trouve que c'est un probleme de (int8;int64....double precision..)
mais j'arrive pas a résoudre le problème
merci d'avance

[Jean Dumoncel]

Bonjour,

Pour te rendre compte de la taille de la variable que tu essaie de créer :
Comment déterminer approximativement la quantité de mémoire nécessaire pour le stockage d'une variable ?

Par défaut le format des variables créées par zeros est le format double. Tu peux essayer le format single :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
x = zeros(N,1,'single');
y = zeros(N,1,'single');

tu perdras un peu en précision et ce n'est pas sûr que cela suffise pour des N trop grand.

[lopppez]

bonsoir;
ça marche pas j'ai toujours le message d'erreur suivant:

??? Error using ==> zeros
Out of memory. Type HELP MEMORY for your options.

Error in ==> histogra at 13
x=zeros(n+1,1,'single');

[Gakusei]

Bonjour,

J'ai regardé très vaguement votre poste. Finalement vous avez une matrice avec des valeurs proche de zeros voir nul.
Une idée pour ne pas prendre trop d'espace mémoire:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

doc sparse

J'espère que c'est en rapport.

[lopppez]

salut

ça marche pas et bloque lorsque je prend n plus grand

[Invité]

Est-ce que tu te rends compte de la taille de la variable que tu demandes?
As-tu fait le calcul avec le lien de magelan? Ta mémoire physique n'est pas extensible.
Regarde ce que retourne

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

memory

[Gakusei]

il donne bien l'histogramme pour des n petit MAIS lorsque j'augmente le nombre de par exemple 10.000.000.000 il donne l'erreur suivante

En fait j'ai regardé un peu plus en détail ton problème. Je ne vois pas l'intérêt de créer une matrice de 10.000.00... Car après tu va afficher son histogramme de 10.000.00.. valeurs. Visuellement tu verras pas beaucoup (voir pas de du tout) de différence avec un affichage de 10000 valeurs.

Si en revanche tu veux que tes valeurs de 'y' aillent de 0 à 10000... pour l'affichage avec un pas large, il te suffit d'en prélever dans ta matrice une sur 1000 par exemple.

Une idée d'algorithme:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
N = 10000000000 ;% nombre d'itération
pas=10000;
% pre-allocation pas obligatoire
%x = zeros(N/pas,1);
%y = zeros(N/pas,1);
% des entier ecrasé au fur et à mesure
xi = 0.0004;
yi = 0.01;
cmpt=1;
% Récurrence
for i = 1:N-1;
    xi = 1+yi-1.4*xi^2 ;
    yi = 0.3*xi;
    if mod(i,pas)==0
       x(cmpt)=xi;
       y(cmpt)=yi;
       cmpt=cmpt+1;
    end
end
y
hist(y)

Cependant, si tu veux une valeur à un 'i' donnée de ta fonction, il existe un autre moyen.
On peut remarquer que ta suite est une suite récurrente d'ordre 2 donc tu peux calculer sont polynôm caractéristique et faire la démarche décrite sur ce lien:
http://fr.wikipedia.org/wiki/Suite_r..._lin%C3%A9aire

[Invité]

il te suffit d'en prélever dans ta matrice une sur 1000 par exemple.

Où est l'intérêt? Autant prendre N = numel(1:pas:10000000000); ce qui revient au même.
Sinon une solution serait de construire l'histogramme manuellement dans la boucle.

Remarque:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
xi_tmp = xi;
xi = 1+yi-1.4*xi^2 ;
yi = 0.3*xi_tmp;

Sinon on prend la "nouvelle" valeur

[Gakusei]

Où est l'intérêt? Autant prendre N = numel(1:pas:10000000000); ce qui revient au même.

En fait non car on doit quand même créer la suite et pour la créer il faut les éléments précédents. Quand je disais prélever 1 sur 1000, c'était prélever 1 élément de la suite tous les 1000 itérations. D'ailleurs c'est ce que fait le script que j'ai mis. Avec en effet l'erreur que tu as fais remarqué sur xi_tmp, autant pour moi.

Par contre grossière erreur de ma part pour la définition de la suite, ce n'est pas une suite linéaire d'ordre 2 car x(i)^2 apparait.

[Invité]

Oui j'avais compris cela, mais prendre N = numel(1:pas:10000000000); échantillons à la suite ou espacés de 1000 itérations revient pour moi au même: on aura toujours un tableau de taille dépendant de N.
Bon après, je ne connais pas assez l'attracteur de Hénon pour te dire si au niveau statistique cela changera quelque chose, mais ce que je veux dire c'est que au final, si on augmente trop N, on aura toujours un problème de mémoire.