bvp_hett_jac.m

计算动力学系统的分岔图

% Jacobian of boundary value problem
%
% ============================================

function result = BVP_HetT_jac(x,p,YS,YU,J)

global hetTds

ups = reshape(x,hetTds.nphase,hetTds.npoints);
%p = num2cell(p);
% Allocate space for sparse jacobian
ks=hetTds.nphase*hetTds.npoints+(hetTds.nphase-hetTds.nu)*hetTds.nu+(hetTds.nphase-hetTds.ns)*hetTds.ns;
result = spalloc(ks,ks,ks/2);
n=hetTds.nphase;
N=hetTds.npoints;
x1=ups(:,1); xN=ups(:,N);
% Component 1 (the initial fixed point)
% ===========

A1=hetT_jac(x1,p,J);
result(1:n, 1:n) = A1-eye(n);
%Derivatives w.r.t active parameter
%jp=hetT_jacp(x1,p,J);
%result(1:n, ks+1) = jp;
%result(1,35),result(2,35),pause

% Component 2 (the iteration conditions)
% ===========
for j=3:N-1
 
result((j-2)*n+1:(j-1)*n, (j-2)*n+1:(j-1)*n)=hetT_jac(ups(:,j-1),p,J);
result((j-2)*n+1:(j-1)*n, (j-1)*n+1:j*n)=-eye(n);
end

%Derivatives w.r.t active parameter
%for j=3:N-1
%jp=hetT_jacp(ups(:,j-1),p,J);
%result((j-2)*n+1:(j-1)*n, ks+1)=jp;%jp(:,hetTds.ActiveParams);
%end
%result(:,35),pause

% Component 3 (the final fixed point)
% ===========

AN=hetT_jac(xN,p,J);
result((N-2)*n+1:(N-1)*n,(N-1)*n+1:N*n) = AN-eye(n);
%Derivatives w.r.t active parameter
%jp=hetT_jacp(xN,p,J);
%result((N-2)*n+1:(N-1)*n, ks+1) =jp;
%result(:,35),pause
%result(29:30,31:32),pause

% Component 5 % UNSTABLE
% Ricatti blocks from unstable eigenspace
% F(Y_U)=R22Y_U-Y_UR11+E21-Y_UR12Y_U;
% ===========

% D1=R_22Y_U
Q0U = hetTds.Q0;
[R11, R12, E21, R22] = het_RicattiCoeff(Q0U,A1,hetTds.nu);

% derivatives of D1 = R22 * YU w.r.t Y_U)ij
l=n*(N-1); h=N*n;
for j=1:n-hetTds.nu
    for i=1:hetTds.nu
         idx1=l+i+(j-1)*hetTds.nu;
         idx2=h+1+(i-1)*(n-hetTds.nu); 
         idx3=h+(n-hetTds.nu)+(i-1)*(n-hetTds.nu);
        result(idx1,idx2:idx3)=result(idx1,idx2:idx3)+R22(j,1:n-hetTds.nu);
    end
end


% derivatives of D2 = YU * R11 w.r.t Y_U)ij

for j=1:n-hetTds.nu
    for s=1:hetTds.nu
        idx1=l+1+(j-1)*hetTds.nu;
        idx2=l+hetTds.nu+(j-1)*hetTds.nu;
        idx3=h+j+(s-1)*(n-hetTds.nu);
       % result(idx1:idx2,idx3)=result(idx1:idx2,idx3)-(R11(1:hetTds.nu,i))';
       result(idx1:idx2,idx3)=result(idx1:idx2,idx3)-(R11(s,1:hetTds.nu))';
    end
end


% derivatives of D3 = YU * R12*YU w.r.t (Y_U)ij
D31=-hetTds.YU*R12;
for j=1:n-hetTds.nu
    for i=1:hetTds.nu
        idx1=l+i+(j-1)*hetTds.nu;
        idx2=h+1+(i-1)*(n-hetTds.nu);
        idx3=h+n-hetTds.nu+(i-1)*(n-hetTds.nu);
        result(idx1,idx2:idx3)=result(idx1,idx2:idx3)+D31(j,:);
    end
end
D32=-R12*hetTds.YU;
for j=1:n-hetTds.nu
    for s=1:hetTds.nu
        idx1=l+i+(j-1)*hetTds.nu;
        idx2=l+hetTds.nu+(j-1)*hetTds.nu;
        idx3=h+j+(s-1)*(n-hetTds.nu);
        result(idx1:idx2,idx3)=result(idx1:idx2,idx3)+D31(j,:);
    end
end
%result(31,33),pause

% Component 6
% derivatives of F(YU) w.r.t  x1 

hess=HetT_hess(x1,p,J);
l=n*(N-1); Q0=hetTds.Q0;

