方程如下:
function dy=accele(t,y,r1)
dy=zeros(2,1);
f=0.8; a=0.001; K=2000; L=0.1;
r0=3; w=3; r=r0+r1*sin(w*t);
k11=(r^2-1-(f^2+a*r)*(pi)^2)*(pi)^2;
k12=3*K^2*(pi)^4/8;
k13=3*L*(pi)^4/4;
dy(1)=y(2);
dy(2)=k11*y(1)-k12*(y(1))^3-k13*(y(1))^2*y(2);
分叉图程序:clear
A=0:0.01:3;
options=odeset('RelTol',1e-7);
hold on
for j=1:length(A)
r1=A(j);
tt=2*pi/3;
[t,y]=ode45(@accele,[0:tt/100:80*tt],[0,0],options,r1);
plot(A(j),y(8000:100:end,1),'k')
end
不知道哪里出了问题,还请多多帮助,谢谢。
|