|
|
马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。
您需要 登录 才可以下载或查看，没有账号？我要加入
x
大家好，我有一个做稳定流形和不稳定流形图的程序，可是老是运行不出来，请高手指点。
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
%--------------------------------------------------------------------------
|
|
发表于 2012-12-16 14:30
提升卡
变色卡
显身卡
发表于 2012-12-16 15:35
楼主
发表于 2012-12-17 10:33
