heteroclinic.asv

计算动力学系统的分岔图

function out = heteroclinic
%
% heteroclinic curve definition file for a problem in mapfile
% 

global hetmds cds
    out{1}  = @curve_func;
    out{2}  = @defaultprocessor;
    out{3}  = @options;
    out{4}  = @jacobian;
    out{5}  = @hessians;
    out{6}  = @testf;
    out{7}  = [];%@userf;
    out{8}  = @process;
    out{9}  = @singmat;
    out{10} = [];%@locate;
    out{11} = @init;
    out{12} = @done;
    out{13} = @adapt;
return


%----------------------------------------------------
function func = curve_func(arg)
  
  [x,YS,YU,p] = rearr(arg);
  func = BVP_Het(x,YS,YU,p);

 
%-----------------------------------------------------
function varargout = jacobian(varargin)
  global hetds  cds  
  J= hetds.niteration;
  [x,YS,YU,p] = rearr(varargin{1});
  varargout{1} = BVP_Het_jac(hetds.func,x,p,YS,YU,J);
%  varargout{1},'varargout{1}',pause
  %BigJac = varargout{1};
%-----------------------------------------------------

function varargout = hessians(varargin)

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

function varargout = defaultprocessor(varargin)
global hetds opt cds
   
  [x,YS,YU,p] = rearr(varargin{1});
  v = rearr(varargin{2});  
 
  p1 = num2cell(p);
  hetds.YS = YS;
  hetds.YU = YU;
 
  if nargin > 2
    % set data in special point structure
    s = varargin{3};
    varargout{3} = s;
  end
%   % all done succesfully
   varargout{1} = 0;
   varargout{2} = hetds.npoints';

%-------------------------------------------------------
  
function option = options
global hetds cds
  % Check for symbolic derivatives in mapfile
% hetds,pause, hetds.Jacobian,pause
     symjac  = ~isempty(hetds.Jacobian);
     symhes  = ~isempty(hetds.Hessians);
     symder  = ~isempty(hetds.Der3);
% %   
     symord = 0; 
      if symjac, symord = 1; end
      if symhes, symord = 2; end
    
%    if symder, symord = 3; end
  %if higher>2, symord = higher; end

  option = contset;
 option = contset(option, 'SymDerivative', symord);
 %option = contset(option, 'Workspace', 1);
%  option = contset(option, 'Locators', zeros(1,13));
  symjacp = ~isempty(hetds.JacobianP); 
%   symhes  = ~isempty(hetds.HessiansP);
%   symordp = 0;
   if symjacp, symordp = 1; end
   if symhes,  symordp = 2; end
   option = contset(option, 'SymDerivativeP', symordp);
  
  cds.symjac  = 1;
  cds.symhess = 1;
  
%------------------------------------------------------  
  
function [out, failed] = testf(id, x0, v)
global hetds cds        

[x,YS,YU,p] = rearr(x0);
p = n2c(p);
ndim = cds.ndim;
J=contjac(x0);%eig((J(:,1:ndim-1))+eye(ndim-1)),
failed = [];
for i=id
  lastwarn('');
  
  switch i
     
