Bifurcation diagram plot command?

Question

I am trying to produce a bifurcation diagram.
I am new to the bifurcation world so i need a little help
I have been able to produce the file that runs, however i cant plot the graph.

    for ilop=20000:100:30000
    global  ms Ixx ks1 ks2 coff1 coff2 lf lr pha w Fo pc lp tf
    
        close all
    clear
    %clc   %CLC clears the command window and homes the cursor
    
    load stiffness
    
    %System parameters .
    Fo=0.1;   w=6.9;
    
    ms=730; Ixx=2460; ks1=17500; ks2=17500; 
    coff1=784; coff2=784; 
    lf=1.011; lr=1.803;
    pha=0.3*pi;
    
    ks1=ks1+100;
    
    pc=0;  lp=1;
    xs0=[0  0  0  0  0  0];   stsz=0.001;
    
    
    i0=0;i1=0;i2=0;i3=0;i4=0;i5=0;i6=0;ig=0;   
    
    
    
    n=6000; tf=((2*pi/w)*n);            st=2000;
    
    Istep=2*pi/(w*st);  lep=0:Istep:(tf-Istep); lp1=length(0: (2*pi/(w*st)):(2*pi/w)-(2*pi/(w*st)));  
    
    lop=round(length(lep)/lp1);  
    
    
    xs0=[0.0  0.0  0.0  0.0  0.0  0.0]; 
    
    
    
    %Solver parameters setting
    
    options=odeset('AbsTol',1*exp(-12),'RelTol',1*exp(-12),'InitialStep',0.0000000001,'MaxStep',0.001);
     
    [t x]=ode45(@NLQCM1,[0:(2*stsz):tf],[xs0(1) xs0(2) xs0(3) xs0(4) xs0(5) xs0(6) ]);
    %ode23, ode45, ode113, ode15s, ode23s, ode23t, ode23tb
     
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    
    %round(Nperiod); 
    rtf=round(3.5*length(t)/5);
    xT=[t x];
    Xp=xT(rtf:length(t),: );
    itl=0;
    
    for tL=  1   : lp1   :   length(Xp)-lp1
        itl=itl+1;
        se=[ks1 Xp(tL,:)];    
        xpm(itl,:)=se;
    end
    %Xt=[t x];
    load HCPM
    Xtot=[Xtot; xpm(:,1:6)];
    
    save HCPM Xtot ;
    save Stiffness ks1 ks2;
    end


The function file is

  function [xs F]=HM(t,x)
   
  global i ms Ixx ks1 ks1n ks2 ks2n coff1 coff2  lf lr  w Fo 
   
  %state space vector
  x1=x(1); x2=x(2); x3=x(3); x4=x(4);  x5=x(5); x6=x(6); 
  
  i+1;
  
  xs=[x1 x2  x3  x4  x5 x6]; 
  
  
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %State space model equations
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       u=Fo*sin(w*t); 
       
       Fo=0.05;   w=6.9;
  
  ms=730; Ixx=2460; ks1=17500; ks2=17500; ks1n=34000; ks2n=34000;
  
  coff1=784; coff2=784; 
  
  lf=1.011; lr=1.803;pha=0.3*pi;  
  
      xg1=Fo*sin(w*t); xg2=Fo*sin(w*t-pha); xg1d=Fo*w*cos(w*t); xg2d=Fo*w*cos(w*t-pha);
    
      Xr=x(1)+lr*x(3);
      Xf=x(1)-lf*x(3);
      Xrd=x(2)+x(4)*lr;
      Xfd=x(2)-x(4)*lf;
       
       %if(i>=1)
       %    i=i+1
       %F(i)=u ;
      % end
        
   
       xs(1)=x2;
       
       xs(2)=(-(coff2*Xrd+coff1*Xfd+(ks2*Xr+ks2n*Xr^3)+(ks1*Xf+ks1n*Xf^3))/ms)+((coff2*xg2d+coff1*xg1d+(ks2*xg2+ks2n*xg2^3)+(ks1*xg1+ks1n*xg1^3))/ms);
       
       xs(3)=x4;
       
       xs(4)=(-(coff2*Xrd*lr-coff1*Xfd*lf+(ks2*Xr+ks2n*Xr^3)*lr+(-ks1*Xf-ks1n*Xf^3)*lf)/Ixx)+((coff2*xg2d*lr-coff1*xg1d*lf+(ks2*xg2+ks2n*xg2^3)*lr+(-ks1*xg1-ks1n*xg1^3)*lf)/Ixx);
          
       xs(5)=Fo*cos(w*t); 
       
       xs(6)=Fo*sin(w*t);
      xs=xs';
      
   end

Bifurcation diagram plot command?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

Bifurcation diagram plot command?

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论