hcli_jac.m

这是国外用的研究分岔的完整的M程序

for j=1:nb_par
  J(nml_n_1:nml+2*n,nml_n_1+j)=...
    J(nml_n_1:nml+2*n,nml_n_1+j)+sys_deri(xx2,par,[],free_par(j),[]);
end;


% eigenvalues and eigenvector equations

for k=1:s1
  res(nml+2*n+(k-1)*n+1:nml+2*n+k*n,1)=...
    eye(n)*lambda_v(k)*v(:,k)-sys_deri(xx1,par,0,[],[])*v(:,k);
  for i=1:nb_tau
    res(nml+2*n+(k-1)*n+1:nml+2*n+k*n,1)=...
      res(nml+2*n+(k-1)*n+1:nml+2*n+k*n,1)...
        -sys_deri(xx1,par,i,[],[])*exp(-lambda_v(k)*tau(i))*v(:,k);
  end;
end;

% entries for the v_k variables

for k=1:s1
  J(nml+2*n+(k-1)*n+1:nml+(k+2)*n,...
    nml_n_1+nb_par+2*n+(k-1)*n+1:nml_n_1+nb_par+2*n+k*n)=...
     J(nml+2*n+(k-1)*n+1:nml+2*n+k*n,...
       nml_n_1+nb_par+2*n+(k-1)*n+1:nml_n_1+nb_par+2*n+k*n)+...
         eye(n)*lambda_v(k)-sys_deri(xx1,par,0,[],[]);
  for i=1:nb_tau
    J(nml+2*n+1+(k-1)*n:nml+2*n+k*n,...
      nml_n_1+nb_par+2*n+(k-1)*n+1:nml_n_1+nb_par+2*n+k*n)=...
       J(nml+2*n+1+(k-1)*n:nml+2*n+k*n,...
         nml_n_1+nb_par+2*n+(k-1)*n+1:nml_n_1+nb_par+2*n+k*n)-...
           sys_deri(xx1,par,i,[],[])*exp(-lambda_v(k)*tau(i));
  end;
end; 

% entry for the lambda_v_k variables

for k=1:s1
  J(nml+2*n+(k-1)*n+1:nml+(k+2)*n,nml_n_1+nb_par+2*n+s1*n+k)=...
    J(nml+2*n+(k-1)*n+1:nml+(k+2)*n,nml_n_1+nb_par+2*n+s1*n+k)+v(:,k);
  for i=1:nb_tau
    J(nml+2*n+(k-1)*n+1:nml+(k+2)*n,nml_n_1+nb_par+2*n+s1*n+k)=...
      J(nml+2*n+(k-1)*n+1:nml+(k+2)*n,nml_n_1+nb_par+2*n+s1*n+k)+...
      sys_deri(xx1,par,i,[],[])*v(:,k)*tau(i)*exp(-lambda_v(k)*tau(i));
  end;
end;

% entries for the free parameters

for k=1:s1
  for j=1:nb_par
    J(nml+2*n+(k-1)*n+1:nml+(k+2)*n,nml_n_1+j)=...
       J(nml+2*n+(k-1)*n+1:nml+(k+2)*n,nml_n_1+j)-...
       sys_deri(xx1,par,0,free_par(j),[])*v(:,k);
    for i=1:nb_tau
      J(nml+2*n+(k-1)*n+1:nml+(k+2)*n,nml_n_1+j)=...
        J(nml+2*n+(k-1)*n+1:nml+(k+2)*n,nml_n_1+j)...
        -sys_deri(xx1,par,i,free_par(j),[])*v(:,k)*exp(-lambda_v(k)*tau(i));
      % special case if delay is free parameter  
      if free_par(j)==n_tau(i)
        J(nml+2*n+(k-1)*n+1:nml+(k+2)*n,nml_n_1+j)=...
          J(nml+2*n+(k-1)*n+1:nml+(k+2)*n,nml_n_1+j)...
          +lambda_v(k)*sys_deri(xx1,par,i,[],[])*v(:,k)*...
          exp(-lambda_v(k)*tau(i));
      end;
    end;
  end;
end;

% entry for the steady state point

