% Author: Thomas Lee
% E-mail: lixf1979@126.com
% Corresponding: School of Mathematics, Physics and Software Engineering, Lanzhou Jiaotong University, Lanzhou 730070, China

% if you want to get more information, please refer to one of the published articles of the author :

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

1. Phase Portrait:

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

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

 

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

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

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

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

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

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

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

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

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

 end

     plot(Y(:,1),Y(:,2),'.','markersize',1);

     %得到Y=[Y(:,1),Y(:,2)],分别指代每一次迭代后的

     %数组X的上面两个变量,然后点画图.

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

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

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

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

2.  Bifurcation diagram 1

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

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

 

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

   for p=linspace(0,1.7,1700);

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

       %等分成1700份

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

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

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

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

       %的向量.

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

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

           if k>400

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

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

               Z=[Z,p+x(1)*i];

                %将P的值以及在这个值所对应最后400次迭代

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

           end

       end

   end

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

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

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

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

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

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

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

 

3. Bifurcation diagram 2

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

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

 

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

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

   for p=linspace(0,1.7,1700);

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

       %等分成1700份

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

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

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

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

       %的向量.

       for k=1:500

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

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

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

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

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

           if k>400

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

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

               Z=[Z;x(1)];

               %这个值所对应最后400次迭代

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

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

           end

       end

   end

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

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

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

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

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

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

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

 

   d0=1e-8;%定义极小的扰动,还可以定义更小

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

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

   for p=linspace(0,1.7,1700)

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

       %等分成1700份

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

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

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

       x1=[1;d0;p;.5];%扰动初始值

       for k=1:1000

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

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

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

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

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

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

           d1=sqrt((x(1)-x1(1))^2+(x(2)-x1(2))^2);

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

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

           %选择基准轨线

           if k>500

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

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

               lsum=lsum+log(d1/d0);

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

           end

       end

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

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

       Z=[Z,le];%存储最后一个指数值

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

       end

   end

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

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

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

5 Chaos control to one of the fixed points

5.1 Fixed points

q=0.5;

for p=0:0.001:1.7

x1=1/(p-q+1);

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

 

5.2 Fixed points controller 

 

 

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)-1.222*(x(1)-1/2.2);

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

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

 

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

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

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

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

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

     X=feval(@Lozi,X);%