Lozi混沌映射系统的Matlab实现

% Author: Thomas Lee
k6a8O rr}4`7K Z&h0% E-mail: lixf1979@126.com振动资讯DN,D1p3B/Y*W9?
% Corresponding: School of Mathematics, Physics and Software Engineering, Lanzhou Jiaotong University, Lanzhou 730070, China

*n0[ pz@F0% if you want to get more information, please refer to one of the published articles of the author :

.Fg#op1h:{)^7b0%李险峰等.Lozi混沌映射的线性反馈控制.河北师范大学学报. 2007,31(4):479-483.

p_R].m t^#d&Q01. Phase Portrait:振动资讯2N v~9a(Y)fd W%L

function newx=Lozi(x);振动资讯 sh,KTF_U

% Lozi方程[差分方程]%newx=Lozi([x;y;p;q]);振动资讯W)zI;i0S
S%c

%方程如下:%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);振动资讯a,a^(oYso

newx(3,1)=x(3);振动资讯
J/G1aG)E`

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

X=[1;0.1;1.7;.5];振动资讯 P
U;go`)]bu
x

%前面两个赋初值,后面两个给出了系统中的振动资讯:J!E:Ty2V }oWV

%常数项的值！与newx=Lozi(x)中给出的列向量相对应！振动资讯}qe{`g

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

 for i=1:10000 %循环即迭代次数振动资讯H H2c:UKy

    X=feval(@Lozi,X);%调用主函数振动资讯2?PK
P(m"BNl]

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

    %将每一次迭代后的数组X的上面两个变量存储！

    %因为下面的两个常数项进行迭代时保持原来的数值不便！振动资讯%u@L S.vH;n\_k

 end

    plot(Y(:,1),Y(:,2),'.','markersize',1);振动资讯!c`8RJ{3I

    %得到Y=[Y(:,1),Y(:,2)],分别指代每一次迭代后的振动资讯u#}s0EM

    %数组X的上面两个变量,然后点画图.振动资讯+{GeHI1{3n

    title('Lozi映射图')%加上图像标题

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

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

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

zyXga2k/@ Wd!e"O02.  Bifurcation diagram 1

function newx=Lozi(x);振动资讯x(|Oe,O0IndL

% 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;振动资讯U2~(B'g(}*X(b;xp

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

newx(3,1)=x(3);振动资讯oCw.x,H(\}

newx(4,1)=x(4);振动资讯}%@5F$CsX-kT

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

  for p=linspace(0,1.7,1700);振动资讯)zu k#Pr+XB

      %将分岔参数用linspace在区间[0,1.7]振动资讯qd\Tl?

      %等分成1700份

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

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

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

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

      %的向量.振动资讯&dI r+v$]Q}s

      for k=1:500 %循环即迭代次数振动资讯2u3H L;S/`

          x=Lozi(x); %调用主函数振动资讯@\/cWL'?^n:C

          if k>400振动资讯'S*]{tQL

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

              %前400个点为迭代过程中的瞬态响应！振动资讯b2de^/?#z?

              Z=[Z,p+x(1)*i];振动资讯7n:o3Vd9pKT

               %将P的值以及在这个值所对应最后400次迭代振动资讯|!f2j \U
v

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

          end

      end

  end振动资讯 E%kT$G v!x,q7x K

  plot(Z,'.','markersize',1)振动资讯
bcK|q0}

   %得到Z=[Z(:,1),Z(:,2)],分别指代P的值以及在这个

   %值所对应最后400次迭代后的数组X的第一个变量x,然后点画图.

  title('Lozi映射分岔图')%加上图像标题

  xlabel('p'),ylabel('x')

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

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

$[},I*SqmH0 

Wc*Jq#~ ?0]*yV03. Bifurcation diagram 2振动资讯r \.{0W#xxy

function newx=Lozi(x);

% Lozi方程[差分方程]%newx=Lozi([x;y;p;q]);振动资讯s[}{Q;q%['b*}

%方程如下:%x(k+1)=-p*x(k)+y(k)+1;% y(k+1)=q*x(k);振动资讯'_^*e#c0bl0c G

%定义列向量newx及函数Lozi(x)！振动资讯Xtns mzd p

newx(1,1)=-x(3)*abs(x(1))+x(2)+1;振动资讯B;skF*G1R;TXO

newx(2,1)=x(4)*x(1);振动资讯F7v |.a!_1KX

newx(3,1)=x(3);振动资讯;]|/e/U)[

newx(4,1)=x(4);振动资讯0U-p P+_%V0`MJ-f

 P=[];%给出一个空数组，以便存储P的值！振动资讯pu0l*t;W]

  Z=[];%给出一个空数组，以便存储得出的迭代解！振动资讯/F'W\u6w'@bL

  for p=linspace(0,1.7,1700);振动资讯j"QW)c
b:Q*~ I

      %将分岔参数用linspace在区间[0,1.7]振动资讯\ gGiT3d c

      %等分成1700份振动资讯;@,}:L5D5G%VNV

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

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

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

      %向量相对应,其中值得注意的是p是连续变化振动资讯k}bf6_9i$zkWN

      %的向量.振动资讯_Ir H'Y6b%A

      for k=1:500

%循环即迭代次数,但是遗憾的是这里的迭代次数太少,振动资讯u3G5B OI%{$F7@y

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

%就越长,但是为了得到较为精确的迭代解,建议适当的增加迭代次数

%估计在10000左右,而这时候取最后的1000个点即可！振动资讯.w:`2RO[*ouB\

          x=Lozi(x); %调用主函数振动资讯$v)Yz}R3O4@

          if k>400振动资讯,x&NoU"I

              %取迭代后的最后一百个点,这里认为振动资讯1jO&UEuY"A

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

              Z=[Z;x(1)];

              %这个值所对应最后400次迭代振动资讯
