连续时间混沌系统MATLAB程序和SIMULINK模型.doc

1、第6章 连续时间混沌系统本章讨论连续时间混沌系统的基本特点与分析方法，主要包括混沌数值仿真和硬件实验方法简介、混沌系数平衡点的计算、平衡点的分类与性质、相空间中的轨道、几类典型连续混沌系统的介绍、混沌机理的分析方法、用特征向量空间法寻找异宿轨道、Lorenz系统及混沌机理定性分析、Lorenz映射、Poincare截面、Chua系统及其混沌机理定性分析、时间序列与相空间重构等内容。6.1 混沌数值仿真和硬件实验方法简介混沌的数值仿真主要包括MATLAB编程、SIMULINK模块构建、EWB仿真以及其他一些相关的软件仿真或数值计算等方法，从而获取混沌吸引子的相图、时域波形图、李氏指数、分叉图和功

2、率谱等。混沌的硬件实验主要包括模拟/数字电路设计与硬件实验、现场可编程门阵列器件（FPGA）、数字信号处理器（DSP）等硬件实现方法来产生混沌信号。本节仅对各种数值仿真方法作简单介绍。1） 混沌系统的MATLAB数值仿真该方法主要根据混沌系统的状态方程来编写MATLAB程序。现举二例来说明这种编程方法。（1） 已知Lorenz系统的状态方程为dx/dt=-a(x-y)dy/dt=bx-xz-ydz/dt=-cz+xy式中a=10,b=30,c=8/3。MATLAB仿真程序如下： %*Function dxdt=lorenz(t,x) %除符号dxdt外，还可用其他编程者习惯的有意义的符号A=1