for i=1:n
    D=Q0'*hess(:,:,i)*Q0;
    D1=D(1:hetTds.nu,1:hetTds.nu);
    D2=D(1:hetTds.nu,hetTds.nu+1:n);
    D3=D(hetTds.nu+1:n,1:hetTds.nu);
    D4=D(hetTds.nu+1:n,hetTds.nu+1:n);
    for j=1:n-hetTds.nu
        for s=1:hetTds.nu;
            idx=l+s+(j-1)*hetTds.nu;
            result(idx,i)=result(idx,i)+D4(j,:)*hetTds.YU(:,s)-hetTds.YU(j,:)*D1(:,s)+D3(j,s);
            for k=1:hetTds.nu
                result(idx,i)=result(idx,i)-hetTds.YU(j,k)*(D2(k,:)*YU(:,s));
            end
        end           
    end
end
%result(31,1),result(31,2),pause

% Component 7
% derivatives of F(Y_U)=R22Y_U-Y_UR11+E21-Y_UR12Y_U; w.r.t  the active parameter 

% hessp=Het_hessp(x1,n2c(p),J);
% l=n*(N-1); Q0=hetTds.Q0;
% D=Q0'*hessp*Q0;
% D1=D(1:hetTds.nu,1:hetTds.nu);
% D2=D(1:hetTds.nu,1:hetTds.nu+1:n);
% D3=D(hetTds.nu+1:n,1:hetTds.nu);
% D4=D(hetTds.nu+1:n,hetTds.nu+1:n);
% for j=1:n-hetTds.nu
%     for s=1:hetTds.nu;
%         idx=l+s+(j-1)*hetTds.nu;
%         result(idx,ks+1)=result(idx,ks+1)+D4(j,:)*hetTds.YU(:,s)-hetTds.YU(j,:)*D1(:,s)+D3(j,s);
%         for k=1:hetTds.nu
%             result(idx,ks+1)=result(idx,ks+1)-hetTds.YU(j,k)*(D2(k,:)*YU(:,s));           
%         end           
%     end
% end
%result(:,35),pause

% Component  6% STABLE
% Ricatti blocks from stable eigenspace
% F(Y_S)=R22Y_S-Y_SR11+E21-YR12Y_S;
% =========
% D1=R_22Y_S
Q1S = hetTds.Q1;
[R11, R12, E21, R22] = het_RicattiCoeff(Q1S,AN,hetTds.ns);
% derivatives of D1 = R22 * Y_S w.r.t (YS)ij
l=n*(N-1)+(n-hetTds.nu)*hetTds.nu; h=N*n+(n-hetTds.nu)*hetTds.nu;
for j=1:n-hetTds.ns
    for i=1:hetTds.ns
        idx1=l+i+(j-1)*hetTds.ns;
        idx2=h+1+(i-1)*(n-hetTds.ns);
        idx3=h+(n-hetTds.ns)+(i-1)*(n-hetTds.ns);
        result(idx1,idx2:idx3)=result(idx1,idx2:idx3)+R22(j,:);
    end
end
% derivatives of D2= Y_S * R11 w.r.t (YU)ij

for j=1:n-hetTds.ns
    for s=1:hetTds.ns
        idx1=l+i+(j-1)*hetTds.ns;
        idx2=h+hetTds.nu+(j-1)*hetTds.ns;
        idx3=h+j+(s-1)*(n-hetTds.ns);
        result(idx1,idx2:idx3)=result(idx1,idx2:idx3)-(R11(s,:))';
    end
end
% derivatives of D3 = YS * R12*YS w.r.t (YS)ij
D31=-hetTds.YS*R12;
for j=1:n-hetTds.ns
    for i=1:hetTds.ns
        idx1=l+i+(j-1)*hetTds.ns;
        idx2=h+1+(i-1)*(n-hetTds.ns);
        idx3=h+n-hetTds.ns+(i-1)*(n-hetTds.ns);
        result(idx1,idx2:idx3)=result(idx1,idx2:idx3)+D31(j,:); 
    end
    
end

D32=-R12*hetTds.YS;
for j=1:n-hetTds.ns
    for s=1:hetTds.ns
        idx1=l+1+(j-1)*hetTds.ns;
        idx2=l+hetTds.ns+(j-1)*hetTds.ns;
        idx3=h+j+(s-1)*(n-hetTds.ns);
        result(idx1:idx2,idx3)=result(idx1:idx2,idx3)+(D32(s,:))';
        
    end
end
%result(32,34),pause