for k=1:s1
  for i=0:nb_tau
    J(nml+2*n+(k-1)*n+1:nml+(2+k)*n,nml_n_1+nb_par+1:nml_n_1+nb_par+n)=...
      J(nml+2*n+(k-1)*n+1:nml+(2+k)*n,nml_n_1+nb_par+1:nml_n_1+nb_par+n)-...
      sys_deri(xx1,par,[0 i],[],v(:,k));
    for j=1:nb_tau
      J(nml+2*n+(k-1)*n+1:nml+(2+k)*n,nml_n_1+nb_par+1:nml_n_1+nb_par+n)=...
        J(nml+2*n+(k-1)*n+1:nml+(2+k)*n,nml_n_1+nb_par+1:nml_n_1+nb_par+n)-...
        sys_deri(xx1,par,[j i],[],v(:,k))*exp(-lambda_v(k)*tau(j));
    end;
  end;             
end;

% entries for the normalization condition: sum(vk_i^2=1);   

for k=1:s1
  res(nml+2*n+s1*n+k,1)=v(:,k)'*v(:,k)-1;
  J(nml+2*n+s1*n+k,nml_n_1+nb_par+2*n+(k-1)*n+1:nml_n_1+nb_par+2*n+k*n)=...
    2*v(:,k)';
end;

% characteristic equations for w_k

for k=1:s2
  res(nml+2*n+s1*(n+1)+(k-1)*n+1:nml+2*n+s1*(n+1)+k*n,1)=...
    lambda_w(k)'*w(:,k)-sys_deri(xx2,par,0,[],[])'*w(:,k);
  for i=1:nb_tau
    res(nml+2*n+s1*(n+1)+(k-1)*n+1:nml+2*n+s1*(n+1)+k*n,1)=...
      res(nml+2*n+s1*(n+1)+(k-1)*n+1:nml+2*n+s1*(n+1)+k*n,1)-...
      (sys_deri(xx2,par,i,[],[])*exp(-lambda_w(k)*tau(i)))'*w(:,k);
  end;
end;

% entries for the w_k variables

for k=1:s2
  J(nml+2*n+s1*(n+1)+(k-1)*n+1:nml+2*n+s1*(n+1)+k*n,...
    nml_n_1+nb_par+2*n+s1*(n+1)+(k-1)*n+1:nml_n_1+nb_par+2*n+s1*(n+1)+k*n)=...
    J(nml+2*n+s1*(n+1)+(k-1)*n+1:nml+2*n+s1*(n+1)+k*n,...
    nml_n_1+nb_par+2*n+s1*(n+1)+(k-1)*n+1:...
    nml_n_1+nb_par+2*n+s1*(n+1)+k*n)+...
    eye(n)*lambda_w(k)'-(sys_deri(xx2,par,0,[],[]))';

  for j=1:nb_tau
    J(nml+2*n+s1*(n+1)+(k-1)*n+1:nml+2*n+s1*(n+1)+k*n,...
      nml_n_1+nb_par+2*n+s1*(n+1)+(k-1)*n+1:...
      nml_n_1+nb_par+2*n+s1*(n+1)+k*n)=...
      J(nml+2*n+s1*(n+1)+(k-1)*n+1:nml+2*n+s1*(n+1)+k*n,...
      nml_n_1+nb_par+2*n+s1*(n+1)+(k-1)*n+1:...
      nml_n_1+nb_par+2*n+s1*(n+1)+k*n)-...
     (sys_deri(xx2,par,j,[],[])*exp(-lambda_w(k)*tau(j)))';
  end;
end; 

% entry for the lambda_w_k variables

for k=1:s2
  J(nml+2*n+s1*(n+1)+(k-1)*n+1:...
    nml+2*n+s1*(n+1)+k*n,nml_n_1+nb_par+2*n+s1*(n+1)+s2*n+k)=...
    J(nml+2*n+s1*(n+1)+(k-1)*n+1:...
    nml+(k+2)*n+s1*(n+1),nml_n_1+nb_par+2*n+s1*(n+1)+s2*n+k)+w(:,k);
  for j=1:nb_tau
    J(nml+2*n+s1*(n+1)+(k-1)*n+1:...
      nml+(k+2)*n+s1*(n+1),nml_n_1+nb_par+2*n+s1*(n+1)+s2*n+k)=...
      J(nml+2*n+s1*(n+1)+(k-1)*n+1:...
      nml+(k+2)*n+s1*(n+1),nml_n_1+nb_par+2*n+s1*(n+1)+s2*n+k)+...
      (sys_deri(xx2,par,j,[],[])*exp(-lambda_w(k)*tau(j)))'*w(:,k)*tau(j);
  end;