3、0;B=30;C=8/3;dxdt=zeros(3,1);dxdt(1)=-A*(x(1)-x(2);dxdt(2)=B*x(1)-x(1).*x(3)-x(2);dxdt(3)=x(1)*x(2)-C*x(3);%*options=odeset(RelTol,1e-6,AbsTol, 1e-6 1e-6 1e-6);t0=0 200;x0=0.02,0.01,0.03;t,x=ode45(lorenz,t0,x0,options);%*n=length(t)n1=round(n/2)%n1=1;%*figure(1);plot(t(n1:n,1),x(n1:n,1);xlabel(t,fon

4、tsize,20,fontname,times new roman,FontAngle,italic);ylabel(x,fontsize,20,fontname,times new roman,FontAngle,italic);figure(2);plot(x(n1:n,1),x(n1:n,3);xlabel(x,fontsize,20,fontname,times new roman,FontAngle,italic);ylabel(z,fontsize,20,fontname,times new roman,FontAngle,italic);%*根据上述MATLAB程序，得Loren

5、z系统的时域波形图和混沌吸引子相图的数值仿真结果如图6-1所示。 图6-1 Lorenz系统的时域波形图和混沌吸引子相图的MATLAB数值仿真结果（2）已知Chua系统的状态方程为dx=ay-f(x)dy=x-y+zdz=-by式中a=10,b=15,m0=-1/7,m1=2/7,f(x)=m1*x+0.5(m0-m1)|x+1|-|x-1|为三分段非线性函数MATLAB仿真如下：function dxdt=chua(t,x)m0=-1/7;m1=2/7;a=10;b=15;%*dxdt=zeros(3,1);f=m1*x(1)+0.5*(m0-m1)*(abs(x(1)+1)-abs(x(1

6、)-1);dxdt(1)=a*(x(2)-f);dxdt(2)=x(1)-x(2)+x(3);dxdt(3)=-b*x(2);%*options=odeset(RelTol,1e-6,AbsTol,1e-6 1e-6 1e-6);t0=0 5e+2;x0=0.01 0.02 0.03;t,x=ode45(chua,t0,x0,options);%*n=length(t)n1=round(n/2)%*figure(1);plot(t(n1:n),x(n1:n,1);xlabel(t,fontsize,20,fontname,times new roman,FontAngle,italic);yl

7、abel(x,fontsize,20,fontname,times new roman,FontAngle,italic);figure(2);plot(x(n1:n,1),x(n1:n,2);xlabel(x,fontsize,20,fontname,times new roman,FontAngle,italic);ylabel(y,fontsize,20,fontname,times new roman,FontAngle,italic);%*根据上述MATLAB程序，得Chua系统的时域波形图和混沌吸引子相图的数值仿真结果如图6-2所示。 图6-2 Chua系统的时域波形图和混沌吸引子

8、相图的MATLAB数值仿真结果2） 混沌系统的SIMULINK仿真该方法主要是根据混沌系统的状态方程，将其转换成积分方程，利用模块和积分算子画出SIMULINK的模块化仿真图。为保证计算的精确度，又不使仿真时间过长，应对仿真图中几个重要参数进行设置。第一个参数是仿真时间：第二个参数是相对误差，通常设为；第三个参数是绝对误差，通常设为，现举二例来说明这种编程方法。（1） 已知Lorenz系统的状态方程仍如（6-1）式，将其转换成积分方程：dx/dt=-a(x-y) dy/dt=bx-xz-y dz/dt=-cz+xy 注意 ，SIMULINK仿真中的微分子算子为S，积分算子为,故得SIMULIN

9、K的仿真如图6-3所示，设其文件名为“simulink_lorenz”,再利用文件名为“y_simulink_lorenz”的程序运行“simulink_lorenz”。程序如下：t,x=sim(simulink_lorenz,200);n=length(t)n1=round(n/2)%*figure(1);plot(t(n1:n,1),x(n1:n,1);xlabel(t,fontsize,20,fontname,times new roman,FontAngle,italic);ylabel(x,fontsize,20,fontname,times new roman,FontAngle,

10、italic);figure(2);plot(x(n1:n,1),x(n1:n,3);xlabel(x,fontsize,20,fontname,times new roman,FontAngle,italic);ylabel(z,fontsize,20,fontname,times new roman,FontAngle,italic);%*运行结果仍如图6-1所示。图6-3 Lorenz系统的SIMULINK仿真图（2）已知Chua系统的状态方程仍如（6-2）式，得SIMULINK的仿真图如图6-4所示。设其文件名为“simulink_chua”,再利用文件名为“y_simulink_ch

11、ua”的程序运行“simulink_chua”。程序如下：%*%global a;a=10;t,x=sim(simulink_chua,500);n=length(t)n1=round(n/2)%*figure(1);plot(t(n1:n,1),x(n1:n,1);xlabel(t,fontsize,20,fontname,times new roman,FontAngel,italic);ylabel(x,fontsize,20,fontname,times new roman,FontAngel,italic);figure(2);plot(x(n1:n,1),x(n1:n,2);xla

12、bel(x,fontsize,20,fontname,times new roman,FontAngel,italic);ylabel(y,fontsize,20,fontname,times new roman,FontAngel,italic);运行结果仍如图6-2所示。 （a）SIMULINK主框图 （b）三分段线性函数f的SIMULINK子框图图6-4 Chua系统的SIMULINK仿真图3） 连续混沌系统离散化的MATLAB数值仿真当用DSP和FPGA等现代数字器件来产生混沌信号时，首先需要将连续混沌系统作离散化处理。离散化和数字化处理方法主要有三种，利用这些离散化的方法，可将状态方

13、程变成差分方程，这些方法将在后续章节中详细介绍。这里采用了一种较简单的Euler算法。（1）已知Lorenz系统的状态方程仍如（6-1）式，根据Euler算法得对应的差分方程为 将其整理成标准的差分方程形式： 根据（6-3）式，得MATLAB仿真程序如下：clear all;hold off;%*x(1)=0.02;y(1)=0.01;z(1)=0.03;a=10;b=30;c=8/3;T=5e-3;%可调取样时间N=5e4;%可调迭代次数for i=1:Nit(i)=i;x(i+1)=(1-a*T)*x(i)+a*T*y(i);y(i+1)=T*b*x(i)-T*x(i)*z(i) +(1-

14、T)*y(i);z(i+1)=(1-T*c)*z(i)+T*x(i)*y(i);end%*n1=N/2;n2=N;figure(1)plot(t(n1:n2),x(n1:n2),b-);xlabel(t_n,fontsize,20,fontname,times new roman,FontAngle,italic);ylabel(x_n,fontsize,20,fontname,times new roman,FontAngle,italic);figure(2);plot(x(n1:n2),x(n1:n2),b-);xlabel(x_n,fontsize,20,fontname,times

15、new roman,FontAngle,italic);ylabel(z_n,fontsize,20,fontname,times new roman,FontAngle,italic);%*运行结果仍如图6-1所示。(2)已知Chua系统的状态方程仍如（6-2）式，根据Euler算法，得对应的差分方程为 将其整理成标准的差分方程形式： 式中a=10,b=15,f(x(n)=m1x(n)+0.5(m0-m1)|x(n)+1|-|x(n)-1|为离散化后的三分段线性函数，其中参数m0=-1/7,m1=2/7.根据（6-4）式，得MATLAB仿真程序如下： clear all; hold off;

16、 %* x(1)=0.02;y(1)=0.01;z(1)=0.03;m0=-1/7;m1=2/7;a=10;b=15;T=1e-2;%N=5e4;%for i=1:Nit(i)=i;f(i)=m1*x(i)+0.5*(m0-m1).*(abs(x(i)+1)-abs(x(i)-1);x(i+1)= x(i)+a*T*y(i)-T*a*f(i);y(i+1)=T *x(i)+(1-T)*y(i)+T*z(i);z(i+1)= -T*b*y(i)+z(i);end%*n1=N/2;n2=N;figure(1)plot(t(n1:n2),x(n1:n2),b-);xlabel(t_n,fontsiz

17、e,20,fontname,times new roman,FontAngle,italic);ylabel(x_n,fontsize,20,fontname,times new roman,FontAngle,italic);figure(2);plot(x(n1:n2),y(n1:n2),b-);xlabel(x_n,fontsize,20,fontname,times new roman,FontAngle,italic);ylabel(y_n,fontsize,20,fontname,times new roman,FontAngle,italic);%*运行结果仍如图6-2所示。4）

18、 离散混沌系统的SIMULINK仿真9仍以Lorenz系统和Chua系统为例来说明。(1)Lorenz系统离散化后的SIMULINK仿真。为方便计，将（6-3）式表示为如下的标准形式： 式中A=1-Ta,B=Ta,C=Tb,D=T,E=1-T,F=1-Tc,G=T,得离散化后Lorenz系统的SIMULINK仿真如图6-5所示。注意，由于是离散系统，故应将连续系统中的算子1/S换成离散系统中的算子1/Z。设其文件为“simulink_discrete_lorenz”的程序运行“simulink_discrete_lorenz”。程序如下：图6-5 离散化后Lorenz系统的SIMULINK仿真

19、图%*a=10;b=30;c=8/3;T=5e-3;%可调取样时间%global A;global B;global D;global E;global F;global G;A=1-T*a;B=T*a;C=T*b;D=T;E=1-T*c;G=T;%*t,x=sim(simulink_discrete_lorenz,5e4);n=length(t)n1=round(n/2);%*figure(1);plot(t(n1:n,1),x(n1:n,1);xlabel(t,fontsize,20,fontname,times new roman,FontAngel,italic);ylabel(x,f

20、ontsize,20,fontname,times new roman,FontAngel,italic);figure(2);plot(x(n1:n,1),x(n1:n,2);xlabel(x,fontsize,20,fontname,times new roman,FontAngel,italic);ylabel(z,fontsize,20,fontname,times new roman,FontAngel,italic);%*运行结果仍如图6-1所示。（2）Chua系统离散化后的SIMULINK仿真。为方便计，将（6-4）式表示为如下的标准形式： 式中，A=1,B=Ta,C=Ta,D=

21、T,E=1-T,F=T,G=Tb,H=1,T=1e-2,a=10,b=15,m0= -1/7,m1=2/7,得离散化后 Lorenz系统的SIMULINK仿真图如6-6所示。(a)SIMULINK主框图 （b）三分段线性函数f(x(n)的SIMULINK子框图图6-6 离散化后Chua系统的SIMULINK仿真图设其文件名为“simulink_discrete_chua”,利用文件名为“y_simulink_discrete_chua”的程序运行“simulink_discrete_chua”。程序如下：%*a=10;b=15;T=1e-2;%可调取样时间%global A;global B;

22、global C;global D;%global E;global F;global G;global H;A=1;B=T*a;C=T*a;D=T;E=1-T;G=T*b;H=1;%*t,x=sim(simulink_discrete_lorenz,5e4);n=length(t)n1=round(n/2);%*figure(1);plot(t(n1:n,1),x(n1:n,1);xlabel(t,fontsize,20,fontname,times new roman,FontAngel,italic);ylabel(x,fontsize,20,fontname,times new rom

23、an,FontAngel,italic);figure(2);plot(x(n1:n,1),x(n1:n,2);xlabel(x,fontsize,20,fontname,times new roman,FontAngel,italic);ylabel(z,fontsize,20,fontname,times new roman,FontAngel,italic);%*运行结果仍如图6-2所示。以上我们分别以Lorenz系统和chua系统为典型实例，对四种数值仿真方法进行了详细的说明。这四种仿真方法分别是混沌系统的MATLAB数值仿真、混沌系统的SIMULINK仿真、连续混沌系统离散化的MATLAB数值仿真和离散混沌系统的SIMULINK仿真。有关用模拟电路、DSP、FPGA产生混沌的问题将在后续章节中给出，此处不再讨论。