Logistic 混沌模型的Matlab实现

上一篇 / 下一篇 2007-10-30 10:45:51

!BX$K6Imn0% Author: Thomas Lee振动资讯iF+[r$H&a

1h(i8Q3Nf?"w7l0% E-mail: lixf1979@126.com % Corresponding: School of Mathematics, Physics and Software Engineering, Lanzhou Jiaotong University, Lanzhou 730070, China

%Progran Logistic bifurcation 1.m振动资讯{2{;z*ZV}

x=0.6;

u=2.6:0.001:4;振动资讯G+E1MeulmF3Lb

for j=1:300

x=u.*(x-x.^2);

end振动资讯C'a"ng zP q*r iG

for i=1:200

x=u.*(x-x.^2);

%if i>=1900

%plot(u,x,'k.','markersize',1);振动资讯wWP7jS$f/j#p0O

hold on;

%end%振动资讯N(VV(k-~ap$]

end振动资讯)P^~1k$B X N/N

%Progran Logistic bifurcation 2.m振动资讯(IT.Ru|'Ix F

x=0.1;振动资讯,S.L$b%mLd MF'Vs

u=2.6:0.001:4;振动资讯fqZyK.OBy~

for j=1:300

x=u.*(x-x.^2);

end振动资讯3bFwD'N:I`C

for i=1:2000振动资讯Jr3P @4a$Q

x=u.*(x-x.^2);振动资讯({1h(~{dR

if i>=1900

plot(u,x,'k.','markersize',1);

hold on;

end

end振动资讯7pV#gm9s3e

振动资讯%X1W-i"zQ(D

振动资讯9z K\P qC

%Progran Logistic Phasespace.m振动资讯2m(`pQ]*`

% solve and plot logistic map:振动资讯U%^7Y3S"X*T})v

% x(n+1)=r*x(n)*(1-x(n))

% number of steps:

N=100;

% initial conditions:

x(1)=0.2;

% parameter:振动资讯Gk5_7Fa.{*s-D

% r=3.8; % chaos振动资讯8ycg5G.P a[

% r=2.0; % steady state振动资讯 Mq5qZz%M.~7{V9q,vZ

% r=3.2; % period 2

r=3.5; % period 4?振动资讯(cB7^]:Sg

for n=1:N

x(n+1)=r*x(n)*(1-x(n));振动资讯 J8eNmZ[/I

end

plot(x,'k')振动资讯J'v2F3[:G1So"`

title('logistic map')

xlabel('n')

ylabel('x(n)')振动资讯 v |Pk6bb0t?#{

figure振动资讯3@ ~lTl k1W4w

plot(x(2:N),x(1:N-1),'k.','markersize',5)

title('logistic map')

xlabel('x(n+1)')

ylabel('x(n)')

主要运行结果：振动资讯'cgkl%bmC*R_

%r=3.8:产生的是混沌图相振动资讯.x%lX&|&v

%r=2.0时产生的是稳定状态振动资讯q)YSKwu

%r=3.2时，产生的是二周期的振动资讯:ew8[$Nc

%r=3.5时，产生的是四周期的

振动资讯r7| g8ON5Jn&^?

%Progran Logistic Lyapunov 1.m振动资讯,a(R;y4fny

%computes lyapunov exponent

r = 4.0;振动资讯#OCVV\#R4w

%%%%% first converge to attractor振动资讯\k9HFL.v(D

x = rand; %choose random initial data振动资讯-~G"[WPW6? nl

for i=1:200;振动资讯4AjE+xW

x = r*x*(1-x); %apply logistic map振动资讯 I*]u-A:`1T

end;

%%%%% find (1/n) Sum ln(f'(x_i)) - limit should be lyapunov exponent

sum = 0;

for i=1:500;振动资讯Fsf"d*w)Sd$J

x = r*x*(1-x); %logistic map

sum = sum + log(abs(r - 2*r*x));振动资讯7vlqD%\Y$EO

sum/i; %print out average振动资讯~!w$?9iO3X

average(i) = sum/i; %collect into vector for later plotting

end;

plot(average,'k'); %plots sucessive averages so a trend can be observed

%Progran Logistic Lyapunov 2.m

%computes lyapunov exponent

r =2.0:0.001: 4.0;

%%%%% first converge to attractor振动资讯%c9o*jm{,n

x = 0.2; %choose random initial data

for i=1:200;

x = r.*x.*(1-x); %apply logistic map振动资讯yh:WMZ-n

end;振动资讯2lrLn@]O.ke

%%%%% find (1/n) Sum ln(f'(x_i)) - limit should be lyapunov exponent振动资讯'e"FbL7K2X+h

sum = 0;振动资讯B.??"tP'SH6RG@

for i=1:50000;

x = r.*x.*(1-x); %logistic map振动资讯U N+J[]9E/s#\.eM

sum = sum + log(abs(r - 2*r.*x));

sum/i; %print out average振动资讯 N'e$]JS;M

average= sum/i;

plot(r,average,'k'); %plots sucessive averages so a trend can be observed

%collect into vector for later plotting振动资讯I4M0l+z"j(YgT~]

end;

hold on

line(r,0)

振动资讯3WV6rxz4s2I

%Progran Logistic Cobweb 1.m

function logistic_cobweb(r,N)

% m-file to test iterative maps: implements cobweb construction, including

% transients, for logistic map振动资讯9g7t-{-]4b%~

yy = 0:0.05:1; ff = r.*yy.*(1-yy);振动资讯Z"b*RB:uo f

figure(1) ; clf ; 振动资讯.EW4}3ty8K }2y

plot(yy,ff,'-',yy,yy,'--')

hold on

x = 0.1; % Initial value

Xn = [];振动资讯8dE Fz"w:p'q

%r = 2; % Bifurcation parameter

x1 = r*x*(1-x);振动资讯kZ6q$j7J$E _;t6\

Xn = [Xn x1];振动资讯;X$S)`$\"I/CN-yB#l

plot([x x],[0 x1],'r-',[x x1],[x1 x1],'r-')振动资讯fUK$X;|[!?

x = x1;振动资讯 @ mAz:A,t R[-H4m-`d

for i = 2:N振动资讯F0Aj.j9G"t?2S-s

x1 = r*x*(1-x);

Xn = [Xn x1];

plot([x x],[x x1],'r-',[x x1],[x1 x1],'r-')振动资讯?,IUX*Rm

x = x1;振动资讯PPJ)Hct{

end

hold off振动资讯Cu"_mKy&Yz

〉〉logistic_cobweb(4,200)振动资讯K9g3x{Z

振动资讯u7~u]l/h*V

%Progran Logistic Cobweb 2.m振动资讯I2r6|!S"ty4~d*W_/^j B

%% Renormalisation of logistic map f振动资讯"M9_-YZu&y3` dY-C

r=input('value of r: ');振动资讯 n1`:P)gYglU

clf reset

%% The left-hand plot shows f and f^2

subplot(2,1,1)振动资讯"jYE4WYlsz

x0=0:0.025:1;

%% y=f(x0) and z=f^2(x0)振动资讯q0w#V0A"Z:}:\Vg

y=r.*x0.*(1-x0);

z=r.*y.*(1-y);振动资讯5AI)We^`B)G'z

%figure

plot(x0,y,'k')振动资讯D$M6S xjP

hold on

plot(x0,x0,'k')

plot(x0,z,'k', 'LineWidth', 2)

axis equal

%% Now find where f^2(nxl)=1-(1/r)

%% using the interval halving method振动资讯Y Z/n8i)B

nxl=1-(1/r);

g=inline('(r^2*x*(1-x)*(1-(r*x*(1-x))))-1+(1/r)','x','r');振动资讯0y y.n$Im:r

xl=1-(1/r)+0.00001; xr=1;

d=1;

while d>0.0001振动资讯pTD%ZynO

m=(xl+xr)/2;

if g(xl,r)*g(m,r)<0

xr=m;

else振动资讯7ByI5g(ri/?W/q

xl=m;振动资讯NN;ceY1JOM)z

end

d=(xr-xl);振动资讯5q@x1N QD#X2_r%^

end

nxr=m;

%% Draw the little `invariant box'

plot([nxr nxr],[nxr nxl],'k')振动资讯!D3R/tmq/~2d

plot([nxr nxl],[nxl nxl],'k')

plot([nxl nxr],[nxr nxr],'k')振动资讯Bo)Y6P"q8Pn]W

plot([nxl nxl],[nxr nxl],'k')

%xlim([0 1])振动资讯ob K8A4tS!i

%ylim([0 1])

%% On the right-hand side draw an expanded version of little box振动资讯"\$uj9Hr Zof

subplot(2,1,2)

x2=linspace(nxl,nxr,100);

y2=r.*x2.*(1-x2); y2=r.*y2.*(1-y2);振动资讯BJr&@:G~B-[

plot(x2,y2,'k','LineWidth', 2);

hold on

plot(x2,x2,'k');

axis equal

%xlim([nxl nxr])振动资讯)lza8E!j-b9ej:WL)Z

%ylim([nxl nxr])振动资讯X n8PD8}1db6C

振动资讯dv;n?p3O;z@

全部脚印 不留脚印留下脚印:

导入论坛引用链接收藏分享给好友管理举报

TAG:

引用删除 liliangbiao / 2008-01-15 21:49:27: 原帖由lijh12于2008-01-08 15:05:12发表
不错，很不错

Thanks!

引用删除 lijh12 / 2008-01-08 15:05:12: 不错，很不错

引用删除 liliangbiao / 2007-11-02 16:43:35: You are welcome!
Any constructive comments, I would like!

引用删除 shenyongjun / 2007-11-01 12:29:21: 错了，重来，应该是5分

引用删除 shenyongjun / 2007-11-01 12:28:44: 评 1 分

引用删除 zhlong / 2007-10-30 14:14:38: 评 5 分

查看全部评论

关于作者

liliangbiao

Logistic 混沌模型的Matlab实现

全部脚印 不留脚印 留下脚印:

关于作者

全部脚印不留脚印留下脚印: