经过近期的研究发现,目前对于系统单参数分岔图的计算共有以下的几种方法: 1)最大值法 即对系统微分方程(组)进行求解,对求解的结果用getmax函数进行取点,并绘图。 2)Poincare截面法 对系统参数的每一次取值,绘制其Poincare截面,进而得到其分岔图。 这种方法需要注意的是,自治系统的Poincare截面是选取一超平面,平面上点的分布即构成一Poincare截面,非自治系统的Poincare截面则是根据系统激励的频率进行取点并绘图。 本帖将以Lorenz系统为例,对这两种方法进行比较 首先对第二种方法进行阐述。 编程如下(matlab) Lorenz系统: function dy = Lorenz(t,y) % Lorenz系统 % 系统微分方程: % dx/dt = -a(x-y) % dy/dt = x(r-z)-y % dz/dt = xy-bz % a=y(4) % r=y(5) % b=y(6) dy=zeros(6,1); dy(1)=-y(4)*(y(1)-y(2)); dy(2)=y(1)*(y(5)-y(3))-y(2); dy(3)=y(1)*y(2)-y(6)*y(3); dy(4)=0; dy(5)=0; dy(6)=0; 随r的分岔图求解程序:——按照x=y平面取截面 function Lorenz_bifur_r Z=[]; for r=linspace(1,500,1000); % 舍弃前面迭带的结果,用后面的结果画图 [T,Y]=ode45('Lorenz',[0,1],[1;1;1;16;r;4]); [T,Y]=ode45('Lorenz',[0,50],Y(length(Y),:)); Y(:,1)=Y(:,2)-Y(:,1); % 对计算结果进行判断,如果点满足x=y,则取点 for k=2:length(Y) f=k-1; if Y(k,1)<0 if Y(f,1)>0 y=Y(k,2)-Y(k,1)*(Y(f,2)-Y(k,2))/(Y(f,1)-Y(k,1)); Z=[Z r+abs(y)*i]; end else if Y(f,1)<0 y=Y(k,2)-Y(k,1)*(Y(f,2)-Y(k,2))/(Y(f,1)-Y(k,1)); Z=[Z r+abs(y)*i]; end end end end plot(Z,'.','markersize',1) title('Lorenz映射分岔图') xlabel('r'),ylabel('|y| where x=y') 结果如图1所示。 ![]() getmax法——取最大值法 function [Xmax] = getmax(y) a=length(y); j=1; for i=(a-1)/2:a b=(y(i,1)-y(i-2,1))/2; c=(y(i,1)+y(i-2,1))/2-y(i-1,1); if y(i-2,1)<=y(i-1,1)&y(i-1,1)>=y(i,1)&c==0 Xmax(j)=y(i-1,1); j=j+1; elseif y(i-2,1)<=y(i-1,1)&y(i-1,1)>=y(i,1) Xmax(j)=y(i-1,1)-b^2/(4*c); j=j+1; end end function Lorenz_bifur_r_getmax % 最大值法求解分岔图 clear all t0=[0 100];%积分时间 %bifurcation for r=linspace(1,500,1000); %r的变化精度 [t,y]=ode45('Lorenz',t0,[1;1;1;16;r;4]); [Xmax]=getmax(y(:,1)); plot(r,Xmax,'b','markersize',1) hold on clear Xmax end 计算结果如图2所示! ![]() 最后上传一下参考计算机仿真第22卷第12期上一篇文章“李雅普诺夫指数的研究与仿真”中Lorenz系统的分岔图计算结果,大家比较一下即可看出孰优孰劣了! ![]() 附件: |
GMT+8, 2025-4-3 18:44 , Processed in 0.090063 second(s), 26 queries , Gzip On.
Powered by Discuz! X3.4
Copyright © 2001-2021, Tencent Cloud.