Lozi混沌映射系统的Matlab实现

% Author: Thomas Lee振动资讯Gz+K9BpE
% E-mail: lixf1979@126.com振动资讯$vy9MbdlZ1\)sXL2r
% Corresponding: School of Mathematics, Physics and Software Engineering, Lanzhou Jiaotong University, Lanzhou 730070, China

PDF\8r'f3p0% if you want to get more information, please refer to one of the published articles of the author :

6d0[%}fM0`.C`
|/p M+K0%李险峰等.Lozi混沌映射的线性反馈控制.河北师范大学学报. 2007,31(4):479-483.振动资讯1b@|}_m

1. Phase Portrait:振动资讯pT??r*g

function newx=Lozi(x);

% Lozi方程[差分方程]%newx=Lozi([x;y;p;q]);振动资讯Tc/|)y-X+Vf J

%方程如下:%x(k+1)=-p*x(k)+y(k)+1;% y(k+1)=q*x(k);

%定义列向量newx及函数Lozi(x)！

newx(1,1)=-x(3)*abs(x(1))+x(2)+1;

newx(2,1)=x(4)*x(1);振动资讯y M*@ z-_6z;r

newx(3,1)=x(3);

newx(4,1)=x(4);

X=[1;0.1;1.7;.5];振动资讯 x{/t7\&@5t

%前面两个赋初值,后面两个给出了系统中的

%常数项的值！与newx=Lozi(x)中给出的列向量相对应！

 Y=[]; %给出一个空数组，以便存储得出的迭代解！

 for i=1:10000 %循环即迭代次数

    X=feval(@Lozi,X);%调用主函数

    Y(i,:)=X(1:2,1);

    %将每一次迭代后的数组X的上面两个变量存储！振动资讯z4hK.t~.L

    %因为下面的两个常数项进行迭代时保持原来的数值不便！

 end振动资讯C-`C7[W]Y

    plot(Y(:,1),Y(:,2),'.','markersize',1);振动资讯~[B4W6O;K!bx:ga

    %得到Y=[Y(:,1),Y(:,2)],分别指代每一次迭代后的振动资讯BU n
G(H5T VPq

    %数组X的上面两个变量,然后点画图.振动资讯!g3K
nsHUD9PA5`3a

    title('Lozi映射图')%加上图像标题振动资讯bj:P#@r
ai`

    xlabel('x'),ylabel('y')

    %分别给对应的向量随处的坐标命名

%另外还可以使用其它的语句对图像进行进一步的调整振动资讯8vh3N&S`,ilL

2.  Bifurcation diagram 1振动资讯-[aQ,ju%Z

function newx=Lozi(x);

% Lozi方程[差分方程]%newx=Lozi([x;y;p;q]);

%方程如下:%x(k+1)=-p*x(k)+y(k)+1;% y(k+1)=q*x(k);振动资讯 PhIGoq0v4s

%定义列向量newx及函数Lozi(x)！振动资讯.k.ge&d}

newx(1,1)=-x(3)*abs(x(1))+x(2)+1;振动资讯6aF+[ BH4J2z8GK3d

newx(2,1)=x(4)*x(1);

newx(3,1)=x(3);振动资讯\-jr lXE5cE

newx(4,1)=x(4);

 振动资讯1A*s3Chpb)z*`#F

  Z=[];%给出一个空数组，以便存储得出的迭代解！振动资讯.vb t;Fj+YOU-U;Vx/E

  for p=linspace(0,1.7,1700);振动资讯wj&I?Jl

      %将分岔参数用linspace在区间[0,1.7]

      %等分成1700份

      x=[1;0;p;.5];

      %前面两个赋初值,后面两个给出了系统中的振动资讯$?'h9N}$}3r3U^B

      %常数项的值！与newx=Lozi(x)中给出的列

      %向量相对应,其中值得注意的是p是连续变化振动资讯|l%[u1cX8t:l

      %的向量.

      for k=1:500 %循环即迭代次数

          x=Lozi(x); %调用主函数

          if k>400

              %取迭代后的最后一百个点,这里认为振动资讯/]R,o!g)?m ~

              %前400个点为迭代过程中的瞬态响应！

              Z=[Z,p+x(1)*i];振动资讯ZQ*T9C5{'o
D0X|.U2K

               %将P的值以及在这个值所对应最后400次迭代振动资讯8H#|3H]!X+S0m2`

               %后的数组X的第一个变量x存储为一个二维列向量组！

          end振动资讯5u+VJb1^1C"^^Z#c

      end振动资讯a hum2YG!ed

  end

  plot(Z,'.','markersize',1)

   %得到Z=[Z(:,1),Z(:,2)],分别指代P的值以及在这个振动资讯!Z%DE)g;k M p

   %值所对应最后400次迭代后的数组X的第一个变量x,然后点画图.振动资讯3Xf~;P*zq.GG ]$s

  title('Lozi映射分岔图')%加上图像标题振动资讯1x3U-N:X!q6jn

  xlabel('p'),ylabel('x')振动资讯'psT%}*kt)i;b

  %分别给对应的向量随处的坐标命名振动资讯(x0v7O [p-VW+Vd

  %另外还可以使用其它的语句对图像进行进一步的调整

&O.l'vS*r F3s$s0 

3. Bifurcation diagram 2振动资讯AiW$T+^j"mh

function newx=Lozi(x);振动资讯K,a)z(\5P-z+{1H

% Lozi方程[差分方程]%newx=Lozi([x;y;p;q]);振动资讯.o j~T!N6T

%方程如下:%x(k+1)=-p*x(k)+y(k)+1;% y(k+1)=q*x(k);

%定义列向量newx及函数Lozi(x)！

newx(1,1)=-x(3)*abs(x(1))+x(2)+1;振动资讯UhxZ/n,Ze v
yC

newx(2,1)=x(4)*x(1);振动资讯Lv.R|nO)~S.g8^

newx(3,1)=x(3);

newx(4,1)=x(4);

 P=[];%给出一个空数组，以便存储P的值！振动资讯(K{8u)[H7a

  Z=[];%给出一个空数组，以便存储得出的迭代解！振动资讯'M7Fk rJr
iV

  for p=linspace(0,1.7,1700);振动资讯 LV9C%V"q*w!Zrc

      %将分岔参数用linspace在区间[0,1.7]

      %等分成1700份

      x=[1;0;p;.5];振动资讯Z}8Rmz'W

      %前面两个赋初值,后面两个给出了系统中的

      %常数项的值！与newx=Lozi(x)中给出的列振动资讯E5cX2|wbZ;W-P:P

      %向量相对应,其中值得注意的是p是连续变化

      %的向量.振动资讯1`P
h xcUlYH

      for k=1:500

%循环即迭代次数,但是遗憾的是这里的迭代次数太少,

%因为如果太多的话,由于上面的p点太多,运行的时间

%就越长,但是为了得到较为精确的迭代解,建议适当的增加迭代次数振动资讯8ZZV+QDR8J

%估计在10000左右,而这时候取最后的1000个点即可！振动资讯(a2lK'gL*[Me

          x=Lozi(x); %调用主函数

          if k>400振动资讯k^(LCN3mm

              %取迭代后的最后一百个点,这里认为振动资讯0z&m1A PH,K W g4MM

              %前400个点为迭代过程中的瞬态响应！振动资讯$\dLx4w2eDZ

              Z=[Z;x(1)];

              %这个值所对应最后400次迭代振动资讯:zWK k Rj})N?+Sf

               %后的数组X的第一个变量x存储为一个二维列向量组！振动资讯T8Y8Iq)G-ks F)b.m

              P=[P;p]; %将P的值存储到空变量！

          end

      end振动资讯p!k }c8QZ'jTm
x&X.}

  end振动资讯8M7ou:w
j
S,e

  plot(P,Z,'.','markersize',1)%点画图振动资讯j7Kt0Zp
q%|3L6O:e

  title('Lozi映射分岔图')%加上图像标题振动资讯"TG `DA^

  xlabel('p'),ylabel('x')振动资讯#cQ K:Zec-D'cw_)^

  %分别给对应的向量随处的坐标命名振动资讯,I6| e2d&?(QAp7M4ms+M

  %另外还可以使用其它的语句对图像进行进一步的调整

4 Lyapunov-exponent spectrum振动资讯V6Itk!Ba)f
i V

function newx=Lozi(x);振动资讯)Nl K
R5a[.`Z!ib

% Lozi方程[差分方程]%newx=Lozi([x;y;p;q]);

%方程如下:%x(k+1)=-p*x(k)+y(k)+1;% y(k+1)=q*x(k);

%定义列向量newx及函数Lozi(x)！

newx(1,1)=-x(3)*abs(x(1))+x(2)+1;

newx(2,1)=x(4)*x(1);

newx(3,1)=x(3);振动资讯n/[,Qk&K+N5Ed

newx(4,1)=x(4);

  d0=1e-8;%定义极小的扰动,还可以定义更小振动资讯%P"G!G+}"h@3qm:GD

  Z=[];%给出一个空数组，以便存储得出的迭代解！

  P=[];%给出一个空数组，以便存储P的值！振动资讯3G0?[[S6f&E W

  for p=linspace(0,1.7,1700)振动资讯5F-k)Gd"A8{-DU9F

      %将分岔参数用linspace在区间[0,1.7]

      %等分成1700份

      le=0;%定义指数的初始值为0.

      lsum=0;%定义和的初始值为0.振动资讯?2Wu-u7i:hkXE

      x=[1;0;p;.5];%初始值

      x1=[1;d0;p;.5];%扰动初始值振动资讯/i0ChPUCH1]

      for k=1:1000

          %循环即迭代次数,但是遗憾的是这里的迭代次数太少,振动资讯VTu7i$P.N[R

%因为如果太多的话,由于上面的p点太多,运行的时间振动资讯{j ]k)fBI

%就越长,但是为了得到较为精确的迭代解,建议适当的增加迭代次数振动资讯Xw z7~.P3j r)ok

%估计在10000左右,而这时候取最后的1000个点即可！

          x=Lozi(x);%调用初始值下的主函数

          x1=Lozi(x1);%调用扰动初始值下的主函数

          d1=sqrt((x(1)-x1(1))^2+(x(2)-x1(2))^2);振动资讯-se#s5w$_C

          %定义两条规线在欧氏空间中的距离

          x1=x+(d0/d1)*(x1-x);振动资讯9iH;ABmpJJ2Gg

          %选择基准轨线振动资讯b9n.ei"h1j1u
B[

          if k>500振动资讯dTJ|R*^6~~

              %取迭代后的最后五百个点,这里认为

              %前500个点为迭代过程中的瞬态响应！振动资讯z4WGd3Zh

              lsum=lsum+log(d1/d0);振动资讯'VnGM6vj$X$Cs

              %求后500积分后的特征值的和振动资讯5WS'@(c7Z foH

          end振动资讯-J!or&Z"~4M5h

      end振动资讯6a&V e QJTf a2m

      le=lsum/(k-500);%求平均值得到指数值振动资讯'zuO*^1Y|L

      if k==1000 %取最后一个点

      Z=[Z,le];%存储最后一个指数值振动资讯o~D)iP3?%p

      P=[P,p];%存储p值(以对应指数)

      end振动资讯k,q$j%qfP!KLOQJ2X4W

  end振动资讯3h;I|od{.R-Afl

  plot(P,Z,'-')%线画图

  title('Lozi最大Lyapunov指数图')

  xlabel('p'),ylabel('lyapunov')

:]Ox6i:A2CKB05 Chaos control to one of the fixed points

)V-X['n-^TTD05.1 Fixed points

q=0.5;

for p=0:0.001:1.7

x1=1/(p-q+1);振动资讯+g2y Ek*|B'c9v

x2=sign(x1);

K=max(q-1-p*x2,q-1+p*x2);

plot(p,K,'r')

hold on

K1=0.25*(p^2*(x2)^2+4*q);

plot(p,K1,'b')

end振动资讯B7W6axw*L

 振动资讯c"S3t-U'~J'c

5.2 Fixed points controller 振动资讯/z'F V7k1x&D)~w H

function newx=Lozi(x);振动资讯vS
ch_)Q7HB_

% Lozi方程[差分方程]%newx=Lozi([x;y;p;q]);振动资讯(ChhT;`1iO

%方程如下:%x(k+1)=-p*x(k)+y(k)+1;% y(k+1)=q*x(k);

%定义列向量newx及函数Lozi(x)！振动资讯;Fn b(aP.VTI

newx(1,1)=-x(3)*abs(x(1))+x(2)+1;

newx(2,1)=x(4)*x(1)-1.222*(x(1)-1/2.2);振动资讯/F)z _;l(g

newx(3,1)=x(3);振动资讯"h"U L7hm

newx(4,1)=x(4);振动资讯?tl$Jdu'f k
U9f0h

 振动资讯2s`Hb6T:r.m.r~E Dbk

 X=[1;0.1;1.7;.5];振动资讯a%{,w;Su p

%前面两个赋初值,后面两个给出了系统中的振动资讯0V.F.R:}]&F~

%常数项的值！与newx=Lozi(x)中给出的列向量相对应！振动资讯"ef.E_5Nu C_}

 Y=[]; %给出一个空数组，以便存储得出的迭代解！

 for i=1:10000 %循环即迭代次数

    X=feval(@Lozi,X);%

Lozi Map

bifurcation diagram 1

bifurcation diagram 2

Lyapunov-exponent spectrum

control chaos to one of the fixed points