Lozi混沌映射系统的Matlab实现
% 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(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=;
%前面两个赋初值,后面两个给出了系统中的
%常数项的值！与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=,分别指代每一次迭代后的
%数组X的上面两个变量,然后点画图.
title('Lozi映射图')%加上图像标题
xlabel('x'),ylabel('y')
%分别给对应的向量随处的坐标命名
%另外还可以使用其它的语句对图像进行进一步的调整
2.Bifurcation diagram 1
function newx=Lozi(x);
% Lozi方程[差分方程]%newx=Lozi();
% 方程如下:%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在区间
%等分成1700份
x=;
%前面两个赋初值,后面两个给出了系统中的
%常数项的值！与newx=Lozi(x)中给出的列
%向量相对应,其中值得注意的是p是连续变化
%的向量.
for k=1:500 %循环即迭代次数
x=Lozi(x); %调用主函数
if k>400
%取迭代后的最后一百个点,这里认为
%前400个点为迭代过程中的瞬态响应！
Z=;
%将P的值以及在这个值所对应最后400次迭代
%后的数组X的第一个变量x存储为一个二维列向量组！
end
end
end
plot(Z,'.','markersize',1)
%得到Z=,分别指代P的值以及在这个
%值所对应最后400次迭代后的数组X的第一个变量x,然后点画图.
title('Lozi映射分岔图')%加上图像标题
xlabel('p'),ylabel('x')
%分别给对应的向量随处的坐标命名
%另外还可以使用其它的语句对图像进行进一步的调整
3. Bifurcation diagram 2
function newx=Lozi(x);
% Lozi方程[差分方程]%newx=Lozi();
% 方程如下:%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在区间
%等分成1700份
x=;
%前面两个赋初值,后面两个给出了系统中的
%常数项的值！与newx=Lozi(x)中给出的列
%向量相对应,其中值得注意的是p是连续变化
%的向量.
for k=1:500
%循环即迭代次数,但是遗憾的是这里的迭代次数太少,
%因为如果太多的话,由于上面的p点太多,运行的时间
%就越长,但是为了得到较为精确的迭代解,建议适当的增加迭代次数
%估计在10000左右,而这时候取最后的1000个点即可！
x=Lozi(x); %调用主函数
if k>400
%取迭代后的最后一百个点,这里认为
%前400个点为迭代过程中的瞬态响应！
Z=;
%这个值所对应最后400次迭代
%后的数组X的第一个变量x存储为一个二维列向量组！
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(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在区间
%等分成1700份
le=0;%定义指数的初始值为0.
lsum=0;%定义和的初始值为0.
x=;%初始值
x1=;%扰动初始值
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=;%存储最后一个指数值
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(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=;
%前面两个赋初值,后面两个给出了系统中的
%常数项的值！与newx=Lozi(x)中给出的列向量相对应！
Y=[]; %给出一个空数组，以便存储得出的迭代解！
for i=1:10000 %循环即迭代次数
X=feval(@Lozi,X);%
附件:
http://forum.vibunion.com/data/attachment/album/2008/07/69372_200807191019561jJNL.thumb.jpg
Lozi Map[时间:2008-7-19 10:19]
http://forum.vibunion.com/data/attachment/album/2008/07/69372_2008071910245119sU8.thumb.jpg
bifurcation diagram 1[时间:2008-7-19 10:24]
http://forum.vibunion.com/data/attachment/album/2008/07/69372_200807191023531YtoH.thumb.jpg
bifurcation diagram 2[时间:2008-7-19 10:23]
http://forum.vibunion.com/data/attachment/album/2008/07/69372_200807191026331N678.thumb.jpg
Lyapunov-exponent spectrum[时间:2008-7-19 10:26]
http://forum.vibunion.com/data/attachment/album/2008/07/69372_200807191030261D0xO.thumb.jpg
control chaos to one of the fixed points[时间:2008-7-19 10:30]
% 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=,分别指代每一次迭代后的
%数组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*') 太好了，谢谢分享！
我在应用第4 Lyapunov-exponent spectrum于一个超混沌中时，出现错误提示：Warnning: Divide by zero.
不知道会是哪项不合适？ 有可能出现Warnning: Divide by zero，但是一般不影响！建议使用avoid warning: Divide by zero这种语句来实现！ 正在学习中，谢谢分享！！！ 谢谢分享。。。
页:
[1]