|
|
马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。
您需要 登录 才可以下载或查看，没有账号？我要加入
x
请问下高手：我的这程序怎么画不出分插图，望高手等指导下，多谢了！！！
function dy=Untitled0(t,y)
global Alpha;
global Beta;
global Deta;
global Gamma;
global Eta;
global v0;
global c0;
global c1;
global c2;
Alpha=1.2;Beta=1.3;Deta=0.8;Gamma=3.0;Eta=0.25;
c0=0.01;c1=0.002;c2=0.002;v0=[0.02:0.01:0.05];
if ((abs(y(1)+Alpha*(y(1)-y(3))+c0*(y(2)-y(4))+c1*y(2))<1.0&&abs(y(2)-v0)<=1.0e-4) && (abs(y(3)+Alpha*(y(3)-y(1))+c0*(y(4)-y(2))+c2*y(4))<Beta && abs(y(4)-v0)<=1.0e-4))
dy=zeros(4,1);
dy(1)=y(2);
dy(2)=0;
dy(3)=y(4);
dy(4)=0;
elseif (abs(y(1)+Alpha*(y(1)-y(3))+c0*(y(2)-y(4))+c1*y(2))>1.0||abs(y(2)-v0)~=1.0e-4)&&(abs(y(3)+Alpha*(y(3)-y(1))+c0*(y(4)-y(2))+c2*y(4))>Beta||abs(y(4)-v0)~=1.0e-4)
dy=zeros(4,1);
dy(1)=y(2);
dy(2)=-y(1)-Alpha*(y(1)-y(3))-c0*(y(2)-y(4))-c1*y(2)-sign(y(2)-v0)*((1.0-Deta)/(1.0+Gamma*fabs(y(2)-v0))+Deta+Eta*(y(2)-v0)^2);
dy(3)=y(4);
dy(4)=-y(3)-Alpha*(y(3)-y(1))-c0*(y(4)-y(2))-c2*y(4)-sign(y(4)-v0)*((1.0-Deta)/(1.0+Gamma*fabs(y(4)-v0))+Deta+Eta*(y(4)-v0)^2);
elseif (abs(y(1)+Alpha*(y(1)-y(3))+c0*(y(2)-y(4))+c1*y(2))<1.0&&abs(y(2)-v0)<=1.0e-4) && (abs(y(3)+Alpha*(y(3)-y(1))+c0*(y(4)-y(2))+c2*y(4))>Beta||abs(y(4)-v0)~=1.0e-4)
dy=zeros(4,1);
dy(1)=y(2);
dy(2)=0;
dy(3)=y(4);
dy(4)=-y(3)-Alpha*(y(3)-y(1))-c0*(y(4)-y(2))-c2*y(4)-sign(y(4)-v0)*((1.0-Deta)/(1.0+Gamma*fabs(y(4)-v0))+Deta+Eta*(y(4)-v0)^2);
elseif (abs(y(1)+Alpha*(y(1)-y(3))+c0*(y(2)-y(4))+c1*y(2))>1.0||abs(y(2)-v0)~=1.0e-4)&& (abs(y(3)+Alpha*(y(3)-y(1))+c0*(y(4)-y(2))+c2*y(4))<Beta&&abs(y(4)-v0)<=1.0e-4)
dy=zeros(4,1);
dy(1)=y(2);
dy(2)=-y(1)-Alpha*(y(1)-y(3))-c0*(y(2)-y(4))-c1*y(2)-sign(y(2)-v0)*((1.0-Deta)/(1.0+Gamma*fabs(y(2)-v0))+Deta+Eta*(y(2)-v0)^2);
dy(3)=y(4);
dy(4)=0;
end
clear;
global F ;
v0=[0.02:0.01:0.05];
k=0;
YY2=[];
for F=v0
y0=[0.01 0.01 0.01 0.01];
F;
k=k+1;
% discard the first 60 periodic data;
%除去前面60个周期的数据，并将最后的结果作为下一次积分的初值
tspan=[0:0.01:100];
[t,Y]=ode45(@Untitled0,tspan,y0);
y0=Y(end,:);
j=1;
for i=60:200
tspan=[0:0.01:100];
[t,Y]=ode45(@Untitled0,tspan,y0);
YY1(k,j)=Y(end,1); % get the omega data from every period end
j=j+1; %取出每一个周期内的第一个解的最后一个值。
y0=Y(end,:);
end
end
bifdata=YY1(:,end-51:end);
plot(v0,bifdata,'k.','markersize',1);
|
|
发表于 2010-11-16 15:38
提升卡
变色卡
显身卡
发表于 2010-11-17 16:39
楼主
