Logistic 混沌模型的Matlab实现

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

@'S_0E.z9p8{(w0% Author: Thomas Lee振动资讯r+N9W"?7@ ba0u0Lps;]

% E-mail: lixf1979@126.com % Corresponding: School of Mathematics, Physics and Software Engineering, Lanzhou Jiaotong University, Lanzhou 730070, China 振动资讯Vb0Ut_nA6Dq

%Progran Logistic bifurcation 1.m

x=0.6;

u=2.6:0.001:4;振动资讯3kl*V;US

for j=1:300振动资讯@+pK){y3w

x=u.*(x-x.^2);振动资讯+}^aMF%~2R~

end振动资讯#mKj8^G.i;hB

for i=1:200

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

%if i>=1900 振动资讯Zx#V7hwD*|8N

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

hold on;

%end%振动资讯;y/f;K/is

end振动资讯8hB.aqU,c

振动资讯n.dB7J4G/x s)[$a S

%Progran Logistic bifurcation 2.m振动资讯F/[qv3^"h?6PugA%cp

x=0.1;

u=2.6:0.001:4;振动资讯'uM,\?Ei~#D*Jv2X,w

for j=1:300

x=u.*(x-x.^2);振动资讯QSvY;m0v

end振动资讯&DSfPs|

for i=1:2000振动资讯&VWftA!Qe*b&}K.Q3^-|

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

if i>=1900

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

hold on;

end

%Progran Logistic Phasespace.m

% solve and plot logistic map:振动资讯YAFh1{-]+\E

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

% number of steps:

N=100;振动资讯bFV"r8Fb@

% initial conditions:

x(1)=0.2;

% parameter:

% r=3.8; % chaos

% r=2.0; % steady state

% r=3.2; % period 2振动资讯C/j.|7j_

r=3.5; % period 4?振动资讯 kt$G B;\M7Cz0TB\|"]

for n=1:N

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

end振动资讯N3DP*gA ~9e"Y%XeCm

plot(x,'k')振动资讯+oYE4n(w`ken

title('logistic map')振动资讯#@&y0O*ox$@

xlabel('n')

ylabel('x(n)')

figure

plot(x(2:N),x(1:N-1),'k.','markersize',5)振动资讯cP8yx/HB'N_)i

title('logistic map')振动资讯FZ&ptG4j;?

xlabel('x(n+1)')振动资讯-N9lZD^nb?#i`2^

ylabel('x(n)')

主要运行结果：振动资讯x k~,\+H;lu `

%r=3.8:产生的是混沌图相

%r=2.0时产生的是稳定状态

%r=3.2时，产生的是二周期的

%r=3.5时，产生的是四周期的振动资讯Q4^CH;e{L/] L{

振动资讯8{"z!w!RV2F V

%Progran Logistic Lyapunov 1.m

%computes lyapunov exponent 振动资讯3U-WZ2n7dBQh~c

r = 4.0;

%%%%% first converge to attractor

x = rand; %choose random initial data振动资讯^!]"m4k(R+KX&Do:P

for i=1:200;振动资讯y1Q~c g)[!B

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

end;振动资讯a(jg){}3G

%%%%% find (1/n) Sum ln(f'(x_i)) - limit should be lyapunov exponent振动资讯aW&T6M4G%tq!a

sum = 0;振动资讯{wYJ,OFC+o*h7z

for i=1:500;

x = r*x*(1-x); %logistic map振动资讯0\J eQ8P V8EK v#G

sum = sum + log(abs(r - 2*r*x));振动资讯 a0va3`cKH

sum/i; %print out average

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

end;

plot(average,'k'); %plots sucessive averages so a trend can be observed振动资讯S(C2wj+uDt5Y}H

振动资讯2R r[$c9S

%Progran Logistic Lyapunov 2.m振动资讯*wl&K:xx-B M

%computes lyapunov exponent

r =2.0:0.001: 4.0;振动资讯$AF6lt \)x7wvr

%%%%% first converge to attractor

x = 0.2; %choose random initial data

for i=1:200;

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

end;

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

sum = 0;

for i=1:50000;振动资讯Z)~*oqq)q^P

x = r.*x.*(1-x); %logistic map振动资讯8MI$X7]Z$Qbmz

