umanifold.m

计算动力学系统的分岔图

            
function M = Umanifold(p0,Arc,alphamax, alphamin,deltamax,deltamin,deltaalphamax,deltaalphamin,deltak,pp,epsb,lamb,del,v,n)
 %Compute 1-D unstable manifold           
global M              
%pp is vector of parameters
p1=p0+del*v;
M=[p0,p1];p = n2c(pp);
preimage=p0-del/lamb;
pleft=p0;lindex=1;
pright=p1;rindex=2;
arcl=norm(p1-p0);
while( arcl < Arc)
    [pcan,preimage,pi,pj,lindex,rindex]=Search_Line(M, pleft,pright,preimage,lindex,rindex, deltak,deltamin,alphamax,pp,epsb,deltamax,n);
    
    pbar=M(:,end)+((norm(M(:,end)-M(:,end-1))/norm(M(:,end)-pcan)))*(M(:,end)-pcan);
    alphak=norm(pbar-M(:,end-1))/norm(M(:,end)-M(end-1));
    
    deltak=norm(pcan-M(:,end)); 
    
  
  if ((alphak<alphamax)&& (deltak*alphak<deltaalphamax))   
        
        M(:,end+1)=pcan;
        arcl=arcl+deltak;
        pleft=preimage;
        pright=pj;        
   
   elseif ((alphak<alphamax)&& (deltak*alphak<deltamax))
    
      deltak=2*deltak;
  elseif deltak > deltamin
      deltak=deltak/2;
           
  end %if    
 
end %while

%-----------------------------------------------------------------------
function [pcan,preimage,pi,pj,lindex, rindex]= Search_Line(M, pleft,pright,preimage,lindex,rindex, deltak,deltamin,alphamax,pp,epsb,deltamax,n);

[S,pleft,pright,lindex,rindex,deltak]=Find_intersection_On_circle(M, pleft,pright,preimage,lindex,rindex,deltak,deltamin,alphamax,pp,epsb,deltamax,n);
[pcan,preimage,pleft,pright] =Bisection(M,pleft,pright,preimage,deltak,epsb,pp,alphamax,n) ;
 pi=pleft;pj=pright;
%--------------------------------------------------------------------------

function [S,pleft,pright,lindex,rindex,deltak]=Find_intersection_On_circle(M, pleft,pright,preimage,lindex,rindex,deltak,deltamin,alphamax,pp,epsb,deltamax,n)
                  
  [S,pleft,pright,lindex,rindex]=Side_of_line(M, pleft,pright,preimage,lindex,rindex,deltak,deltamin,alphamax,pp,epsb,deltamax,n);
  j=0;
  while (S==0)&&(j<100)   
      j=j+1;
      if M(:,rindex)~=M(:,end)        
         lindex=rindex;rindex=rindex+1;
         pleft=M(:,lindex); pright=M(:,rindex); 
      elseif deltak>deltamin
          deltak=deltak/2;
      end
      [S,pleft,pright,lindex,rindex]=Side_of_line(M, pleft,pright,preimage,lindex,rindex,deltak,deltamin,alphamax,pp,epsb,deltamax,n);
  end   

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

function [S,pleft,pright,lindex,rindex] =Side_of_line(M, pleft,pright,preimage,lindex,rindex,deltak,deltamin,alphamax,pp,epsb,deltamax,n)

p = n2c(pp);
fp=norm(fun_eval(0,preimage,p{:},n)-M(:,end));  
fr=norm(fun_eval(0,pright,p{:},n)-M(:,end));
if ((fp<deltak)&&(fr>deltak))|((fp>deltak)&&(fr<deltak))
%The line segment L that is mapped by f,  intersects the circle with center pk and radius deltak 
    S=1;
else
    S=0;
end

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

function [pcan,preimage,pleft,pright] =Bisection(M,pleft,pright,preimage,deltak,epsb,pp,alphamax,n)  
 pleft=preimage;
 ptry=(pleft+pright)/2;
 p=n2c(pp);
 ft=norm(fun_eval(0,ptry,p{:},n)-M(:,end));
 fl=norm(fun_eval(0,preimage,p{:},n)-M(:,end));
 fr=norm(fun_eval(0,pright,p{:},n)-M(:,end));
  if (abs(fl)> abs(fr))
      temp=fl;
      fl=fr;
      fr=temp;
  end     

i=0;
 while (ft>(1+epsb)*deltak) && (ft<(1-epsb)*deltak)&&(i<10) 
     i=i+1;
  if (ft>= deltak)
     % f(ptry) is on the same side as f(pend)     
      pright=ptry;
      ptry=(pleft+pright)/2; 
     
  else   
     pleft=ptry;
     ptry=(pleft+pright)/2;
  end  
  
 end %while 

 p=n2c(pp);ftry=fun_eval(0,ptry,p{:},n);
 pcan=ftry;preimage=ptry;
  %------------------------------------------------------------------------  

  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

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