Yo&X$cW#k

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

              P=[P;p]; %将P的值存储到空变量！振动资讯aX
G,\2u
Je

          end

      end振动资讯.t/P/o(DwW^|/sa5F

  end

  plot(P,Z,'.','markersize',1)%点画图

  title('Lozi映射分岔图')%加上图像标题振动资讯;T.|y!WM9G4N

  xlabel('p'),ylabel('x')振动资讯qjqN t

  %分别给对应的向量随处的坐标命名振动资讯d&AB r"O:VT

  %另外还可以使用其它的语句对图像进行进一步的调整振动资讯m)R Zm$U:t

Ii;|j-O!~2{%T)Dg04 Lyapunov-exponent spectrum

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);

%定义列向量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);振动资讯'|6bv*|eOd

newx(4,1)=x(4);振动资讯"D.p2b1A)kr,o#I

 振动资讯1K3o1e%SA3L*K&e

  d0=1e-8;%定义极小的扰动,还可以定义更小振动资讯K1i;W"U? {`"bH(S

  Z=[];%给出一个空数组，以便存储得出的迭代解！振动资讯:J$B3`wa{ E\cd

  P=[];%给出一个空数组，以便存储P的值！

  for p=linspace(0,1.7,1700)振动资讯,s\9e}$|%gL$V

      %将分岔参数用linspace在区间[0,1.7]振动资讯"O'mN`(Ws7Z{1x:k

      %等分成1700份

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

      lsum=0;%定义和的初始值为0.

      x=[1;0;p;.5];%初始值振动资讯w(YY[-VIf

      x1=[1;d0;p;.5];%扰动初始值振动资讯2GN,m{Kv UP

      for k=1:1000振动资讯O)K:|)]Y
}!V

          %循环即迭代次数,但是遗憾的是这里的迭代次数太少,振动资讯2_5`l;n6]d&RF

%因为如果太多的话,由于上面的p点太多,运行的时间振动资讯-cM'R,a$Z3h

%就越长,但是为了得到较为精确的迭代解,建议适当的增加迭代次数

%估计在10000左右,而这时候取最后的1000个点即可！振动资讯km`6A-z^5Q
a

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

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

          d1=sqrt((x(1)-x1(1))^2+(x(2)-x1(2))^2);振动资讯1ehnVhm7Z

          %定义两条规线在欧氏空间中的距离振动资讯4b3N-rf-m%Ki

          x1=x+(d0/d1)*(x1-x);

          %选择基准轨线振动资讯,u/VE;_u/un.U

          if k>500振动资讯 I`;^?L

              %取迭代后的最后五百个点,这里认为振动资讯k!@5X:o!oX$gc

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

              lsum=lsum+log(d1/d0);

              %求后500积分后的特征值的和

          end振动资讯.p!d1G2qB0XA]T

      end振动资讯5R2W&nw ZC5FJU

      le=lsum/(k-500);%求平均值得到指数值

      if k==1000 %取最后一个点振动资讯|3A^L2W7O,E|6_

      Z=[Z,le];%存储最后一个指数值振动资讯"a@B~5w[

      P=[P,p];%存储p值(以对应指数)振动资讯,t%B#S*xR}3BzQ

      end振动资讯#xh)vZ2i,vy

  end振动资讯H"k};dm+A.r"p:`

  plot(P,Z,'-')%线画图振动资讯{
?.MG@9Dj

  title('Lozi最大Lyapunov指数图')振动资讯&C["{DH3zf

  xlabel('p'),ylabel('lyapunov')振动资讯 j:DeB$vuh*e

5 Chaos control to one of the fixed points振动资讯.mM7^|-T(L8W

5.1 Fixed points

q=0.5;振动资讯&[S4~q~1`

for p=0:0.001:1.7振动资讯!I[9d+H!g

x1=1/(p-q+1);振动资讯9@3n*A.N&UJ#cVHfv

x2=sign(x1);振动资讯b9jo3e W,bw

K=max(q-1-p*x2,q-1+p*x2);振动资讯;ZgcaL$G:O

plot(p,K,'r')

hold on振动资讯cey cC-`H

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

plot(p,K1,'b')

end

 振动资讯*GuuI^9Pm~

5.2 Fixed points controller 振动资讯M0hys*tv I

 振动资讯'h+f,hFPz

function newx=Lozi(x);

% Lozi方程[差分方程]%newx=Lozi([x;y;p;q]);振动资讯)YuA"s @

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

%定义列向量newx及函数Lozi(x)！振动资讯$uP jd/MD}/[Prg

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);振动资讯 S7q8X4m7g+O

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

newx(4,1)=x(4);振动资讯Nf5t+H G*U8U

 X=[1;0.1;1.7;.5];

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

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

 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