请教一个分岔程序
本帖最后由 relou 于 2011-7-12 21:38 编辑function s=rk45(t,q,flag,w)
%global m c k ks e g mb cb kc d f
m=0.5;
c=120;
k=4e4;
ks=40e4;
e=0.5e-3;
g=9.8;
mb=1;
cb=108;
kc=80e4;
d=1e-3;
f=0.05;
r=sqrt((q(1)-q(5))^2+(q(3)-q(7))^2);
%r=sqrt(q(1)^2+q(3)^2);
if r<d
fx=0;
fy=0;
else
%fx=-(q(1)-q(5))*kc*(1-d/r)+(q(3)-q(7))*f*kc*(1-d/r);
fx=-kc*(1-d/r)*(q(1)-f*q(3));
%fy=-(q(3)-q(7))*kc*(1-d/r)-(q(1)-q(5))*f*kc*(1-d/r);
fy=-kc*(1-d/r)*(f*q(1)+q(3));
end
s=zeros(8,1)
s(1)=q(2);
s(2)=fx/(m*w^2)-c*q(2)/(m*w)-k*q(1)/(m*w^2)-ks*(q(1)^2+q(3)^2)*q(1)/(m*w^2)+e*cos(t);
s(3)=q(4);
s(4)=fy/(m*w^2)-c*q(4)/(m*w)-k*q(3)/(m*w^2)-ks*(q(1)^2+q(3)^2)*q(3)/(m*w^2)+e*sin(t)-g/(w^2);
s(5)=q(6);
s(6)=-fx/(mb*w^2)-cb*q(6)/(mb*w)-kc*q(5)/(mb*w^2);
s(7)=q(8);
s(8)=-fy/(mb*w^2)-cb*q(8)/(mb*w)-kc*q(7)/(mb*w^2)-g/(w^2);
end
function nonlinear
q0=
w=linspace(2200,3001,100)
options=odeset('RelTol',1e-6,'AbsTol',1e-6)
y=[]
for n=1:length(w)
T=2*pi
ts=
=ode45('rk45',ts,q0,options,w(n))
q0=z(end,:)
%figure(1)
subplot(2,1,1)
plot(w(n),z(90000:200:end,1),'r.','markersize',1)
hold on
subplot(2,1,2)
plot(w(n),z(90000:200:end,3),'k.','markersize',1)
hold on
end
%figure(2)
%subplot(2,1,1)
%plot(tt(50000:end),z(50000:end,1),'k-')
%%hold on
y=z(end,:)
save 'end_data.txt'y-ascii
save 'time.txt'tt-ascii
我参考刘长春、吴敬东的《非线性转子-机匣系统的动力学特性》的文章写的一个程序，但是结果就和原文章的对不上，请各位前辈指点迷津，不胜感激 给几章对比下吧 补发文章中的一些参数和图片。 你这分岔图采用的是什么方法？我也在做，跟参考的文献也是出入很大，能大家讨论下吗
权限有限，看不到图，估计是初值有变吧，很多文献上给的程序和文中的图是对不上的
页:
[1]