《分岔图的总结.doc》由会员分享,可在线阅读,更多相关《分岔图的总结.doc(15页珍藏版)》请在得力文库 - 分享文档赚钱的网站上搜索。
1、混沌动力学行为研究的程序说明1.分岔图:一个耦合发电机系统:耦合系统函数:function dx=ouhe1(t,x)dx(1,1)=-x(4)*x(1)+x(2)*(x(3)+x(5);dx(2,1)=-x(4)*x(2)+x(1)*(x(3)-x(5);dx(3,1)=x(3)-x(1)*x(2);dx(4,1)=0;dx(5,1)=0;分岔图程序:clear;clc;Z=;for a=linspace(0.5,10.5,500); %舍弃前面迭带的结果,用后面的结果画图:即0.5-10.5分为500点 T,Y=ode45(ouhe1,1,1;1;1;2;a); T,Y=ode45(ouh
2、e1,20,Y(length(Y),:); Y(:,1)=Y(:,2)-Y(:,1); for k=2:length(Y) f=k-1; if Y(k,1)0 y=Y(k,2)-Y(k,1)*(Y(f,2)-Y(k,2)/(Y(f,1)-Y(k,1); Z=Z a+abs(y)*i; end else if Y(f,1)100 lsum=lsum+log(d1/d0); endend le=lsum/(i-100)最大李雅谱指数谱程序:仍对上述耦合系统clear; %可调参数区间和步长,如a的第5和36行clc;LE1=;d0=1e-8;for a=linspace(0.5,10.5,300)
3、; le=0; lsum=0; x=1; y=1; z=1; x1=1; y1=1; z1=1+d0; for i=1:150 T1,Y1=ode45(ouhe1,0,1,x;y;z;2;a); T2,Y2=ode45(ouhe1,0,1,x1;y1;z1;2;a); n1=length(Y1); n2=length(Y2); x=Y1(n1,1); y=Y1(n1,2); z=Y1(n1,3); x1=Y2(n2,1); y1=Y2(n2,2); z1=Y2(n2,3); d1=sqrt(x-x1)2+(y-y1)2+(z-z1)2); x1=x+(d0/d1)*(x1-x); y1=y+(
4、d0/d1)*(y1-y); z1=z+(d0/d1)*(z1-z); if i50 lsum=lsum+log(d1/d0); end end le=lsum/(i-50); LE1=LE1 le;end a=linspace(0.5,10.5,300); plot(a,LE1,-); title(largest Lyapunov exponents of ouhe1); xlabel(parameter a),ylabel(largest Lyapunov exponents); grid 双参数空间的最大李雅谱指数谱程序:仍对以上耦合系统clear; %可调参数区间和步长clc;glob
5、al u aN1=linspace(0,0,200);N2=linspace(0,0,400);for I=1:200 u=1.5+I*0.025; d0=1e-8; for L=1:400 a=0.5+L*0.025; le=0; lsum=0; x=1; y=1; z=1; x1=1; y1=1; z1=1+d0; for i=1:150 T1,Y1=ode45(ouhe1,0,1,x;y;z;u;a); T2,Y2=ode45(ouhe1,0,1,x1;y1;z1;u;a); n1=length(Y1); n2=length(Y2); x=Y1(n1,1); y=Y1(n1,2); z=
6、Y1(n1,3); x1=Y2(n2,1); y1=Y2(n2,2); z1=Y2(n2,3); d1=sqrt(x-x1)2+(y-y1)2+(z-z1)2); x1=x+(d0/d1)*(x1-x); y1=y+(d0/d1)*(y1-y); z1=z+(d0/d1)*(z1-z); if i50 lsum=lsum+log(d1/d0); end end le=lsum/(i-50); LE1(I,L)=le; N2(L)=a; end N1(I)=u;end X,Y=meshgrid(N1,N2); Z=LE1; pcolor(X,Y,Z); %画伪彩图 %colormap jet,s
7、hading interp %连续变化的变异饱和色图,表面画伪彩图 %contourf(X,Y,Z)%画等高线 title(largest Lyapunov exponents of ouhe1) xlabel(parameter itu) ylabel(parameter ita) zlabel(最大李雅普指数delta,FontSize,12)4.李雅谱指数程序:仍对以上耦合系统clear;clc;x=1;y=1;z=1;h=0.003;a=3.5;u=2;V=eye(3);S=V;b1=0;k=50000;for i=1:k x1=x+h*(-u*x+y*(z+a); y1=y+h*(-
8、u*y+x*(z-a); z1=z+h*(z-x*y); x=x1;y=y1;z=z1; J=-u a+z y z-a -u x -y -x 1; J=eye(3)+h*J; B=J*V*S; V,S,U=svd(B); a_max=max(diag(S); S=(1/a_max)*S; b1=b1+log(a_max);endLyapunov=(log(diag(S)+b1)/(k*h)李雅谱指数谱程序:仍对以上耦合系统clear;%奇异值分解法计算ouhe系统的李雅普诺夫指数谱clc;Z1=;Z2=;Z3=;x=1;y=1;z=1;h=0.002;u=2;%a=3;k=10000;for
9、a=linspace(0.5,10.5,1000); V=eye(3); S=V; b1=0; lp=0; for i=1:k x1=x+h*(-u*x+y*(z+a); y1=y+h*(-u*y+x*(z-a); z1=z+h*(z-x*y); x=x1;y=y1;z=z1; J=-u a+z y z-a -u x -y -x 1; J=eye(3)+h*J; B=J*V*S; V,S,U=svd(B); a_max=max(diag(S); S=(1/a_max)*S; b1=b1+log(a_max); end lp=(log(diag(S)+b1)/(k*h); Z1=Z1 lp(1)
10、; Z2=Z2 lp(2); Z3=Z3 lp(3);end a=linspace(0.5,10.5,1000); plot(a,Z1,-,a,Z2,-,a,Z3,-); title(Lyapunov exponents of ouhe); xlabel(parameter a),ylabel(lyapunov exponents);grid on6.三维、二维以及时间序列图:耦合系统的函数程序:function dx=ouhe(t,x)dx=zeros(3,1);a=3;u=2;dx(1)=-u*x(1)+x(2)*(x(3)+a);dx(2)=-u*x(2)+x(1)*(x(3)-a);d
11、x(3)=x(3)-x(1)*x(2);图程序:clc;clear;t,x=ode45(ouhe,0 1000,1 -1 1);subplot(2,2,1); plot(x(:,1),k,markersize,0.5);xlabel(t/ms,fontsize,12,fontName,times new Roman,fontweight,bold,color,k);ylabel(x,fontsize,12,fontName,times new Roman,fontweight,bold,color,k);subplot(2,2,2); plot(x(:,2),k,markersize,0.5)
12、;xlabel(t/ms,fontsize,12,fontName,times new Roman,fontweight,bold,color,k);ylabel(y,fontsize,12,fontName,times new Roman,fontweight,bold,color,k);subplot(2,2,3); plot(x(:,3),k,markersize,0.5);xlabel(t/ms,fontsize,12,fontName,times new Roman,fontweight,bold,color,k);ylabel(z,fontsize,12,fontName,time
13、s new Roman,fontweight,bold,color,k);subplot(2,2,4); plot3(x(:,1),x(:,2),x(:,3),k,markersize,0.5);grid on;xlabel(x,fontsize,12,fontName,times new Roman,fontweight,bold,color,k);ylabel(y,fontsize,12,fontName,times new Roman,fontweight,bold,color,k);zlabel(z,fontsize,12,fontName,times new Roman,fontwe
14、ight,bold,color,k);figure(2);plot3(x(:,1),x(:,2),x(:,3),k,markersize,0.5);grid on;xlabel(itxrm_1,fontsize,20,fontName,times new Roman,fontweight,bold,color,k);ylabel(itxrm_2,fontsize,20,fontName,times new Roman,fontweight,bold,color,k);zlabel(itxrm_3,fontsize,20,fontName,times new Roman,fontweight,bold,color,k);第 15 页