变系数系统分岔分析

zizhisuole

主程序是:
clc
clf
Ld=0.6;
r=0.1;
lx=0.5;
p=1;
K1=0.025;
K2=0.03;
m2=10;
Rx=-40;
w=6;  %参数
[T,Y]=ode45('dy2',[0 2*pi/w],[1.001 0 0 0],[],w,lx,r,Ld,p,K1,K2,Rx,m2);
[n1,mm1]=size(Y);
    j(1,:)=[Y(n1,1) Y(n1,2) Y(n1,3) Y(n1,4)];
for i=2:1000
    i
[T,Y]=ode45('dy2',[0 2*pi/w],[Y(n1,1) Y(n1,2) Y(n1,3) Y(n1,4)],[],w,lx,r,Ld,p,K1,K2,Rx,m2);
[n1,mm1]=size(Y);
j(i,:)=[Y(n1,1) Y(n1,2) Y(n1,3) Y(n1,4)];
end;
plot(j(:,1),j(:,2),'.');
xlabel('x1');ylabel('x2');
子程序是：
function dy=dy2(t,y,flag,w,lx,r,Ld,p,K1,K2,Rx,m2)
g=9.81;
k1=K1/(m2*lx*lx);
k2=K2/(m2*lx*lx);
kg=g/lx;
dy=[0,0,0,0];
xg=Ld-r*cos(w*t);
yg=r*sin(w*t);
xgt=r*w*sin(w*t);
ygt=r*w*cos(w*t);
xgtt=r*w^2*cos(w*t);
ygtt=-r*w^2*sin(w*t);
d=sqrt(xg*xg+yg*yg);
k=sqrt(4*lx^2-d^2);
cf1=(xg-yg*k/d)/(2*lx);
cf2=d^2/(2*lx^2)-1;
sf1=(yg+xg*k/d)/(2*lx);
sf2=d*k/(2*lx^2);
cf21=(xg+yg*k/d)/(2*lx);
sf21=(-yg+xg*k/d)/(2*lx);
f1t=-(xgt-ygt*k/d+4*lx^2*yg*(xg*xgt+yg*ygt)/(k*d^3))/(yg+xg*k/d);
f2t=-2*(xg*xgt+yg*ygt)/(d*k);
h1=d*d/(yg^4-xg^4+2*xg*yg*d*k+4*lx*lx*xg^2);
h2=(xg*ygtt-xgt*ygt)*4*lx*lx/d/d-(xg*ygtt-2*xgt*ygt+xgtt*yg)+(xgt^2-ygt^2-xg*xgtt+yg*ygtt)*k/d;
h3=(xg*xgt+yg*ygt)*((xgt*yg-xg*ygt)/d/d-(xg*xgt+yg*ygt)/d/k)-(xgt^2+xg*xgtt+ygt^2+yg*ygtt)*(xg*yg/d/d+yg^2/d/k);
h4=-16*lx*lx*(xg*yg*(d^2-2*lx^2)/d/d-yg^2*(3*lx*lx-d*d)/d/k)*(xg*xgt+yg*ygt)^2/(d^4*k^2);
f1tt=h1*(h2+4*lx*lx*h3/d/d+h4);
f2tt=-2*(2*(xg*xgt+yg*ygt)^2*(d*d-2*lx*lx)/d^2/k^2+xgt^2+xg*xgtt+ygt^2+yg*ygtt)/d/k;
a1=y(2);
a2=(12/(4*p+12-9*cf2^2)*(-1/2*kg*((-p-2)*sf1+sf21)+Rx)+3*(2+3*cf2)/(4*p+12-9*cf2^2)*kg*sf21)*y(1)+(12/(4*p+12-9*cf2^2)*(-k1+f2t*sf2)-6*(2+3*cf2)/(4*p+12-9*cf2^2)*f1t*sf2)*y(2)+(12/(4*p+12-9*cf2^2)*(-1/2*f2t*(f2t-2*f1t)*cf2+1/2*kg*sf21+sf2*f1tt-1/2*sf2*f2tt)+6*(2+3*cf2)/(4*p+12-9*cf2^2)*(-1/2*f1t^2*cf2-1/2*kg*sf21-1/2*sf2*f1tt+1/3*Rx))*y(3)+(12/(4*p+12-9*cf2^2)*(k2-(1/2*f2t-f1t)*sf2-1/2*f2t*sf2)-6*(2+3*cf2)/(4*p+12-9*cf2^2)*k2)*y(4)+12/(4*p+12-9*cf2^2)*(-1/4*kg*((-p-2)*cf1-cf21)*y(1)^2-1/2*kg*cf21*y(1)*y(3)+f2t*cf2*y(2)*y(3)+sf2*y(2)*y(4)+(1/4*f2t*(f2t-2*f1t)*sf2+1/4*kg*cf21+1/2*cf2*f1tt-1/4*cf2*f2tt)*y(3)^2+2*(1/2*(-1/2*f2t+f1t)*cf2-1/4*f2t*cf2)*y(3)*y(4)-1/2*sf2*y(4)^2+sf2*y(3)*((12/(4*p+12-9*cf2^2)*(-1/2*kg*((-p-2)*sf1+sf21)+Rx)+3*(2+3*cf2)/(4*p+12-9*cf2^2)*kg*sf21)*y(1)+(12/(4*p+12-9*cf2^2)*(-k1+f2t*sf2)-6*(2+3*cf2)/(4*p+12-9*cf2^2)*f1t*sf2)*y(2)+(12/(4*p+12-9*cf2^2)*(-1/2*f2t*(f2t-2*f1t)*cf2+1/2*kg*sf21+sf2*f1tt-1/2*sf2*f2tt)+6*(2+3*cf2)/(4*p+12-9*cf2^2)*(-1/2*f1t^2*cf2-1/2*kg*sf21-1/2*sf2*f1tt+1/3*Rx))*y(3)+(12/(4*p+12-9*cf2^2)*(k2-(1/2*f2t-f1t)*sf2-1/2*f2t*sf2)-6*(2+3*cf2)/(4*p+12-9*cf2^2)*k2)*y(4))-1/2*sf2*y(3)*((6*(2+3*cf2)/(4*p+12-9*cf2^2)*(-1/2*kg*((-p-2)*sf1+sf21)+Rx)+6*(p+4+3*cf2)/(4*p+12-9*cf2^2)*kg*sf21)*y(1)+(6*(2+3*cf2)/(4*p+12-9*cf2^2)*(-k1+f2t*sf2)-12*(p+4+3*cf2)/(4*p+12-9*cf2^2)*f1t*sf2)*y(2)+(6*(2+3*cf2)/(4*p+12-9*cf2^2)*(-1/2*f2t*(f2t-2*f1t)*cf2+1/2*kg*sf21+sf2*f1tt-1/2*sf2*f2tt)+12*(p+4+3*cf2)/(4*p+12-9*cf2^2)*(-1/2*f1t^2*cf2-1/2*kg*sf21-1/2*sf2*f1tt+1/3*Rx))*y(3)+(6*(2+3*cf2)/(4*p+12-9*cf2^2)*(k2-(1/2*f2t-f1t)*sf2-1/2*f2t*sf2)-12*(p+4+3*cf2)/(4*p+12-9*cf2^2)*k2)*y(4)))+6*(2+3*cf2)/(4*p+12-9*cf2^2)*(-1/4*kg*cf21*y(1)^2+1/2*kg*cf21*y(1)*y(3)-1/2*sf2*y(2)^2-f1t*cf2*y(2)*y(3)+(1/4*f1t^2*sf2-1/4*kg*cf21-1/4*cf2*f1tt)*y(3)^2-1/2*sf2*y(3)*((12/(4*p+12-9*cf2^2)*(-1/2*kg*((-p-2)*sf1+sf21)+Rx)+3*(2+3*cf2)/(4*p+12-9*cf2^2)*kg*sf21)*y(1)+(12/(4*p+12-9*cf2^2)*(-k1+f2t*sf2)-6*(2+3*cf2)/(4*p+12-9*cf2^2)*f1t*sf2)*y(2)+(12/(4*p+12-9*cf2^2)*(-1/2*f2t*(f2t-2*f1t)*cf2+1/2*kg*sf21+sf2*f1tt-1/2*sf2*f2tt)+6*(2+3*cf2)/(4*p+12-9*cf2^2)*(-1/2*f1t^2*cf2-1/2*kg*sf21-1/2*sf2*f1tt+1/3*Rx))*y(3)+(12/(4*p+12-9*cf2^2)*(k2-(1/2*f2t-f1t)*sf2-1/2*f2t*sf2)-6*(2+3*cf2)/(4*p+12-9*cf2^2)*k2)*y(4)));
a3=y(4);
a4=(6*(2+3*cf2)/(4*p+12-9*cf2^2)*(-1/2*kg*((-p-2)*sf1+sf21)+Rx)+6*(p+4+3*cf2)/(4*p+12-9*cf2^2)*kg*sf21)*y(1)+(6*(2+3*cf2)/(4*p+12-9*cf2^2)*(-k1+f2t*sf2)-12*(p+4+3*cf2)/(4*p+12-9*cf2^2)*f1t*sf2)*y(2)+(6*(2+3*cf2)/(4*p+12-9*cf2^2)*(-1/2*f2t*(f2t-2*f1t)*cf2+1/2*kg*sf21+sf2*f1tt-1/2*sf2*f2tt)+12*(p+4+3*cf2)/(4*p+12-9*cf2^2)*(-1/2*f1t^2*cf2-1/2*kg*sf21-1/2*sf2*f1tt+1/3*Rx))*y(3)+(6*(2+3*cf2)/(4*p+12-9*cf2^2)*(k2-(1/2*f2t-f1t)*sf2-1/2*f2t*sf2)-12*(p+4+3*cf2)/(4*p+12-9*cf2^2)*k2)*y(4)+6*(2+3*cf2)/(4*p+12-9*cf2^2)*(-1/4*kg*((-p-2)*cf1-cf21)*y(1)^2-1/2*kg*cf21*y(1)*y(3)+f2t*cf2*y(2)*y(3)+sf2*y(2)*y(4)+(1/4*f2t*(f2t-2*f1t)*sf2+1/4*kg*cf21+1/2*cf2*f1tt-1/4*cf2*f2tt)*y(3)^2+2*(1/2*(-1/2*f2t+f1t)*cf2-1/4*f2t*cf2)*y(3)*y(4)-1/2*sf2*y(4)^2+sf2*y(3)*((12/(4*p+12-9*cf2^2)*(-1/2*kg*((-p-2)*sf1+sf21)+Rx)+3*(2+3*cf2)/(4*p+12-9*cf2^2)*kg*sf21)*y(1)+(12/(4*p+12-9*cf2^2)*(-k1+f2t*sf2)-6*(2+3*cf2)/(4*p+12-9*cf2^2)*f1t*sf2)*y(2)+(12/(4*p+12-9*cf2^2)*(-1/2*f2t*(f2t-2*f1t)*cf2+1/2*kg*sf21+sf2*f1tt-1/2*sf2*f2tt)+6*(2+3*cf2)/(4*p+12-9*cf2^2)*(-1/2*f1t^2*cf2-1/2*kg*sf21-1/2*sf2*f1tt+1/3*Rx))*y(3)+(12/(4*p+12-9*cf2^2)*(k2-(1/2*f2t-f1t)*sf2-1/2*f2t*sf2)-6*(2+3*cf2)/(4*p+12-9*cf2^2)*k2)*y(4))-1/2*sf2*y(3)*((6*(2+3*cf2)/(4*p+12-9*cf2^2)*(-1/2*kg*((-p-2)*sf1+sf21)+Rx)+6*(p+4+3*cf2)/(4*p+12-9*cf2^2)*kg*sf21)*y(1)+(6*(2+3*cf2)/(4*p+12-9*cf2^2)*(-k1+f2t*sf2)-12*(p+4+3*cf2)/(4*p+12-9*cf2^2)*f1t*sf2)*y(2)+(6*(2+3*cf2)/(4*p+12-9*cf2^2)*(-1/2*f2t*(f2t-2*f1t)*cf2+1/2*kg*sf21+sf2*f1tt-1/2*sf2*f2tt)+12*(p+4+3*cf2)/(4*p+12-9*cf2^2)*(-1/2*f1t^2*cf2-1/2*kg*sf21-1/2*sf2*f1tt+1/3*Rx))*y(3)+(6*(2+3*cf2)/(4*p+12-9*cf2^2)*(k2-(1/2*f2t-f1t)*sf2-1/2*f2t*sf2)-12*(p+4+3*cf2)/(4*p+12-9*cf2^2)*k2)*y(4)))+12*(p+4+3*cf2)/(4*p+12-9*cf2^2)*(-1/4*kg*cf21*y(1)^2+1/2*kg*cf21*y(1)*y(3)-1/2*sf2*y(2)^2-f1t*cf2*y(2)*y(3)+(1/4*f1t^2*sf2-1/4*kg*cf21-1/4*cf2*f1tt)*y(3)^2-1/2*sf2*y(3)*((12/(4*p+12-9*cf2^2)*(-1/2*kg*((-p-2)*sf1+sf21)+Rx)+3*(2+3*cf2)/(4*p+12-9*cf2^2)*kg*sf21)*y(1)+(12/(4*p+12-9*cf2^2)*(-k1+f2t*sf2)-6*(2+3*cf2)/(4*p+12-9*cf2^2)*f1t*sf2)*y(2)+(12/(4*p+12-9*cf2^2)*(-1/2*f2t*(f2t-2*f1t)*cf2+1/2*kg*sf21+sf2*f1tt-1/2*sf2*f2tt)+6*(2+3*cf2)/(4*p+12-9*cf2^2)*(-1/2*f1t^2*cf2-1/2*kg*sf21-1/2*sf2*f1tt+1/3*Rx))*y(3)+(12/(4*p+12-9*cf2^2)*(k2-(1/2*f2t-f1t)*sf2-1/2*f2t*sf2)-6*(2+3*cf2)/(4*p+12-9*cf2^2)*k2)*y(4)));
dy=[a1;a2;a3;a4];
但是运行之后结果一直不对，希望大家帮忙看一下~谢谢！