求Feigenbaum图源程序

sihaiwuya

那位高手帮一下忙，这是混沌动力学的的一个分岔图，求matlab的源程序。不胜感激，刚学matlab，还望高人指点。

[ 本帖最后由 lxq 于 2006-12-15 14:20 编辑 ]

MVH

是不是logistic模型分叉图

方法一：

1. clear;
   
2. clf;
   
3. u=2.6:0.001:4.0;
   
4. x=0.1;
   
5. for i=1:300
   
6. x=u.*(x-x.^2);
   
7. end
   
8. for j=1:80
   
9. x=u.*(x-x.^2);
   
10. plot(u,x,'k.','markersize',2)
    
11. hold on;
    
12. end
    
13. grid on;

复制代码

方法二：

1. clear;
   
2. niter=2500;
   
3. hold on;
   
4. for alpha=1:0.005:4
   
5. y(1)=rand(1);
   
6. x(1)=alpha;
   
7. for j=1:niter
   
8. y(j+1)=alpha*y(j)*(1-y(j));
   
9. x(j+1)=alpha;
   
10. end
    
11. h=plot(x(200:niter+1),y(200:niter+1),'k.');
    
12. set(h,{'MarkerSize'},{0.5});
    
13. end
    
14. xlabel('--alpha-->');
    
15. title('Bifurcation Diagram of Logistic Map');

复制代码

方法三：

1. function rep = itera(x)
   
2. % calcul de la fonction logistique
   
3. global lamb % lamb est en commun avec le premier prog.
   
4. rep = lamb.*x.*(1-x);
   
5. % Dessine l'arbre de Feigenbaum,
   
6. %c.a.d. l'arbre de bifurcation ddes iterations
   
7. % de l'application logistiquefigure(1); % choisit une fenetre
   
8. clf;hold on; % nettoie l'ecran, fixe les le cadre
   
9. % entre deux affichages pour une surimpression
   
10. a=1;b=4; % le ";" evite l'affichage du resultat
    
11. axis([a,b,0,1]); % [xmin,xmax,ymin,ymax] est un vecteur qui
    
12. % defini le cadre
    
13. global lamb; % met une variable en commun
    
14. N=500;
    
15. for i = 1:N, % noter la virgule
    
16. x = 0.71;
    
17. lamb=a+i/N*(b-a);
    
18. Y = [ ]; % initialise un vecteur (vide)
    
19. for j = 1:200,
    
20. x = itera(x); % itera est une fonction utilisateur
    
21. Y = [Y,x]; % ajoute une composante x au vecteur Y
    
22. end; % fin du for j
    
23. Y=Y(100:200); % extrait les composantes de 100 a 200
    
24. plot(ones(size(Y))*lamb,Y,'r.');
    
25. if rem(i,10)==0, drawnow; end;
    
26. % affiche le dessin toutes les 10 it\'erations
    
27. % (rem = remainder = reste de la division
    
28. % de i par 10)
    
29. end;

复制代码

方法四:

1. r=linspace(2,4,500);
   
2. initial=linspace(0.3,0.7,50);
   
3. trans=200;
   
4. k=1;
   
5. while k<50
   
6. 
   
7. i=1;
   
8. while i<=500
   
9. x(1)=initial(k);
   
10. j=1;
    
11. while j<=(trans-1)
    
12. x(j+1)=r(i)*x(j)*(1-x(j));
    
13. j=j+1;
    
14. end
    
15. a(i)=x(trans);
    
16. b(i)=x(trans-1);
    
17. c(i)=x(trans-2);
    
18. d(i)=x(trans-3);
    
19. i=i+1;
    
20. end
    
21. plot(r,a,'.','Color','red');hold on;
    
22. %set(findobj(gca,'Type','line'),'Color','red');
    
23. plot(r,b,'.','Color','blue');hold on;
    
24. %set(findobj(gca,'Type','line'),'Color','blue');
    
25. plot(r,c,'.','Color','yellow');hold on;
    
26. plot(r,d,'.','Color','green');
    
27. k=k+1;
    
28. end

复制代码

sihaiwuya

谢谢了！

xsj3917

非常感激

zhlong

非常感谢

zhangyandlut

真是感激呀！！！