bvp_het_jac.m

计算动力学系统的分岔图

% 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
ks=hetds.nphase*hetds.npoints+(hetds.nphase-hetds.nu)*hetds.nu+(hetds.nphase-hetds.ns)*hetds.ns;
result = spalloc(ks,ks+1,ks/2);
n=hetds.nphase;
N=hetds.npoints;
x1=ups(:,1); xN=ups(:,N);
% Component 1 (the initial fixed point)
% ===========

A1=hetjac(x1,p,J);
result(1:n, 1:n) = A1-eye(n);
%Derivatives w.r.t active parameter
jp=hetjacp(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)=hetjac(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=hetjacp(ups(:,j-1),p,J);
result((j-2)*n+1:(j-1)*n, ks+1)=jp;%jp(:,hetds.ActiveParams);
end
%result(:,35),pause

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

AN=hetjac(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=hetjacp(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 = hetds.Q0;
[R11, R12, E21, R22] = het_RicattiCoeff(Q0U,A1,hetds.nu);

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


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

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


% derivatives of D3 = YU * R12*YU w.r.t (Y_U)ij
D31=-hetds.YU*R12;
for j=1:n-hetds.nu
    for i=1:hetds.nu
        idx1=l+i+(j-1)*hetds.nu;
        idx2=h+1+(i-1)*(n-hetds.nu);
        idx3=h+n-hetds.nu+(i-1)*(n-hetds.nu);
        result(idx1,idx2:idx3)=result(idx1,idx2:idx3)+D31(j,:);
    end
end
D32=-R12*hetds.YU;
for j=1:n-hetds.nu
    for s=1:hetds.nu
        idx1=l+i+(j-1)*hetds.nu;
        idx2=l+hetds.nu+(j-1)*hetds.nu;
        idx3=h+j+(s-1)*(n-hetds.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=Het_hess(x1,p,J);
l=n*(N-1); Q0=hetds.Q0;

for i=1:n
    D=Q0'*hess(:,:,i)*Q0;
    D1=D(1:hetds.nu,1:hetds.nu);
    D2=D(1:hetds.nu,hetds.nu+1:n);
    D3=D(hetds.nu+1:n,1:hetds.nu);
    D4=D(hetds.nu+1:n,hetds.nu+1:n);
    for j=1:n-hetds.nu
        for s=1:hetds.nu;
            idx=l+s+(j-1)*hetds.nu;
            result(idx,i)=result(idx,i)+D4(j,:)*hetds.YU(:,s)-hetds.YU(j,:)*D1(:,s)+D3(j,s);
            for k=1:hetds.nu
                result(idx,i)=result(idx,i)-hetds.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,p,J);
l=n*(N-1); Q0=hetds.Q0;
D=Q0'*hessp*Q0;
D1=D(1:hetds.nu,1:hetds.nu);
D2=D(1:hetds.nu,1:hetds.nu+1:n);
D3=D(hetds.nu+1:n,1:hetds.nu);
D4=D(hetds.nu+1:n,hetds.nu+1:n);
for j=1:n-hetds.nu
    for s=1:hetds.nu;
        idx=l+s+(j-1)*hetds.nu;
        result(idx,ks+1)=result(idx,ks+1)+D4(j,:)*hetds.YU(:,s)-hetds.YU(j,:)*D1(:,s)+D3(j,s);
        for k=1:hetds.nu
            result(idx,ks+1)=result(idx,ks+1)-hetds.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 = hetds.Q1;
[R11, R12, E21, R22] = het_RicattiCoeff(Q1S,AN,hetds.ns);
% derivatives of D1 = R22 * Y_S w.r.t (YS)ij
l=n*(N-1)+(n-hetds.nu)*hetds.nu; h=N*n+(n-hetds.nu)*hetds.nu;
for j=1:n-hetds.ns
    for i=1:hetds.ns
        idx1=l+i+(j-1)*hetds.ns;
        idx2=h+1+(i-1)*(n-hetds.ns);
        idx3=h+(n-hetds.ns)+(i-1)*(n-hetds.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-hetds.ns
    for s=1:hetds.ns
        idx1=l+i+(j-1)*hetds.ns;
        idx2=h+hetds.nu+(j-1)*hetds.ns;
        idx3=h+j+(s-1)*(n-hetds.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=-hetds.YS*R12;
for j=1:n-hetds.ns
    for i=1:hetds.ns
        idx1=l+i+(j-1)*hetds.ns;
        idx2=h+1+(i-1)*(n-hetds.ns);
        idx3=h+n-hetds.ns+(i-1)*(n-hetds.ns);
        result(idx1,idx2:idx3)=result(idx1,idx2:idx3)+D31(j,:); 
    end
    
end

D32=-R12*hetds.YS;
for j=1:n-hetds.ns
    for s=1:hetds.ns
        idx1=l+1+(j-1)*hetds.ns;
        idx2=l+hetds.ns+(j-1)*hetds.ns;
        idx3=h+j+(s-1)*(n-hetds.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=Het_hess(xN,p,J);
l=n*(N-1)+(n-hetds.nu)*hetds.nu;
h=n*(N-1);
hess=Het_hess(xN,p,J);
for i=1:n
    D=Q1S'*hess(:,:,i)*Q1S;
    D1=D(1:hetds.ns,1:hetds.ns);
    D2=D(1:hetds.ns,1:hetds.ns+1:n);
    D3=D(hetds.ns+1:n,1:hetds.ns);
    D4=D(hetds.ns+1:n,hetds.ns+1:n);
    for j=1:n-hetds.ns
        for s=1:hetds.ns;
            idx=l+s+(j-1)*hetds.ns;
            result(idx,h+i)=result(idx,h+i)+D4(j,:)*hetds.YS(:,s)-hetds.YS(j,:)*D1(:,s)+D3(j,s);
            for k=1:hetds.ns
                result(idx,h+i)=result(idx,h+i)-hetds.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)+hetds.nu*(n-hetds.nu);
Q1=hetds.Q1;
D=Q1'*hessp*Q1;
D1=D(1:hetds.ns,1:hetds.ns);
D2=D(1:hetds.ns,1:hetds.ns+1:n);
D3=D(hetds.ns+1:n,1:hetds.ns);
D4=D(hetds.ns+1:n,hetds.ns+1:n);
for j=1:n-hetds.ns
    for s=1:hetds.ns;
        idx=l+s+(j-1)*hetds.ns;
        result(idx,ks+1)=D4(j,:)*hetds.YS(:,s)-hetds.YS(j,:)*D1(:,s)+D3(j,s);
        for k=1:hetds.nu
            result(idx,ks+1)=result(idx,ks+1)-hetds.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 = hetds.Q0;
    vect = ups(:,2) - x1;
    QU =Q0U*[-YU'; eye(size(YU,1))];
    l=n*(N-1)+(n-hetds.nu)*hetds.nu+(n-hetds.ns)*hetds.ns;
    for i=1:n-hetds.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-hetds.nu)*hetds.nu+(n-hetds.ns)*hetds.ns;
    for i=1:n-hetds.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(hetds.nu); 
 
 for i=1:1:n-hetds.nu
     idx1=h+1+(n-hetds.nu)*(i-1);
     idx2=h+hetds.nu+(n-hetds.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 = hetds.Q1;
    vect = ups(:,N-1) - xN;
    QS =Q1S*[-YS'; eye(size(YS,1))];
    l=n*(N-1)+(n-hetds.nu)*hetds.nu+(n-hetds.ns)*hetds.ns+hetds.ns;
    for i=1:n-hetds.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-hetds.nu)*hetds.nu+(n-hetds.ns)*hetds.ns+hetds.ns;
    for i=1:n-hetds.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-hetds.nu)*hetds.nu;
     for i=1:1:n-hetds.ns
      for j=1:hetds.ns
        result(l+i,h+j+(n-hetds.ns)*(i-1))=-H(j);
      end
     end
    
%    det(result(:,1:106)),pause