|
楼主 |
发表于 2010-4-7 00:09
|
显示全部楼层
回复 板凳 咕噜噜 的帖子
主程序:
function xdot=xiongyuan01(t,x,flag,w);
m1=3010200;
m2=115065;
m3=39120;
m4=39120;
m5=131630;
m6=96600;
k1=29287000000;
k2=31803000000;
k3=126450000000;
k4=149020000000;
k5=31803000000;
k6=55576000000;
c1=1620;
c2=2340;
c3=1000;
c4=900;
c5=1000;
c6=2340;
c7=610;
k7=55576000000;
e=42000000;
p1=1/6;
p2=1/6;
xdot=[x(2);-(-c1+p1*c1*x(2)^2)*x(2)/m1-[-c2+p1*c2*(x(2)-x(4))^2]*(x(2)-x(4))/m1-(k1+p2*k1*x(1)^2)*x(1)/m1-[k2+p2*k2*(x(1)-x(3))^2]*(x(1)-x(3))/m1;
x(4);[-c2+p1*c2*(x(3)-x(1))^2]*(x(2)-x(4))/m2-[-c3+p1*c3*(x(5)-x(3))^2]*(x(4)-x(6))/m2+[k2+p2*k2*(x(3)-x(1))^2]*(x(1)-x(3))/m2-[k3+p2*k3*(x(3)-x(5))^2]*(x(3)-x(5))/m2;
x(6);[-c3+p1*c3*(x(5)-x(3))^2]*(x(4)-x(6))/m3-[-c4+p1*c4*(x(7)-x(5))^2]*(x(6)-x(8))/m3+[k3+p2*k3*(x(5)-x(3))^2]*(x(3)-x(5))/m3-[k4+p2*k4*(x(5)-x(7))^2]*(x(5)-x(7))/m3+e*cos(w*t)/m3;
x(8);[-c4+p1*c4*(x(7)-x(5))^2]*(x(6)-x(8))/m4-[-c5+p1*c5*(x(9)-x(7))^2]*(x(8)-x(10))/m4+[k4+p2*k4*(x(7)-x(5))^2]*(x(5)-x(7))/m4-[k5+p2*k5*(x(7)-x(9))^2]*(x(7)-x(9))/m4-e*cos(w*t)/m4;
x(10);[-c5+p1*c5*(x(9)-x(7))^2]*(x(8)-x(10))/m5-[-c6+p1*c6*(x(11)-x(9))^2]*(x(10)-x(12))/m5+[k5+p2*k5*(x(9)-x(7))^2]*(x(7)-x(9))/m5-[k6+p2*k6*(x(9)-x(11))^2]*(x(9)-x(11))/m5;
x(12);[-c6+p1*c6*(x(11)-x(9))^2]*(x(10)-x(12))/m6-(-c7+p1*c7*x(11)^2)*x(12)/m6+[k6+p2*k6*(x(9)-x(11))^2]*(x(9)-x(11))/m6-(k7+p2*k7*x(11)^2)*x(11)/m6]
分岔图程序:
w=1:4:200;
figure
hold on
for j=1:length(w)
T=2*pi/w(j);
[t,x]=ode23('xiongyuan01',[0:T/200:200*T],[0,0,0,0,0,0,0,0,0,0,0,0],[],w(j));
plot(w(j),x(30000:200:end,5),'k.');
hold on
title ('分岔图') ;xlabel ('w') ;ylabel ('x')
end |
|