heteroclinic.asv

计算动力学系统的分岔图

        end
        
        phitemp = BigJacT \ [zeros(homds.ncoords,1);1];
        phinew = phitemp(1:end-1);
        
        psitemp = BigJacT' \ [zeros(homds.ncoords,1);1];
        psinew = psitemp(1:end-1);
        
        OldBigY = [psinew/norm(psinew) phinew/norm(phinew)];
        sigma2 = psinew(end-1);
        
        % Needed small eigenvector
        v1s = VN(:,mu1);
        % Orthogonal complement of stable eigenspace
        Q0S = homds.oldStableQ;
        LSorth = Q0S * [-YS'; eye(size(YS,1))];
        res(13,1) = exp(-D1(mu1) * T) * (v1s' * LSorth * sigma2);
        res(13,1) = sign(res(13,1)) * res(13,1) - cds.options.VarTolerance;
    end
else
    if lambda1 <= homds.nphase
        res(13,1) = -D1(lambda1);
        res(14,1) = D1(lambda1);
    else 
        res(13,1) = -D1(mu1);
        res(14,1) = D1(mu1);
    end
        
    failed = [];
    return
end

if (lambda1 <= homds.nphase) & (~isempty(YU))
    if abs(imag(D1(lambda1))) < cds.options.VarTolerance
        % 14. Inclination-flip with respect to unstable manifold
        if isempty(V)
            [V,D] = eig(A);
            VN = V(:,indlist);
            if homds.nneg == 0
                homds.nneg = homds.nneg + 1;
            elseif homds.nneg == homds.nphase
                homds.nneg = homds.nneg - 1;
            end

            if (isempty(OldBigY)) | (size(OldBigY,2) < 2) | (size(OldBigY,1) ~= homds.ncoords)
                psiold = rand(homds.ncoords,1);
                phiold = rand(homds.ncoords,1);
                BigJacT = [BigJac2 psiold; phiold' 0];
                while cond(BigJacT) > 1e10
                    psiold = rand(homds.ncoords,1);
                    phiold = rand(homds.ncoords,1);
                    BigJacT = [BigJac2 psiold; phiold' 0];
                end            
            else
                psiold = OldBigY(:,1);
                phiold = OldBigY(:,2);   
                BigJacT = [BigJac2 psiold; phiold' 0];         
            end
        
            phitemp = BigJacT \ [zeros(homds.ncoords,1);1];
            phinew = phitemp(1:end-1);
        
            psitemp = BigJacT' \ [zeros(homds.ncoords,1);1];
            psinew = psitemp(1:end-1);
        
            OldBigY = [psinew/norm(psinew) phinew/norm(phinew)];
        end
        
        sigma1 = psinew(end);
        
        % Needed small eigenvector
        v1u = VN(:,lambda1);
        % Orthogonal complement of stable eigenspace
        Q0U = homds.oldUnstableQ;
        LUorth = Q0U * [-YU'; eye(size(YU,1))];
        res(14,1) = exp(D1(lambda1) * T) * (v1u' * LUorth * sigma1);
        res(14,1) = sign(res(14,1)) * res(14,1) - cds.options.VarTolerance;
    end
else
    res(14,1) = D1(mu1);
end

out = res';
failed = [];

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

function [out, failed] = userf(userinf, id, x, v)
global  homds
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(homds.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 homds
[x,x0,p,T,eps0,eps1,YS,YU] = rearr(x);
AP = p(homds.ActiveParams);
failed = 0;
switch(id)    
    case 1
        A = Hom_odejac(x0,num2cell(p));
        D = eig(A);
        nneg = sum(real(D) < 0);

        [D1,indlist] = sort(real(D));
        [val, ind] = min(abs(real(D1)));
    
        if real(D1(ind)) > 0
            lambda1 = ind;
            lambda2 = ind+1;
            lambda3 = ind+2;
            mu1 = ind-1;
            mu2 = ind-2;
            mu3 = ind-3;
        else
            lambda1 = ind+1;
            lambda2 = ind+2;
            lambda3 = ind+3;
            mu1 = ind;
            mu2 = ind-1;
            mu3 = ind-2;
        end
        if (imag(D1(mu1)) == 0) && (imag(D1(lambda1)) == 0)
            fprintf('Neutral saddle, parameters = %g and %g.\n',AP(1),AP(2));
            s.msg  = sprintf('Neutral saddle');
        elseif (imag(D1(mu1)) == 0) || (imag(D1(lambda1)) == 0)
            fprintf('Saddle-focus, parameters = %g and %g.\n',AP(1),AP(2));
            s.msg  = sprintf('Saddle-focus');
        else
            fprintf('Bi-focus, parameters = %g and %g.\n',AP(1),AP(2));
            s.msg  = sprintf('Bi-focus');
        end
    case 2
        fprintf('Double real stable leading eigenvalues, parameters = %g and %g.\n',AP(1),AP(2));
        s.msg  = sprintf('Double real stable leading eigenvalues (saddle - saddle-focus transition)');
    case 3
        fprintf('Double real unstable leading eigenvalues, parameters = %g and %g.\n',AP(1),AP(2));
        s.msg  = sprintf('Double real unstable leading eigenvalues (saddle - saddle-focus transition)');
    case 4
        fprintf('Neutrally-divergent saddle-focus (stable), parameters = %g and %g.\n',AP(1),AP(2));
        s.msg  = sprintf('Neutrally-divergent saddle-focus (stable)');
    case 5
        fprintf('Neutrally-divergent saddle-focus (unstable), parameters = %g and %g.\n',AP(1),AP(2));
        s.msg  = sprintf('Neutrally-divergent saddle-focus (unstable)');
    case 6
        fprintf('Three leading eigenvalues (stable), parameters = %g and %g.\n',AP(1),AP(2));
        s.msg  = sprintf('Three leading eigenvalues (stable)');
    case 7
        fprintf('Three leading eigenvalues (unstable), parameters = %g and %g.\n',AP(1),AP(2));
        s.msg  = sprintf('Three leading eigenvalues (unstable)');
    case 8
        fprintf('Shil''nikov-Hopf (Non-hyperbolic equilibrium), parameters = %g and %g.\n',AP(1),AP(2));
        s.msg  = sprintf('Shil''nikov-Hopf (Non-hyperbolic equilibrium)');
    case 9
        fprintf('Non-central homoclinic-to-saddle-node (Non-hyperbolic equilibrium), parameters = %g and %g.\n',AP(1),AP(2));
        s.msg  = sprintf('Non-central homoclinic-to-saddle-node (Non-hyperbolic equilibrium)');
    case 10
        fprintf('Bogdanov-Takens point, parameters = %g and %g.\n',AP(1),AP(2));
        s.msg  = sprintf('Bogdanov-Takens point');
    case 11
        fprintf('Orbit-flip with respect to the stable manifold, parameters = %g and %g.\n',AP(1),AP(2));
        s.msg  = sprintf('Orbit-flip with respect to the stable manifold');
    case 12
        fprintf('Orbit-flip with respect to the unstable manifold, parameters = %g and %g.\n',AP(1),AP(2));
        s.msg  = sprintf('Orbit-flip with respect to the unstable manifold');
    case 13
        fprintf('Inclination-flip with respect to the stable manifold, parameters = %g and %g.\n',AP(1),AP(2));
        s.msg  = sprintf('Inclination-flip with respect to the stable manifold');
    case 14
        fprintf('Inclination-flip with respect to the unstable manifold, parameters = %g and %g.\n',AP(1),AP(2));
        s.msg  = sprintf('Inclination-flip with respect to the unstable manifold');
end

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

function [S,L] = singmat

  S = 8 * (ones(14,14) - eye(14,14));
  S(1,8) = 1;
%   S(4,[1 9]) = 1;
  S(8,[9 10]) = 1;
  S(9,8) = 0;
  S(9,10) = 1;
  L = ['NS ';'DRS';'DRU';'NDS';'NDU';'3LS';'3LU';'SH ';'NCH';'BT ';'OFS';'OFU';'IFS';'IFU'];

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

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 homds cds

cds.adapted = 1;

YU = homds.YU;
YS = homds.YS;

[x,v] = Hom_adapt_mesh(x,v);

Q0S = homds.oldStableQ;
QbS1 = Q0S(:,1:homds.nneg);
QbS2 = Q0S(:,homds.nneg+1:end);
if ~isempty(YS)
    [U1,S1,R1] = svd(QbS1 + QbS2*YS , 0);
    [U2,S2,R2] = svd(QbS2 - QbS1*YS', 0);
    Q1S = [U1*R1', U2*R2'];
else
    Q1S = Q0S;
end

Q0U = homds.oldUnstableQ;
QbU1 = Q0U(:,1:homds.npos);
QbU2 = Q0U(:,homds.npos+1:end);
if ~isempty(YU)
    [U1,S1,R1] = svd(QbU1 + QbU2*YU , 0);
    [U2,S2,R2] = svd(QbU2 - QbU1*YU', 0);
    Q1U = [U1*R1', U2*R2']; 
else
    Q1U = Q0U;
end
% homds.oldStableQ = Q1S;
% homds.oldUnstableQ = Q1U;

res = 1;

%----------------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ---------------------------------------------------------------
 
function [x,YS,YU,eps1,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;
eps1=x1(idx+1,1);
%eps1=x1(idx+2,1);
x1(end-1:end+1,1),pause

p(hetds.ActiveParams) = x1(end-1:end,1);

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

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

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