init_hett_hett.asv

计算动力学系统的分岔图

function [x0,v0]= init_HetT_HetT(mapfile, X,nphase,nu,ns, p,ap,J)

% Initializes a fold continuation from a LP point
%
global cds hetTds
% check input
if size(ap,2)~=2
    errordlg('Two active parameter are needed for a heteroclinic tangency curve continuation');
end
% initialize hetTds
 %X=[x1,..., xn2]^T, xi,i=1,...,n2 are  mesh points 
n=size(X,1);
hetTds.nphase=nphase;
k=nu*(hetTds.nphase-nu)+ns*(hetTds.nphase-ns);
hetTds.npoints=(n-k)/hetTds.nphase;
hetTds.Niterations=J;
hetTds.sizep = size(p,1);
hetTds.mapfile = mapfile;
func_handles = feval(hetTds.mapfile);
hetTds.func = func_handles{2};
hetTds.Jacobian  = func_handles{3};
hetTds.JacobianP = func_handles{4};
hetTds.Hessians  = func_handles{5};
hetTds.HessiansP = func_handles{6};
hetTds.Der3      = func_handles{7};
hetTds.Niterations=J;
siz = size(func_handles,2);
if siz > 9
    j=1;
    for i=10:siz
        hetTds.user{j}= func_handles{i};
        j=j+1;
    end
end

hetTds.ActiveParams = ap;
hetTds.P0 = p;
hetTds.p0(ap) = p(hetTds.ActiveParams);
p=n2c(p);
%if size(varargin,1)>0,hetTds.BranchParams=varargin{1};else hetTds.BranchParams=[];end 
cds.curve = @heteroclinicT;
hetTds.nu=nu;
hetTds.ns=ns;
[x,YS,YU,p] = rearr(X);
hetTds.YS=YS;
hetTds.YU=YU;
x0=[x;YU;YS;(hetTds.p0(ap))'];
x1=x(1:2,1);xN=x(end-1:end,1);
%x0(hetTds.nphase+1:hetTds.nphase+2,:) = hetTds.P0(ap);
A1= hetT_jac(x1,p,J);
AN = hetT_jac(xN,p,J);
[QU, eigvlU, dimU] = Het_computeBase(A1,1,hetTds.nu);
[QS,eigvlS, dimS] = Het_computeBase(AN,0,hetTds.ns);
hetTds.Q0 = QU;
hetTds.Q1 = QS;
cds.ndim = length(x0)+2;
curvehandles = feval(cds.curve);
cds.curve_func = curvehandles{1};
cds.curve_options = curvehandles{3};
cds.curve_jacobian = curvehandles{4};
cds.curve_hessians = curvehandles{5}; 
cds.options = feval(cds.curve_options);
cds.options = contset(cds.options,'Increment',1e-5);
jac =BVP_HetT_jac(x,p,YS,YU,J);
% calculate eigenvalues
[U,S,V]=svd(full(jac));
hetTds.b=U(:,end);
hetTds.c=V(:,end);
v0=[];
rmfield(cds,'options'); 
size(x0),pause,x0),pause


% ---------------------------------------------------------------
function [x,YS,YU,p] = rearr(x1)
% Rearranges x1 into all of its components
global hetTds
ap=hetTds.ActiveParams;
x = x1(1:hetTds.nphase*hetTds.npoints,1);
p = hetTds.P0;
idx=hetTds.npoints*hetTds.nphase;
ju=hetTds.nphase-hetTds.nu;
js=hetTds.nphase-hetTds.ns;
YU = reshape(x1(idx+1:idx+ju*hetTds.nu,1),hetTds.nphase-hetTds.nu,hetTds.nu);
idx = idx + ju*hetTds.nu;
YS = reshape(x1(idx+1:idx+js*hetTds.ns,1),hetTds.nphase-hetTds.ns,hetTds.ns);
idx = idx + js*hetTds.ns;x1(end,1);
p(hetTds.ActiveParams) = p(ap);