smanifold.m

计算动力学系统的分岔图

function SM = Smanifold(p0,p1,Arc,alphamax, alphamin,deltamax,deltamin,dalphamax,dalphamin, deltak,pp,epsb,n)
%Compute 1-D stable manifold of a map without inverse           
         
%pp is vector of parameters

global SM 
SM=[p0,p1];
lindex=1;
rindex=2;
arcl=norm(p1-p0);
while( arcl < Arc)
    [pcan,indexi,indexj]=Search_Cicle(SM,SM(:,lindex),lindex,SM(:,rindex),rindex,deltak, dalphamax,pp,epsb,n);    
    pbar=SM(:,end)+((norm(SM(:,end)-SM(:,end-1))/norm(SM(:,end)-pcan)))*(SM(:,end)-pcan);
    alphak=2*asin(norm(pbar-SM(:,end-1))/(2*norm(SM(:,end)-SM(:,end-1))));
    deltak=norm(pcan-SM(:,end));
    
        
    if ((alphak<alphamax)&& (deltak*alphak<dalphamax))
       SM(:,end+1)=pcan;
       arcl=arcl+deltak;
       lindex=indexi;
       rindex=indexj;
       
    elseif ((alphak< alphamin)&&((alphak*deltak)< dalphamin))  
           deltak=2*deltak;
    else   
             deltak=deltak/2;
     end 
  
end 

%--------------------------------------------------------------------------

function [pcan,indexi, indexj]= Search_Cicle(SM, pleft,lindex,pright,rindex,deltak,alphamax,pp,epsb,n)
                                              
[pcan,tau]=Find_Point_On_Line(SM,pleft,pright,deltak,alphamax,pp,epsb,n);
R=0;
while ((tau<0)|(tau>1))&& (R<1000)   
  if (tau <0) %move backward
     rindex=lindex;
     lindex=lindex-1;
     pleft=SM(:,lindex);
     pright=SM(:,rindex);
     
  else  %move forward
        lindex=rindex;
        rindex=rindex+1;
        pright=SM(:,rindex);
        pleft=SM(:,lindex);
               
  end 
  R=R+1;
  [pcan,tau]=Find_Point_On_Line(SM,pleft,pright,deltak,alphamax,pp,epsb,n);    
end  
if R > 999  
   debug('detection of a sharp point!');
   stop
   
end

 indexi=lindex; indexj=rindex;
%--------------------------------------------------------------------------

function [pcan, tau]= Find_Point_On_Line(SM,pleft,pright,deltak,alphamax, pp,epsb,n)
                
  S=Side_of_line(SM,pleft,pright,deltak,alphamax,pp,n);
  
  while (S==0)
      if alphamax<0.5 
          alphamax=2*alphamax;
      else
          deltak=2*deltak;
      end
      
      S=Side_of_line(SM,pleft,pright,deltak,alphamax,pp,n);   
  end   
  p = n2c(pp);
 [ptry,pi,pj,tau] =Bisection(SM,pleft,pright,deltak,alphamax,epsb,pp,alphamax,n);
 pcan=ptry;
 
 
% -------------------------------------------------------------------------

function S =Side_of_line(SM,pleft,pright,deltak,alphamax,pp,n) 
p = n2c(pp);                                    
L=-pright+pleft;
tstart=-alphamax;tend=alphamax;
VV=(deltak*(SM(:,end)-SM(:,end-1)))/norm(SM(:,end)-SM(:,end-1));
VV1=[-VV(2,:); VV(1,:)]; 
pstart=SM(:,end)+VV*cos(tstart)+VV1*sin(tstart);
pend=SM(:,end)+VV*cos(tend)+VV1*sin(tend);
Vstart=fun_eval(0,pstart,p{:},n)-pright;
Vend=fun_eval(0,pend,p{:},n)-pright;

W=[-L(2,:); L(1,:)];

if (W'*Vstart)*(W'*Vend)>=0
    % f(pstart) and f(pend) do not lie on opposite sides of L
    S=0;
else
    S=1;
end
 
    %  --------------------------------------------------------------------

function [ptry,  pleft, pright,tau] =Bisection(SM,pleft,pright,deltak,alphamax,epsb,pp,amax,n)  
 p=n2c(pp);  
 L=-pright+pleft;
 W=[-L(2,:); L(1,:)];
 tstart=-alphamax;tend=alphamax;
 ttry=(tstart+tend)/2.0; 
 VV=deltak*(SM(:,end)-SM(:,end-1))/norm(SM(:,end)-SM(:,end-1));
 VV1=[-VV(2,:); VV(1,:)];
 ptry=SM(:,end)+VV*cos(ttry)+VV1*sin(ttry);
 pstart=SM(:,end)+VV*cos(tstart)+VV1*sin(tstart);
 pend=SM(:,end)+VV*cos(tend)+VV1*sin(tend);
 V=fun_eval(0,ptry,p{:},n)-pright;
 Vstart=fun_eval(0,pstart,p{:},n)-pright;
 Vend=fun_eval(0,pend,p{:},n)-pright;
 
i=0;
 while (abs( V'*W)> epsb)&&(i<150) 
     i=i+1;
  if ((Vend'*W)*(V'*W)>=0)
     % f(ptry) is on the same side as f(pend)
      tend=ttry;
      pend=SM(:,end)+VV*cos(tend)+VV1*sin(tend);
      ttry=(tstart+tend)/2.0;
      ptry=SM(:,end)+VV*cos(ttry)+VV1*sin(ttry); 
      
  else
     % f(ptry) is on the same side as f(pstart)
     tstart=ttry;
     pstart=SM(:,end)+VV*cos(tstart)+VV1*sin(tstart);
     ttry=(tstart+tend)/2.0;     
     ptry=SM(:,end)+VV*cos(ttry)+VV1*sin(ttry);
          
  end  
  Vend=fun_eval(0,pend,p{:},n)-pright;
  V=fun_eval(0,ptry,p{:},n)-pright;  
 end %while 
 
fptry=fun_eval(0,ptry,p{:},n);
tau=-((fptry-pleft)'*(pleft-pright))/(norm(pleft-pright)^2);

  %------------------------------------------------------------------------   

  function map = fun_eval(t,kmrgd,alpha,beta,R,S,n)
     

   for i=1:n       
    map=[kmrgd(2); alpha-beta*kmrgd(1)-kmrgd(2)^2+R*kmrgd(1)*kmrgd(2)+S*kmrgd(2);];
    kmrgd(1)=map(1);
    kmrgd(2)=map(2);
   end
%--------------------------------------------------------------------------