end;
% entries for the free parameters

for k=1:s2
  for j=1:nb_par
    J(nml+2*n+s1*(n+1)+(k-1)*n+1:nml+(2+s1+k)*n+s1,nml_n_1+j)=...
      J(nml+2*n+s1*(n+1)+(k-1)*n+1:nml+(2+s1+k)*n+s1,nml_n_1+j)-...
      (sys_deri(xx2,par,0,free_par(j),[]))'*w(:,k);
    for i=1:nb_tau
      J(nml+2*n+s1*(n+1)+(k-1)*n+1:nml+s1*(n+1)+(k+2)*n,nml_n_1+j)=...
        J(nml+2*n+s1*(n+1)+(k-1)*n+1:nml+(k+2)*n+s1*(n+1),nml_n_1+j)...
        -(sys_deri(xx2,par,i,free_par(j),[])*exp(-lambda_w(k)*tau(i)))'*w(:,k);
      % special case if delay is free parameter  
      if free_par(j)==n_tau(i)
        J(nml+2*n+s1*(n+1)+(k-1)*n+1:nml+2*n+s1*(n+1)+k*n,nml_n_1+j)=...
        J(nml+2*n+s1*(n+1)+(k-1)*n+1:nml+2*n+s1*(n+1)+k*n,nml_n_1+j)...
        +lambda_w(k)*(sys_deri(xx2,par,i,[],[])*exp(-lambda_w(k)*tau(i)))'*w(:,k);
      end;
    end;
  end;
end;

% entry for the steady state point
In=eye(n);
for k=1:s2
  for i=0:nb_tau
    for j=1:n
      J(nml+2*n+(k-1)*n+s1*(n+1)+j,nml_n_1+nb_par+n+(1:n))=...
        J(nml+2*n+(k-1)*n+s1*(n+1)+j,nml_n_1+nb_par+n+(1:n))-...
          w(:,k)'*sys_deri(xx2,par,[i 0],[],In(:,j));
    end;
    for j=1:nb_tau
      for l=1:n
        J(nml+2*n+(k-1)*n+s1*(n+1)+l,nml_n_1+nb_par+n+(1:n))=...
          J(nml+2*n+(k-1)*n+s1*(n+1)+l,nml_n_1+nb_par+n+(1:n))-...
          w(:,k)'*sys_deri(xx2,par,[i j],[],In(:,l))*...
          exp(-lambda_w(k)*tau(j));
      end ;
    end;
  end; 
end;            

% entries for the normalization condition: sum(wk_i^2=1);   

for k=1:s2
  res(nml+2*n+s1*(n+1)+s2*n+k,1)=w(:,k)'*w(:,k)-1;
  J(nml+2*n+s1*n+s1+s2*n+k,...
    nml_n_1+nb_par+2*n+s1*(n+1)+(k-1)*n+1:nml_n_1+nb_par+2*n+s1*(n+1)+k*n)=...
    2*w(:,k)';
end;

% integral orthogonality conditions

for k=1:s2
  res(nml+2*n+s1*n+s1+s2*n+s2+k,1)=w(:,k)'*(profile(:,ml+1)-x2);
end;

% first term of derivative to w_k 
for k=1:s2
  J(nml+2*n+s1*n+s1+s2*n+s2+k,...
    nml_n_1+nb_par+2*n+s1*n+s1+(k-1)*n+1:nml_n_1+nb_par+2*n+s1*n+s1+k*n)=...
    J(nml+2*n+s1*n+s1+s2*n+s2+k,...
    nml_n_1+nb_par+2*n+s1*n+s1+(k-1)*n+1:...
    nml_n_1+nb_par+2*n+s1*n+s1+k*n)+...
    (profile(:,ml+1)-x2)';
end;

% derivative to y(t_l)

for k=1:s2 
  J(nml+2*n+s1*n+s1+s2*n+s2+k,nml+1:nml+n)=w(:,k)';
end;