% Component 7
% derivatives of F(Y_S) w.r.t  xN
hess=HetT_hess(xN,p,J);
l=n*(N-1)+(n-hetTds.nu)*hetTds.nu;
h=n*(N-1);
hess=HetT_hess(xN,p,J);
for i=1:n
    D=Q1S'*hess(:,:,i)*Q1S;
    D1=D(1:hetTds.ns,1:hetTds.ns);
    D2=D(1:hetTds.ns,1:hetTds.ns+1:n);
    D3=D(hetTds.ns+1:n,1:hetTds.ns);
    D4=D(hetTds.ns+1:n,hetTds.ns+1:n);
    for j=1:n-hetTds.ns
        for s=1:hetTds.ns;
            idx=l+s+(j-1)*hetTds.ns;
            result(idx,h+i)=result(idx,h+i)+D4(j,:)*hetTds.YS(:,s)-hetTds.YS(j,:)*D1(:,s)+D3(j,s);
            for k=1:hetTds.ns
                result(idx,h+i)=result(idx,h+i)-hetTds.YS(j,k)*(D2(k,:)*YS(:,s));
            end
        end           
    end
end
%result(32,31:32),pause
% Derivatives of F(Y_S)=R22Y_S-Y_SR11+E21-Y_SR12Y_S w.r.t the active parameter 

% hessp=Het_hessp(xN,p,J); 
% l=n*(N-1)+hetTds.nu*(n-hetTds.nu);
% Q1=hetTds.Q1;
% D=Q1'*hessp*Q1;
% D1=D(1:hetTds.ns,1:hetTds.ns);
% D2=D(1:hetTds.ns,1:hetTds.ns+1:n);
% D3=D(hetTds.ns+1:n,1:hetTds.ns);
% D4=D(hetTds.ns+1:n,hetTds.ns+1:n);
% for j=1:n-hetTds.ns
%     for s=1:hetTds.ns;
%         idx=l+s+(j-1)*hetTds.ns;
%         result(idx,ks+1)=D4(j,:)*hetTds.YS(:,s)-hetTds.YS(j,:)*D1(:,s)+D3(j,s);
%         for k=1:hetTds.nu
%             result(idx,ks+1)=result(idx,ks+1)-hetTds.YS(j,k)*(D2(k,:)*YS(:,s));    
%         end           
%     end
% end
%result,pause


% Component  7 (First vectors along unstable eigenspaces)
% ===========
% derivatives w.r.t x1 
    Q0U = hetTds.Q0;
    vect = ups(:,2) - x1;
    QU =Q0U*[-YU'; eye(size(YU,1))];
    l=n*(N-1)+(n-hetTds.nu)*hetTds.nu+(n-hetTds.ns)*hetTds.ns;
    for i=1:n-hetTds.nu
       result(l+i, 1:n )=-QU(:,i)';
    end
   %result(33,1:2),pause
    
    
   % derivatives w.r.t x2
    vect = ups(:,2) - x1;
    l=n*(N-1)+(n-hetTds.nu)*hetTds.nu+(n-hetTds.ns)*hetTds.ns;
    for i=1:n-hetTds.nu
       result(l+i, n+1:2*n )=QU(:,i)';
    end
   
    %result(33,3:4),pause
    
 % derivatives w.r.t  components of YU_{(nu+i)}   
 H=vect'*Q0U;h=n*N;h=n*N;H=H(hetTds.nu); 
 
 for i=1:1:n-hetTds.nu
     idx1=h+1+(n-hetTds.nu)*(i-1);
     idx2=h+hetTds.nu+(n-hetTds.nu)*(i-1);
     result(l+i,idx1:idx2)=result(l+i,idx1:idx2)-H;
 end

%result(33,33),pause

 % Component  8 (Last vectors along stable eigenspaces)
% ===========
% derivatives w.r.t xN-1

    Q1S = hetTds.Q1;
    vect = ups(:,N-1) - xN;
    QS =Q1S*[-YS'; eye(size(YS,1))];
    l=n*(N-1)+(n-hetTds.nu)*hetTds.nu+(n-hetTds.ns)*hetTds.ns+hetTds.ns;
    for i=1:n-hetTds.ns
       result(l+i, n*(N-2)+1:n*(N-1))=result(l+i, n*(N-2)+1:n*(N-1))+QS(:,i)'; 
    end
   
    
   % derivatives w.r.t xN    
    l=n*(N-1)+(n-hetTds.nu)*hetTds.nu+(n-hetTds.ns)*hetTds.ns+hetTds.ns;
    for i=1:n-hetTds.nu
        
       result(l+i, n*(N-1)+1:N*n )=-QS(:,i)';  %reza20
    end
     %result(34,29:32),pause
    
    % derivatives w.r.t  components of YU_{(nu+i)}
     vect = ups(:,N-1) - xN;
     H=vect'*Q1S;h=n*N+(n-hetTds.nu)*hetTds.nu;
     for i=1:1:n-hetTds.ns
      for j=1:hetTds.ns
        result(l+i,h+j+(n-hetTds.ns)*(i-1))=-H(j);
      end
     end
    
%    det(result(:,1:106)),pause