混沌動力學行為研究的程序說明

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('ouhe1',20,Y(length(Y),:));Y(:,1)=Y(:,2)-Y(:,1);for k=2:length(Y)f=k-1; ???if Y(k,1)<0 ?if Y(f,1)>0y=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)<0y=Y(k,2)-Y(k,1)*(Y(f,2)-Y(k,2))/(Y(f,1)-Y(k,1));Z=[Z a+abs(y)*i];endend ??endendplot(Z,'.','markersize',1);title('ouhe映射分岔圖');xlabel('a'),ylabel('|y|')

?

2.功率譜：對上面耦合系統的功率譜的研究

clear allx=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;for i=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([0 4 -3 6]); ?%限制橫、縱坐標范圍%plot(c,abs(X1(1:4096))) ????有時候也可用x的絕對值表示功率大小，沒求對數xlabel('\itf\rm/HZ','fontsize',18,'fontName','times new Roman','fontweight','bold','color','k'); %加橫坐標,\it表傾斜，\rm表復正ylabel('power spectrum/dB','fontsize',18,'fontName','times new Roman','fontweight','bold','color','k'); ?%縱坐標標示

3.最大李雅譜指數程序：仍對上述耦合系統

clear; ?%與ouhe系統的分岔圖（fotran）符合的很好clc;d0=1e-8;le=0;lsum=0;x=1;y=1;z=1;x1=1;y1=1;z1=1+d0;for i=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);if i>100lsum=lsum+log(d1/d0); ??endendle=lsum/(i-100)

最大李雅譜指數譜程序：仍對上述耦合系統

clear; %可調參數區間和步長，如a的第5和36行clc;LE1=[];d0=1e-8;for a=linspace(0.5,10.5,300);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+(d0/d1)*(y1-y);z1=z+(d0/d1)*(z1-z);if i>50lsum=lsum+log(d1/d0); ??endendle=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;global u aN1=linspace(0,0,200);N2=linspace(0,0,400);for I=1:200u=1.5+I*0.025;d0=1e-8;for L=1:400a=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=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 i>50lsum=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); %畫偽彩圖%colormap jet,shading 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*(-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 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);endlp=(log(diag(S))+b1)/(k*h);Z1=[Z1 lp(1)];Z2=[Z2 lp(2)];Z3=[Z3 lp(3)];enda=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 on

6.三維、二維以及時間序列圖：

耦合系統的函數程序：

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);dx(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);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','times 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','fontweight','bold','color','k');figure(2);plot3(x(:,1),x(:,2),x(:,3),'k','markersize',0.5);grid on;xlabel('\itx\rm_1','fontsize',20,'fontName','times new Roman','fontweight','bold','color','k');ylabel('\itx\rm_2','fontsize',20,'fontName','times new Roman','fontweight','bold','color','k');zlabel('\itx\rm_3','fontsize',20,'fontName','times new Roman','fontweight','bold','color','k');

取參數a=3:0.001:4,迭代初值x(0)=0.1，計算Logistic映射的Lyapunov指數，并做可視化呈現

clear all; hold on alpha0=3:0.001:4; N=1000; for j=1:length(alpha0 alpha=alpha0(j); x0=0.1;%初始值 s=0; for ii=1:N df=alpha-2*alpha*x0; s=s+log(abs(df));%lambda疊加 x0=alpha*x0*(1-x0);%x迭代 end Lm(j)=s/N;% 指數 end plot(alpha0,Lm,'r')

N維離散系統動力學Mapping的分岔圖與李雅普諾夫指數計算

clear all;clca=4;b1=1;b2=1;d=0.5;c1=1.08;c2=1;z=0.5;m=0.25;p=0.5;T=0.1;S=0.1;k2=0.15;N1 = 100;N2 = 400;k10 = 0:0.001:0.4;f = zeros(N1+N2,length(k10));f2 = zeros(N1+N2,length(k10));for kk=1:length(k10)x0=2.5;y0=2.35;k1=k10(kk)L1=0;L2=0;% ????s=zeros(2,1);J1=eye(2);for j = 1:N1+N2x1=x0+ k1.*x0.*(a-2*b1.*x0+d.*y0+b1.*(c1+(1-m).*z-S));y1=y0+ k2.*y0.*(a-2*b2.*y0+d.*x0+b2.*(c2-p.*m.*z+T));f(j,kk) =x1;f2(j,kk) =y1;j1=1 + a*k1 + d*k1*y1 + b1*k1 *(c1 - S - 4*x1 + z - m*z);j2=d*k1*x1;j3=d*k2*y1;j4=1 + a*k2 + d *k2 *x1 + b2 *k2 *(c2 + T - 4 *y1 - m *p* z);x0=x1;y0=y1;endendhold onf = f(N1+1:end,:);plot(k10,f,'b.','MarkerSize',1)ylabel('Price');hold onf2 = f2(N1+1:end,:);plot(k10,f2,'r.','MarkerSize',1)xlabel('\mu');ylabel('Price');% plot(k10,Lm2,'r','linewidth',0.5)plot(k10,Lm1,'b')line([0 0.4],[0 0])

=0 臨界點

<0穩定 >0不穩定，分叉

?

總結

以上是生活随笔為你收集整理的混沌动力学行为研究-分叉图的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。