for l=1:nb_tau 
  if (1-tau(l)/T)~=1  % complicated here because of machine precision 
                      % problems with tau=0 and variable
    A=sys_deri(xx2,par,l,[],[]);
    % recall that gauss_c=poly_gau(m);

    tau_T=1-tau(l)/T;
    index=length(t)-m;
    %save tauT tau_T t T tau;
    while tau_T<t(index)
      index=index-m;
    end;
    K=index;
    gauss_c_trans(1:m)=gauss_c*(t(index+m)-tau_T)+tau_T;
    for i=1:m       
      Pa(index-K+i,:)=...
        poly_elg(m,(gauss_c_trans(index-K+i)-tau_T)/(t(index+m)-tau_T));
      gauss_y(1:n,index-K+i)=profile(:,index:index+m)*Pa(index-K+i,:)'-x2;
      gauss_weight(index-K+i)=gauss_abs(i)*T*(t(index+m)-tau_T);
    end;

    for index=K+m:m:ml-m+1
      % transform gauss_c to points in [t(index),t(index+m)]
      gauss_c_trans(index-K+1:index-K+m)=...
        gauss_c*(t(index+m)-t(index))+t(index);
      % evaluate the profile in those transformed gauss points 
      for i=1:m
        Pa(index-K+i,:)=poly_elg(m,(gauss_c_trans(index-K+i)-...
          t(index))/(t(index+m)-t(index)));
        gauss_y(1:n,index-K+i)=profile(:,index:index+m)*Pa(index-K+i,:)'-x2;
        gauss_weight(index-K+i)=gauss_abs(i)*T*(t(index+m)-t(index));
      end;
    end;
    theta=(gauss_c_trans*T-T+tau(l));
   
    for k=1:s2
      for i=1:length(gauss_weight)
        res(nml+2*n+s1*(n+1)+s2*n+s2+k,1)=...
          res(nml+2*n+s1*(n+1)+s2*n+s2+k,1)+...
          gauss_weight(i)*w(:,k)'*exp(-lambda_w(k)*theta(i))* ...
          A*gauss_y(:,i);
      end;   
    end;

    for k=1:s2
      % derivative to w_k
      for i=1:length(gauss_weight)
        J(nml+2*n+s1*n+s1+s2*n+s2+k,...
          nml_n_1+nb_par+2*n+s1*n+s1+(k-1)*n+1:...
          nml_n_1+nb_par+2*n+s1*n+s1+k*n)=...
          J(nml+2*n+s1*n+s1+s2*n+s2+k,...
          nml_n_1+nb_par+2*n+s1*n+s1+(k-1)*n+1:...
          nml_n_1+nb_par+2*n+s1*n+s1+k*n)+...
          gauss_weight(i)*exp(-lambda_w(k)*theta(i))*(A*gauss_y(:,i))';
      end;
    end;

    % derivative to lambda_w_k

    for k=1:s2
      for i=1:length(gauss_weight)
        J(nml+2*n+s1*n+s1+s2*n+s2+k,nml_n_1+nb_par+2*n+s1*n+s1+s2*n+k)=...
          J(nml+2*n+s1*n+s1+s2*n+s2+k,nml_n_1+nb_par+2*n+s1*n+s1+s2*n+k)-...
          gauss_weight(i)*w(:,k)'*theta(i)*...
          exp(-lambda_w(k)*theta(i))*A*gauss_y(:,i);
      end;  
    end;

    % derivative to x2
  
    for k=1:s2
      J(nml+2*n+s1*n+s1+s2*n+s2+k,nml_n_1+nb_par+n+1:...
        nml_n_1+nb_par+2*n)=-w(:,k)';
      for i=1:length(gauss_weight)
        dA1_dx2= sys_deri(xx2,par,[l 0],[],gauss_y(:,i));
        for j=1:nb_tau
          dA1_dx2= dA1_dx2+sys_deri(xx2,par,[l j],[],gauss_y(:,i));
        end;
        J(nml+2*n+s1*n+s1+s2*n+s2+k,nml_n_1+nb_par+n+1:nml_n_1+nb_par+2*n)=...
        J(nml+2*n+s1*n+s1+s2*n+s2+k,nml_n_1+nb_par+n+1:nml_n_1+nb_par+2*n)+...
          gauss_weight(i)*w(:,k)'*exp(-lambda_w(k)*theta(i))*(dA1_dx2-A);
      end;    
    end;

    % derivative to free parameters
 
    for j=1:nb_par
      dA1_deta = sys_deri(xx2,par,l,free_par(j),[]);
      for k=1:s2
        for i=1:length(gauss_weight)
          J(nml+2*n+s1*n+s1+s2*n+s2+k,nml_n_1+j)=...
            J(nml+2*n+s1*n+s1+s2*n+s2+k,nml_n_1+j)+...
            gauss_weight(i)*w(:,k)'*...
            exp(-lambda_w(k)*theta(i))*dA1_deta*gauss_y(:,i);
        end;   
      end;
    end;

    % derivative to T is zero

    % derivatives to representation points 
    for k=1:s2
      for i=1:length(gauss_weight)
        for teller=1:m+1
          J(nml+2*n+s1*n+s1+s2*n+s2+k,...
            (K-1+m*floor((i-1)/m)+teller-1)*n+(1:n))=...
            J(nml+2*n+s1*n+s1+s2*n+s2+k,...
            (K-1+m*floor((i-1)/m)+teller-1)*n+(1:n))+...
            gauss_weight(i)*w(:,k)'*exp(-lambda_w(k)*theta(i))*A*Pa(i,teller); 
        end;
      end;
    end; 

  end;

