稳定流形和不稳定流形图

yanxiangya

大家好，我有一个做稳定流形和不稳定流形图的程序，可是老是运行不出来，请高手指点。
            
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
%--------------------------------------------------------------------------  

ChaChing

没报错吗!?
若有的话, 建议给齐报错讯息

yanxiangya

ChaChing 发表于 2012-12-16 15:35

                               
登录/注册后可看大图

没报错吗!?
若有的话, 建议给齐报错讯息

可是根据报错信息修改出来还是不正确呀，现在都不会修改了

happy

yanxiangya 发表于 2012-12-16 16:55

                               
登录/注册后可看大图

可是根据报错信息修改出来还是不正确呀，现在都不会修改了

请描述具体现象