分岔图的总结

混沌动力学行为研究的程序说明1.分岔图：一个耦合发电机系统：112322133312()()xxxxaxxxxaxxxx耦合系统函数：functiondx=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=[];fora=linspace(0.5,10.5,500);%舍弃前面迭带的结果，用后面的结果画图：即0.5-10.5分为500点[T,Y]=ode45('ouhe1',1,[1;1;1;2;a]);[T,Y]=ode45('ouhe1',20,Y(length(Y),:));Y(:,1)=Y(:,2)-Y(:,1);fork=2:length(Y)f=k-1;ifY(k,1)0ifY(f,1)0y=Y(k,2)-Y(k,1)*(Y(f,2)-Y(k,2))/(Y(f,1)-Y(k,1));Z=[Za+abs(y)*i];endelseifY(f,1)0y=Y(k,2)-Y(k,1)*(Y(f,2)-Y(k,2))/(Y(f,1)-Y(k,1));Z=[Za+abs(y)*i];endendendendplot(Z,'.','markersize',1);title('ouhe映射分岔图');xlabel('a'),ylabel('|y|')0246810-6-4-2024681012x1a2.功率谱：对上面耦合系统的功率谱的研究clearallx=zeros(1,10001);y=zeros(1,10001);z=zeros(1,10001);%数组置零x(1)=1;y(1)=1;z(1)=1;h=0.001;k=10000;a=3;u=2;fori=1:kx(i+1)=x(i)+h*(-u*x(i)+y(i)*(z(i)+a));%欧拉离散y(i+1)=y(i)+h*(-u*y(i)+x(i)*(z(i)-a));z(i+1)=z(i)+h*(z(i)-x(i)*y(i));endX1=fft(z,16384);%对x做傅里叶变换，取8192个点p=X1.*conj(X1)/16384;%求x的模及功率谱密度，单位：dB同样可求y或zc=100*[0:8191]/16384;%取双边，也可取单边c=[0:4095]/0.8192;%title('题目')%顶端题目plot(c,log10(p(1:8192)),'k')%画出左半部分axis([04-36]);%限制横、纵坐标范围%plot(c,abs(X1(1:4096)))有时候也可用x的绝对值表示功率大小，没求对数xlabel('\itf\rm/HZ','fontsize',18,'fontName','timesnewRoman','fontweight','bold','color','k');%加横坐标,\it表倾斜，\rm表复正ylabel('powerspectrum/dB','fontsize',18,'fontName','timesnewRoman','fontweight','bold','color','k');%纵坐标标示05101520253035404550-5-4-3-2-101234f/HZpowerspectrum/dB3.最大李雅谱指数程序：仍对上述耦合系统clear;%与ouhe系统的分岔图（fotran）符合的很好clc;d0=1e-8;le=0;lsum=0;x=1;y=1;z=1;x1=1;y1=1;z1=1+d0;fori=1:500[T1,Y1]=ode45('ouhe1',[0,1],[x;y;z;4.05;7.0]);[T2,Y2]=ode45('ouhe1',[0,1],[x1;y1;z1;4.05;7.0]);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+(d0/d1)*(y1-y);z1=z+(d0/d1)*(z1-z);ifi100lsum=lsum+log(d1/d0);endendle=lsum/(i-100)最大李雅谱指数谱程序：仍对上述耦合系统clear;%可调参数区间和步长，如a的第5和36行clc;LE1=[];d0=1e-8;fora=linspace(0.5,10.5,300);le=0;lsum=0;x=1;y=1;z=1;x1=1;y1=1;z1=1+d0;fori=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+(d0/d1)*(y1-y);z1=z+(d0/d1)*(z1-z);ifi50lsum=lsum+log(d1/d0);endendle=lsum/(i-50);LE1=[LE1le];enda=linspace(0.5,10.5,300);plot(a,LE1,'-');title('largestLyapunovexponentsofouhe1');xlabel('parametera'),ylabel('largestLyapunovexponents');grid024681012-0.200.20.40.60.811.2largestLyapunovexponentsofouhe1parameteralargestLyapunovexponents双参数空间的最大李雅谱指数谱程序：仍对以上耦合系统clear;%可调参数区间和步长clc;globaluaN1=linspace(0,0,200);N2=linspace(0,0,400);forI=1:200u=1.5+I*0.025;d0=1e-8;forL=1:400a=0.5+L*0.025;le=0;lsum=0;x=1;y=1;z=1;x1=1;y1=1;z1=1+d0;fori=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=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);ifi50lsum=lsum+log(d1/d0);endendle=lsum/(i-50);LE1(I,L)=le;N2(L)=a;endN1(I)=u;end[X,Y]=meshgrid(N1,N2);Z=LE1;pcolor(X,Y,Z);%画伪彩图%colormapjet,shadinginterp%连续变化的变异饱和色图，表面画伪彩图%contourf(X,Y,Z)%画等高线title('largestLyapunovexponentsofouhe1')xlabel('parameter\itu')ylabel('parameter\ita')zlabel('最大李雅普指数{\delta}','FontSize',12)234567246810ua01E-30.20000.30000.40000.50000.60000.70000.80000.90001.0001.1001.2004.李雅谱指数程序：仍对以上耦合系统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;fori=1:kx1=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=[-ua+zyz-a-ux-y-x1];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;fora=linspace(0.5,10.5,1000);V=eye(3);S=V;b1=0;lp=0;fori=1:kx1=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=[-ua+zyz-a-ux-y-x1];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);endlp=(log(diag(S))+b1)/(k*h);Z1=[Z1lp(1)];Z2=[Z2lp(2)];Z3=[Z3lp(3)];enda=linspace(0.5,10.5,1000);plot(a,Z1,'-',a,Z2,'-',a,Z3,'-');title('Lyapunovexponentsofouhe');xlabel('parametera'),ylabel('lyapunovexponents');gridon024681012-5-4-3-2-1012Lyapunovexponentsofouheparameteralyapunovexponents6.三维、二维以及时间序列图：耦合系统的函数程序：functiondx=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);dx(3)=x(3)-x(1)*x(2);图程序：clc;clear;[t,x]=ode45('ouhe',[01000],[1-11]);subplot(2,2,1);plot(x(:,1),'k','markersize',0.5);xlabel('t/ms','fontsize',12,'fontName','timesnewRoman','fontweight','bold','color','k');ylabel('x','fontsize',12,'fontName','timesnewRoman','fontweight','bold','color','k');subplot(2,2,2);plot(x(:,2),'k','markersize',0.5);xlabel('t/ms','fontsize',12,'fontName','timesnewRoman','fontwe