end;

% continuity condition in the initial point

res(nml+2*n+s1*n+s1+s2*n+2*s2+1:nml+2*n+s1*(n+1)+s2*n+2*s2+n,1)=...
  profile(:,1)-x1;
for k=1:s1
  res(nml+2*n+s1*n+s1+s2*n+2*s2+1:nml+2*n+s1*(n+1)+s2*n+2*s2+n,1)=...
    res(nml+2*n+s1*n+s1+s2*n+2*s2+1:nml+2*n+s1*(n+1)+s2*n+2*s2+n,1)-...
    epsilon*alpha(k)*v(:,k);
end;

% derivative to x1

J(nml+2*n+s1*(n+1)+s2*(n+2)+1:nml+2*n+s1*(n+1)+s2*(n+2)+n,...
  nml_n_1+nb_par+1:nml_n_1+nb_par+n)=-eye(n);

% derivative to alpha(k)

J(nml+2*n+s1*(n+1)+s2*(n+2)+1:nml+2*n+s1*(n+1)+s2*(n+2)+n,...
  nml_n_1+nb_par+2*n+s1*n+s1+s2*n+s2+1:...
  nml_n_1+nb_par+2*n+s1*n+s1+s2*n+s2+s1)=-epsilon*v;

% derivative to v(:,k)

for k=1:s1

  J(nml+2*n+s1*(n+1)+s2*(n+2)+1:nml+2*n+s1*(n+1)+s2*(n+2)+n,...
    nml_n_1+nb_par+2*n+(k-1)*n+1:nml_n_1+nb_par+2*n+k*n)=...
    -eye(n)*alpha(k)*epsilon;

end;

% derivative to y(0)

J(nml+2*n+s1*(n+1)+s2*(n+2)+(1:n),1:n)=eye(n);

% phase condition:

if ph & non_gauss,
  for l_i=1:l
    index_a=(l_i-1)*m+1;
    for k=1:m
      fac=gauss_abs(k)*(t((l_i-1)*m+1)-t(l_i*m+1));
      dPa=poly_dla(t(index_a:index_a+m),gauss_c(k));
      Pa=poly_lgr(t(index_a:index_a+m),gauss_c(k));
      u_prime_previous=previous.profile(:,index_a:index_a+m)*dPa';
      for q=1:m+1
        J(nml+2*n+s1*(n+1)+s2*(n+2)+n+1,...
          (l_i-1)*m*n+1+(q-1)*n:(l_i-1)*m*n+q*n)= ...
	  J(nml+2*n+s1*(n+1)+s2*(n+2)+n+1,...
          (l_i-1)*m*n+1+(q-1)*n:(l_i-1)*m*n+q*n) + ...
          fac*Pa(q)*u_prime_previous';
      end;
    end;
  end;
end;

if ph
  res(nml+2*n+s1*(n+1)+s2*(n+2)+n+1,1)=0;
end;

% extra condition on the alpha(k)'s (their squares have to sum
%  up to 1 -- with a fixed epsilon, this adds the unknown T)
res(nml+2*n+s1*n+s1+s2*n+s2+s2+n+2,1)=alpha'*alpha-1;
J(nml+2*n+s1*n+s1+s2*n+s2+s2+n+2,...
  nml_n_1+nb_par+2*n+s1*(n+1)+s2*(n+1)+1:...
    nml_n_1+nb_par+2*n+s1*(n+1)+s2*(n+1)+s1)=2*alpha';

return;