bvp_hom_jac.asv

计算动力学系统的分岔图

% Buihomds up jacobian of boundary value problem
%
% ============================================

function result = BVP_Het_jac(mapfile,x,p,YS,YU,J)

global hetds
ups = reshape(x,hetds.nphase,hetds.npoints);
p = num2cell(p);

% Allocate space for sparse jacobian
k=hetds.nphase*hetds.npoints+2*hetds.nu*hetds.ns;
result = spalloc(k,k,k/3);
n=hetds.nphase; N=hetds.npoints;
% Component 1 (the initial fixed point)
% ===========
result(1:n, 1:n) = hetjac(ups(:,1),p)-eye(n);
% Component 2 (the iteration conditions)
% ===========
for j=3:N-1
result((j-2)*n+1:(j-1)*n, n(j-2)*n+1:(j-1)*n)=hetjac(ups(j-1),p);
result((j-2)*n+1:(j-1)*n, n(j-1)*n+1:j*n)=-eye(n);
end
% Component 3 (the final fixed point)
% ===========
result((N-2)*n+1:(N-1)n,(N-2)*n+1:N*n) = hetjac(ups(:,N),p)-eye(n);

% Component 4 (Unstable eigenspaces)
% ===========
% D1=R_22Y_U
Q0U = homds.oldUnstableQ;
Q0S = homds.oldStableQ;
hess = Hom_odehess(x0,p);
hessp = Hom_odehessp(x0,p);
for j=1:homds.nphase+2
    if j<=homds.nphase
        Atemp = hess(:,:,j);
    else
        Atemp = hessp(:,:,j-homds.nphase);
    end    
    [U11, U12, UE21, U22] = ricattiCoeff(Q0U,Atemp,homds.npos);
    tmp = (U22*YU - YU*U11 + UE21 - YU*U12*YU)';
    for i=1:homds.nneg
        result(homds.ncoords+phasepresent+1:homds.ncoords+phasepresent+homds.Ysize,homds.ncoords+j) = tmp(:,i);
    end    
    [S11, S12, SE21, S22] = ricattiCoeff(Q0S,Atemp,homds.nneg);
    tmp = (S22*YS - YS*S11 + SE21 - YS*S12*YS)';
    for i=1:homds.npos
        result(homds.ncoords+phasepresent+homds.Ysize+1:homds.ncoords+phasepresent+2*homds.Ysize,homds.ncoords+j) = tmp(:,i);
    end
end

% UNSTABLE
% ricatti blocks from unstable eigenspace
A=Hom_odejac(x0,p);
derU = [];
[U11, U12, UE21, U22] = ricattiCoeff(Q0U,A,homds.npos);
U12Y = U12 * YU;
YU12 = YU * U12;
for j=1:homds.npos
    for i=1:homds.nneg
        myres = sparse(homds.nneg,homds.npos,0);
        myres(:,j) = myres(:,j) + U22(:,i);
        myres(i,:) = myres(i,:) - U11(j,:);        
        myres(i,:) = myres(i,:) - U12Y(j,:);
        myres(:,j) = myres(:,j) - YU12(:,i);        
        % Derivative is computed to YU(i,j) i=1:homds.nneg j=1:homds.npos
        % must come at column (j-1)*(homds.nneg)+i, 
        % in form T and col-first == row-first
        rows = size(myres,1);
        cols = size(myres,2);
        for k=1:rows
            derU((k-1)*cols+1:k*cols,(j-1)*(homds.nneg)+i) = myres(k,:)';
        end
    end
end

% STABLE
derS = [];
[S11, S12, SE21, S22] = ricattiCoeff(Q0S,A,homds.nneg);
S12Y = S12 * YS;
YS12 = YS * S12;
for i=1:homds.npos
    for j=1:homds.nneg
        myres = sparse(homds.npos,homds.nneg,0);
        myres(:,j) = myres(:,j) + S22(:,i);
        myres(i,:) = myres(i,:) - S11(j,:);
        myres(i,:) = myres(i,:) - S12Y(j,:);
        myres(:,j) = myres(:,j) - YS12(:,i);        
        % Derivative is computed to YS(i,j)
        % must come at column (j-1)*(homds.npos)+i, 
        % in form T and col-first == row-first
        rows = size(myres,1);
        cols = size(myres,2);
        for k=1:rows
            derS((k-1)*cols+1:k*cols,(j-1)*(homds.npos)+i) = myres(k,:)';
        end
    end
end
result(homds.ncoords+phasepresent+1:homds.ncoords+phasepresent+homds.Ysize,end-2*homds.Ysize+1:end-homds.Ysize) = derU;
result(homds.ncoords+phasepresent+homds.Ysize+1:homds.ncoords+phasepresent+2*homds.Ysize,end-homds.Ysize+1:end) = derS;
end

% Component 5 (Last and first vectors along stable and unstable eigenspaces)
% ===========
if homds.npos
    indxs = [1:homds.nphase   homds.ncoords+1:homds.ncoords+homds.nphase];
    vect = ups(:,1) - x0;
    Q1U = Q0U * [-YU'; eye(size(YU,1))];
    vQ = -vect' * Q0U;
    for i=1:homds.nneg
        result(end-2-homds.nphase+i,indxs) = [Q1U(:,end-i+1)'     -Q1U(:,end-i+1)'];
    end
    myeye = zeros(homds.nneg);
    for j=1:homds.nneg
        myeye(j,homds.nneg-j+1) = 1;
    end
    for j=1:homds.npos
        result(end-2-homds.nphase+1:end-2-homds.nphase+homds.nneg,end-2*homds.Ysize+(j-1)*(homds.nneg)+1:end-2*homds.Ysize+j*(homds.nneg)) = vQ(:,j) * myeye;
    end
end

if homds.nneg
    indxs = [homds.ncoords-homds.nphase+1:homds.ncoords+homds.nphase];
    vect = ups(:,end) - x0;
    Q1S = Q0S * [-YS'; eye(size(YS,1))];
    vQ = -vect' * Q0S;
    for i=1:homds.npos
        result(end-2-homds.nphase+homds.nneg+i,indxs) = [Q1S(:,end-i+1)'     -Q1S(:,end-i+1)'];
    end
    myeye = zeros(homds.npos);
    for j=1:homds.npos
        myeye(j,homds.npos-j+1) = 1;
    end
    for j=1:homds.nneg
        result(end-2-homds.nphase+homds.nneg+1:end-2,end-homds.Ysize+(j-1)*homds.npos+1:end-homds.Ysize+j*homds.npos) = vQ(:,j) * myeye;
    end
end

% Component 6 (Distances from endpoints to equilibrium equal to epsilons)
% ===========
    indxs = [homds.phases   homds.ncoords+1:homds.ncoords+homds.nphase];
    val = ups(:,1) - x0;
    val = val' / norm(val);
    result(end-1,indxs) = [val   -val];
    if homds.extravec(2)
        result(end-1,homds.ncoords+homds.nphase+3+homds.extravec(1)) = -1;
    end
    indxs = [homds.ncoords-homds.nphase+1:homds.ncoords+homds.nphase];
    val = ups(:,end) - x0;
    val = val' / norm(val);
    result(end,indxs) = [val   -val];
    if homds.extravec(3)
        result(end,homds.ncoords+homds.nphase+3+sum(homds.extravec(1:2))) = -1;
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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);