sum = sum + log(abs(r - 2*r.*x));振动资讯\D5t@(aZ

sum/i; %print out average振动资讯f&Dysc${l

average= sum/i;

plot(r,average,'k'); %plots sucessive averages so a trend can be observed振动资讯-\Zhy iDi

%collect into vector for later plotting

end;振动资讯M|q ^expX

hold on振动资讯$_M;t?"\Y |m5LsK Q

line(r,0)振动资讯*^ u5c:\X |8Q

振动资讯k@ wKu

%Progran Logistic Cobweb 1.m

function logistic_cobweb(r,N)振动资讯*bJ ]dGz

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

% transients, for logistic map振动资讯`4S[:\:|V.F9^/Iy

yy = 0:0.05:1; ff = r.*yy.*(1-yy);

figure(1) ; clf ;

plot(yy,ff,'-',yy,yy,'--')振动资讯3f!u8^[,{ E

hold on振动资讯:o2Iv`q!Rp

x = 0.1; % Initial value振动资讯K/EcLpL8`5e#k

Xn = [];振动资讯E.f)N"T3st

%r = 2; % Bifurcation parameter振动资讯!U!M+@&m#v4d(B&I2Z

x1 = r*x*(1-x);振动资讯FLf1?m4A.Qs

Xn = [Xn x1];振动资讯B(WVp?)a f"Z

plot([x x],[0 x1],'r-',[x x1],[x1 x1],'r-')

x = x1;振动资讯rrc)QXW#cLRA

for i = 2:N

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

Xn = [Xn x1];

plot([x x],[x x1],'r-',[x x1],[x1 x1],'r-')

x = x1;振动资讯n;i(B{.S!_%[YT

end振动资讯xCTO"H-i'~ c

hold off

〉〉logistic_cobweb(4,200)

振动资讯Dk!~mP$q.a{$n

%Progran Logistic Cobweb 2.m

%% Renormalisation of logistic map f

r=input('value of r: ');

clf reset

%% The left-hand plot shows f and f^2振动资讯u$Q$O#x _-b#A

subplot(2,1,1)

x0=0:0.025:1;振动资讯2UYn\kn

%% y=f(x0) and z=f^2(x0)振动资讯_{YF{!ux;l&@+JH

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

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

%figure

plot(x0,y,'k')振动资讯~9af"g7b#_,O

hold on振动资讯/ugY(cvby

plot(x0,x0,'k')

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

axis equal

%% Now find where f^2(nxl)=1-(1/r)振动资讯M NA+E5M'G:Cv%e0RUQ }O

%% using the interval halving method振动资讯0Y/e8MXw2X8}5d

nxl=1-(1/r);振动资讯T}_-K)s@

g=inline('(r^2*x*(1-x)*(1-(r*x*(1-x))))-1+(1/r)','x','r');

xl=1-(1/r)+0.00001; xr=1;振动资讯 GNE$BjW%t

d=1;振动资讯R9D cCagy

while d>0.0001

m=(xl+xr)/2;振动资讯`,?(ddKaP] ?d

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

xr=m;

else

xl=m;

end振动资讯0cr,fH"eq6rQ7T

d=(xr-xl);振动资讯K^,J y6S1o!nD9VZ

end

nxr=m;振动资讯(]/nyzE*uqZO9uG

%% Draw the little `invariant box'

plot([nxr nxr],[nxr nxl],'k')振动资讯n,Z;_r(QOu

plot([nxr nxl],[nxl nxl],'k')振动资讯$V7G%eH/](Qf

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

plot([nxl nxl],[nxr nxl],'k')振动资讯g6Ia/H*z

%xlim([0 1])振动资讯(YA:AX_/z b

%ylim([0 1])

%% On the right-hand side draw an expanded version of little box

subplot(2,1,2)

x2=linspace(nxl,nxr,100);振动资讯4Muk9{%|?#D(m

y2=r.*x2.*(1-x2); y2=r.*y2.*(1-y2);振动资讯P3ixn@\k-F(t fd

plot(x2,y2,'k','LineWidth', 2);振动资讯S/K"A s{ Q

hold on

plot(x2,x2,'k');振动资讯z I.k#l6v]

axis equal振动资讯_/q8^m+{Uh

%xlim([nxl nxr])

%ylim([nxl nxr])振动资讯z*XK&?2D;i^

振动资讯&MP\L6|3W;q

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

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

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实现

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

关于作者

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