case 1 % LP
     out(1) = v(end);
 case 2 % BP
      B = [J; v'];
      out(2) = det(B);
   
  otherwise
    error('No such testfunction');
  end
  if ~isempty(lastwarn)
    msg = sprintf('Could not evaluate tf %d\n', i);
    failed = [failed i];
  end
  
end

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

function [out, failed] = userf(userinf, id, x, v)
global  hetds cds
dim =size(id,2);
failed = [];
[x,x0,p,T,eps0,eps1,YS,YU] = rearr(x); p = num2cell(p);
for i=1:dim
  lastwarn('');
  if (userinf(i).state==1)
      out(i)=feval(hetds.user{id(i)},0,x0,p{:});
  else
      out(i)=0;
  end
  if ~isempty(lastwarn)
    msg = sprintf('Could not evaluate userfunction %s\n', id(i).name);
    failed = [failed i];
  end
end

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

function [failed,s] = process(id, x, v, s)
global  cds hetds
[x0,YS,YU,p] = rearr(x);
p = n2c(p);
ndim = cds.ndim; 
nphase=hetds.nphase;
n=hetds.niteration;
 % WM: Removed SL array
fprintf('label = %s, x = ', s.label); printv(x)
p1=p;
 switch id     
  case 1 % LP      
     jac =hetjac(x,p,n);
     [V,D]= eig(jac-eye(nphase));
     [Y,i]=min(abs(diag(D)));
     vext=real(V(:,i));
     vext=vext/norm(vext);
     [V,D]= eig(jac'-eye(nphase));
     [Y,i]=min(abs(diag(D)));
     wext=real(V(:,i));      
     wext=wext/(vext'*wext);
     s.msg=sprintf('Limit point\n');
       
  case 2 %BP
      s.msg=sprintf('Branch point\n');  
      s.data.v=v;
end

% Compute eigenvalues for every singularity
J=contjac(x);
if ~issparse(J)
  [v,d]=eig(J(:,1:ndim-1));
else
  opt.disp=0;
  % WM: fixed a bug (incompatability between MatLab 6.0 and 5.5?)
  [v,d]=eigs(J(:,1:ndim-1),min(6,ndim-1),'lm',opt);
end
%d=d+eye(nphase);
s.data.evec = v;
%s.data.eval = diag(d)';

failed = 0;

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

function [S,L] = singmat
global hetds cds
 
% 0: testfunction must vanish
% 1: testfunction must not vanish
% everything else: ignore this testfunction

  S = [  0 8 
         8 1 ] ;

  L = [  'LP  ';'BP  ' ];


  %elseif strcmp(arg, 'locate')

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

function [x,v] = locate(id, x1, v1, x2, v2)
msg = sprintf('No locator defined for singularity %d', id);
error(msg);
    
%----------------------------------------------------------

function varargout = init(varargin)

WorkspaceInit(varargin{1:2});
% all done succesfully
varargout{1} = 0;

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

function varargout = done

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

function [res,x,v] = adapt(x,v)
global hetds cds

res = []; % no re-evaluations needed


cds.adapted = 1;
% 
YU = hetds.YU;
YS = hetds.YS;
QU=hetds.Q0;
QS = hetds.Q1;
[U,S,V]=svd(QU * [ eye(size(YU,1));YU]);
hetds.Q0= U;
hetds.YU=zeros(hetds.nphase-hetds.nu,hetds.nu);

[U,S,V]=svd(QS * [ eye(size(YS,1));YS]);
hetds.Q1= U;
hetds.YS=zeros(hetds.nphase-hetds.ns,hetds.ns);
res = 1;
%hetds.Q0,hetds.Q1,pause
%----------------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ---------------------------------------------------------------
 
function [x,YS,YU,p] = rearr(x1)
% Rearranges x1 into all of its components
global hetds

x = x1(1:hetds.nphase*hetds.npoints,1);
p = hetds.P0;
% eps0 = hetds.eps0;
% eps1 = hetds.eps1;
idx=hetds.npoints*hetds.nphase;%+hetds.nu.hetds.ns;
ju=hetds.nphase-hetds.nu;
js=hetds.nphase-hetds.ns;
YU = reshape(x1(idx+1:idx+ju*hetds.nu,1),hetds.nphase-hetds.nu,hetds.nu);
idx = idx + ju*hetds.nu;
YS = reshape(x1(idx+1:idx+js*hetds.ns,1),hetds.nphase-hetds.ns,hetds.ns);
idx = idx + js*hetds.ns;
p(hetds.ActiveParams) = x1(end,1);

%size(x),pause,YS,YU,eps0,eps1,p,pause
    
% -------------------------------------------------------------

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

function WorkspaceInit(x,v)
global cds hetds
% hetds.cols_p1 = 1:(hetds.ncol+1);
% hetds.cols_p1_coords = 1:(hetds.ncol+1)*hetds.nphase;
% hetds.ncol_coord = hetds.ncol*hetds.nphase;
% hetds.col_coords = 1:hetds.ncol*hetds.nphase;
% hetds.coords = 1:hetds.ncoords;
% hetds.pars = hetds.ncoords+(1:2);
% hetds.tsts = 1:hetds.ntst;
% hetds.cols = 1:hetds.ncol;
% hetds.phases = 1:hetds.nphase;
% hetds.ntstcol = hetds.ntst*hetds.ncol;
% 
% hetds.idxmat = reshape(fix((1:((hetds.ncol+1)*hetds.ntst))/(1+1/hetds.ncol))+1,hetds.ncol+1,hetds.ntst);
% hetds.dt = hetds.msh(hetds.tsts+1)-hetds.msh(hetds.tsts);
% 
% hetds.wp = kron(hetds.wpvec',eye(hetds.nphase));
% hetds.pwwt = kron(hetds.wt',eye(hetds.nphase));
% hetds.pwi = hetds.wi(ones(1,hetds.nphase),:);
% 
% hetds.wi = nc_weight(hetds.ncol)';
% 
% [hetds.bialt_M1,hetds.bialt_M2,hetds.bialt_M3,hetds.bialt_M4]=bialtaa(hetds.nphase);

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

function [x,v,s] = WorkspaceDone(x,v,s)

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

function K = fastkron(c,p,A,B)
t = p:((c+2)*p-1);
K = A(ones(1,p),fix(t/p)).*B(:,rem(t,p)+1);