如何用以下程序作出poincare图或混沌分叉图

shipo507

第一部分：

function df=rhs(t,f,data)
%方程中的系统参数
m=data(1);g=data(2);c=data(3);n=data(4);b=data(5);mf=data(6);M=data(7);D0=data(8);De=data(9);
tau0=data(10);Ke=data(11);K0=data(12);omega=data(13);e=data(14);
%化简后的动力学方程
df(1)=f(2);
df(2)=-(De+D0*(1-f(1)^2-f(3)^2)^(-n))*f(2)/(M*omega)-2*tau0*mf*f(4)*(1-sqrt(f(1)^2+f(3)^2))^b/M-(Ke+K0*(1-f(1)^2-f(3)^2)^(-n)-mf*omega^2*tau0^2*(1-sqrt(f(1)^2+f(3)^2))^(2*b))*f(1)/(M*omega^2)-(tau0*D0*(1-sqrt(f(1)^2+f(3)^2))^b*(1-f(1)^2-f(3)^2)^(-n))*f(3)/(M*omega)+m*e*cos(t)/(M*c);
df(3)=f(4);
df(4)=-(De+D0*(1-f(1)^2-f(3)^2)^(-n))*f(4)/(M*omega)+2*tau0*mf*f(2)*(1-sqrt(f(1)^2+f(3)^2))^b/M-(Ke+K0*(1-f(1)^2-f(3)^2)^(-n)-mf*omega^2*tau0^2*(1-sqrt(f(1)^2+f(3)^2))^(2*b))*f(3)/(M*omega^2)+(tau0*D0*(1-sqrt(f(1)^2+f(3)^2))^b*(1-f(1)^2-f(3)^2)^(-n))*f(1)/(M*omega)+m*e*sin(t)/(M*c)-m*g/(M*c*omega^2);
df=[df(1);df(2);df(3);df(4)];


第二部分：

clc;
clear;
format long
m=50;       Ke=7.2762e6;    De=2e3;
n0=0.079;    m0=-0.25;    g=9.8;
e=0.0002;    R=0.067;      l=0.102;
n=2.5;       b=0.45;     tau0=0.4;   c=0.0003;      nu=50;
z=0.1;       p=2e5;        delta=0.0003;    T=l/nu;     upsilon=1.47e-5;
fid=fopen('poincare.txt','w');
fidd=fopen('time.txt','w');
fiddd=fopen('center.txt','w');
% 转速
omega=1300;
% 初始位置
x0=0;
y0=-m*g/Ke;
X0=x0/c;Y0=y0/c;
% 无量纲转速
nq=length(omega);
for i=1:nq
       s(i)=omega(i)/sqrt(Ke/m);
% Child短密封轴承动力系数
Rv(i)=omega(i)*R*delta/upsilon;
Ra=2*nu*delta/upsilon;
lambda(i)=n0*Ra^m0*(1+(Rv(i)/Ra)^2)^((1+m0)/2);
sigma(i)=lambda(i)*l/delta;
B(i)=2-((Rv(i)/Ra)^2-m0)/((Rv(i)/Ra)^2+1);
E(i)=(1+z)/2/(1+z+2*sigma(i));
u0(i)=2*sigma(i)^2*E(i)*(1-m0)/(1+z+2*sigma(i));
u1(i)=2*sigma(i)^2*(E(i)/sigma(i)+B(i)*(1/6+E(i))/2)/(1+z+2*sigma(i));
u2(i)=sigma(i)*(1/6+E(i))/(1+z+2*sigma(i));
u3(i)=pi*R*p/lambda(i);
% 密封力Muszynska模型系数
% K=K0*(1-e^2)^(-n)
% D=D0*(1-e^2)^(-n)
% tau=tau0*(1-e)^b
% e=sqrt(x^2+y^2)/c
K0(i)=u0(i)*u3(i);  
D0(i)=u1(i)*u3(i)*T;
mf(i)=u2(i)*u3(i)*(T^2);
M(i)=mf(i)+m;
% 系统结构动力学方程中的参数
data=[m,g,c,n,b,mf(i),M(i),D0(i),De,tau0,Ke,K0(i),omega(i),e];
%初始条件
f0=[X0;0;Y0;0];
t0=0;h=0.001;tf=600;
tn=[t0:h:tf];
[t,f]=ode45(@rhs,tn,f0,[],data);
%求出的是无量纲化的位移、速度
%现转化为实际的位移、速度
x=f(:,1);y=f(:,3);
dx=f(:,2);dy=f(:,4);
nn=(tf-t0)/h+1;nnn=300000;%nnn为除去瞬态响应后的第一点数，nn为采集的总点数
nnnn=400000;
for i=nnn:nn
fprintf(fidd,'%6.4f  %12.8f\n',t(i),x(i));hold on;
end
fclose(fidd);
for j=nnn:nn
    fprintf(fiddd,'%12.8f %12.8f\n',x(j),y(j));hold on;
end
fclose(fiddd);

% % 求庞加莱映射点
[A,jj]=max(x(nnn:nnnn));
jj=jj+nnn-1;
TT=2*pi;
  for k=1:nn
    tt(k)=(jj-1)*h+k*TT;
      if tt(k)>tf
          break;
      end
      for kk=nnn:nn
          if tt(k)>(kk-1)*h&tt(k)<kk*h
              l1=(tt(k)-(kk+1)*h)*(tt(k)-kk*h)/(2*h^2);
              l2=(tt(k)-(kk-1)*h)*(tt(k)-(kk+1)*h)/(-h^2);
              l3=(tt(k)-(kk-1)*h)*(tt(k)-kk*h)/(2*h^2);
          xx(k)=l1*x(kk-1)+l2*x(kk)+l3*x(kk+1);
          dx(k)=l1*dx(kk-1)+l2*dx(kk)+l3*dx(kk+1);
          fprintf(fid,'%12.8f %12.8f\n',xx(k),dx(k));hold on;
          end
      end
  end  
end
fclose(fid);