Lozi混沌映射系统的Matlab实现

liliangbiao

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


附件:

                               
登录/注册后可看大图

Lozi Map  [时间:2008-7-19 10:19]

                               
登录/注册后可看大图

bifurcation diagram 1  [时间:2008-7-19 10:24]

                               
登录/注册后可看大图

bifurcation diagram 2  [时间:2008-7-19 10:23]

                               
登录/注册后可看大图

Lyapunov-exponent spectrum  [时间:2008-7-19 10:26]

                               
登录/注册后可看大图

control chaos to one of the fixed points  [时间:2008-7-19 10:30]

liliangbiao

% for the length limit of the journal in my space, the remains of the last Matlab codes are:
% Continue with the X=feval(@Lozi,X);%
X=feval(@Lozi,X);%调用主函数
     Y(i, : )=X(1:2,1);
     %将每一次迭代后的数组X的上面两个变量存储！
     %因为下面的两个常数项进行迭代时保持原来的数值不便！
end
     plot(Y(:,1),Y(:,2))%,'.','markersize',5);
     %得到Y=[Y(:,1),Y(:,2)],分别指代每一次迭代后的
     %数组X的上面两个变量,然后点画图.
     title('Lozi映射图')%加上图像标题
     xlabel('x'),ylabel('y')
     %分别给对应的向量随处的坐标命名
%另外还可以使用其它的语句对图像进行进一步的调整
>> hold on
>> plot(1/2.2,0.5/(1.7-(0.5-1)),'r*')
>> plot(1,0.1,'k*')

happy_7

太好了，谢谢分享！
我在应用第4 Lyapunov-exponent spectrum于一个超混沌中时，出现错误提示：Warnning: Divide by zero.
不知道会是哪项不合适？

liliangbiao

有可能出现Warnning: Divide by zero，但是一般不影响！建议使用avoid warning: Divide by zero这种语句来实现！

kingsir

正在学习中，谢谢分享！！！

kddhhm

谢谢分享。。。