% 连续时间系统混沌化切控制方法
% 《动力系统的混沌化》陈关荣 汪小帆
%
% Example(函数图象):
% [T,Y]=ode45('chaos2',40,[2;2;1]);
% plot3(Y(:,1),Y(:,2),Y(:,3));
% hold on
% [T,Y]=ode45('chaos2',40,[2;2;-1]);
% plot3(Y(:,1),Y(:,2),Y(:,3));
%
function dx=chaos2(t,x);
a=3;b=20;c=-20;
k=4;d=10;
m=4;e=-10;
A= [ a b 0
-b a 0
0 0 c];
if x(3)+norm(x(1:2))>k & x(3)>0
u=k*[-x(1) -x(2) d]';
elseif x(3)-norm(x(1:2))<-m & x(3)<0
u=m*[-x(1) -x(2) e]';
else
u=0;
end
dx=A*x+u;
% 连续时间系统混沌化切控制方法
% 《动力系统的混沌化》陈关荣 汪小帆
%
% Example(函数图象):
% [T,Y]=ode45('chaos1',40,[2;2;1]);
% plot3(Y(:,1),Y(:,2),Y(:,3));
%
function dx=chaos1(t,x);
a=3;
%global a;
b=20;c=-20;
k=4;d=10;
A= [ a b 0
-b a 0
0 0 c];
if x(3)+norm(x(1:2))>k
u=k*[-x(1) -x(2) d]';
else
u=0;
end
dx=A*x+u;
% 求最大Lyapunov指数谱
% 连续时间系统混沌化切控制方法
% 《动力系统的混沌化》陈关荣 汪小帆
%
Z=[];
global a;
d0=1e-8;
for a=linspace(0.1,6,150);
y0=[2;2;1];
y=[2+d0;2;1];
lsum=0;
p=0;
for k=1:60
[T,Y0] = ode45('chaos1',.4,y0);
[T,Y] = ode45('chaos1',.4,y);
y0 = Y0(size(Y0,1),:);
y = Y(size(Y,1),:);
d1=norm(y-y0);
if d1~=0
y=y0+(d0/d1)*(y-y0);
if k>10
lsum=lsum+log(d1/d0);
end
else
p=p+1;
y(1)=y0(1)+d0;
end
end
le=lsum/(k-10-p);
Z=[Z,a+le*i];
end
plot(Z,'-')
title('Lyapunov exponents')
xlabel('parameter a'),ylabel('Maximum lyapunov exponents')
grid on
