|
楼主 |
发表于 2011-10-15 12:32
|
显示全部楼层
回复 19 # freedignity 的帖子
哈哈 现在的程序是这个样的:
function dx=duf(t,X)
global F wd;
r=0.01;
x=X(1);
y=X(2);
psi=X(3);
dx=zeros(3,1);
dx(1)=y;
dx(2)=-r*y+x*(1-x^2)+F*sin(psi);
dx(3)=wd;
clear;
global F wd;
wd=1.0;
range=[0.1:0.0001:0.3];
period=2*pi/wd; %
k=0;
YY2=[];
step=2*pi/100; %步长。
for F=range
y0=[0 0 0];
F
k=k+1;
% discard the first 60 periodic data;
%除去前面60个周期的数据,并将最后的结果作为下一次积分的初值
tspan=[0:step:60*period];
[t,Y]=ode45(@duf,tspan,y0);
y0=Y(end,:);
j=1;
for i=60:200
tspan=[i*period:step:(i+1)*period];
[t,Y]=ode45(@duf,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(range,bifdata,'k.','markersize',1);
我讨论的是矩形薄板,呵呵,把板长及厚度变了,mu就大了。
谢谢你这么耐心的解答我的疑问,共同进步吧
|
|