hcli_jac.m

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

function [J,res]=hcli_jac(c,T,profile,t,deg,par,free_par,ph,lambda_v,...
                               lambda_w,v,w,alpha,epsilon,x1,x2,previous)

% function [J,res]=hcli_jac(c,T,profile,t,deg,par,free_par,ph,lambda_v,...
%                               lambda_w,v,w,alpha,epsilon,x1,x2,previous)
% INPUT:
%	c collocation parameters in [0,1]^m
%	T period 
%	profile profile in R^(n x m*l+1)
%	t representation points in [0,1]^(m*l+1)
%	deg degree piecewise polynomial
%	par current parameter values in R^p
%	free_par free parameters numbers in N^d 
%	ph use phase condition or not (1 or 0)
%       lambda_v unstable eigenvalues of x1 (in R^s1)
%       lambda_w unstable eigenvalues of x2 (in R^s2)
%       v unstable eigenvectors of x1 (in R^n x s1)
%       w unstable eigenvectors of x2 (in R^n x s2)
%       alpha coefficients of initial functionsegment (in R^s1)
%       epsilon global coefficient of initial function segment (in R)
%       x1 steady state solution at t=-inf (in R^n)
%       x2 steady state solution at t=+inf (in R^n)
%       (optional) previous: previous solution, used in phase condition
% OUTPUT: 
%	J   jacobian in
%            R^(n*m*l+3*n+(s1+s2)*n+s1+2*s2+1 x n*m*l+3*n+(s1+s2)*n+2*s1+s2+d)
%	res residual in R^(n*m*l+3*n+(s1+s2)*n+s1+2*s2+1)

% (c) DDE-BIFTOOL v. 2.02, 16/6/2002 

% initializiation of dimensions etc.

if ~exist('previous') | isempty(previous) 
  previous.kind='hcli';
  previous.parameter=par;
  previous.mesh=t;
  previous.degree=deg;
  previous.profile=profile;
  previous.period=T;
  previous.x1=x1;
  previous.x2=x2;
  previous.lambda_v=lambda_v;
  previous.lambda_w=lambda_w;
  previous.v=v;
  previous.w=w;
  previous.alpha=alpha;
  previous.epsilon=epsilon;
end;

if ((length(previous.profile)-1)/deg==floor((length(previous.profile)-1)/deg))
  previous.profile=hcli_eva(previous.profile,previous.mesh,t,deg);
  previous.mesh=t;
else 
  err=([previous.degree deg]);
  error('HCLI_JAC: previous and current solution are not of same degree');
end;

m=deg;
n=size(profile,1); % system dimension

tp_del=nargin('sys_tau');
if tp_del==0,
  n_tau=sys_tau; % delay numbers
  tau=par(n_tau); % delay values
  nb_tau=length(n_tau); % number of delays
else
  error('HCLI_JAC: computing connected orbits is not implemented for equations with state-dependent delays');
end;

l=(length(t)-1)/m; % number of intervals
s1=length(lambda_v); % number of unstable eigenvalues of x1
s2=length(lambda_w); % idem x2
nb_par=length(free_par); % number of free parameters

% check:

if l~=floor(l)
  err=[length(t) m],
  error('HCLI_JAC: t does not contain l intervals of m points!');
end;

if length(c)~=m & ~isempty(c)
  err=[length(c) m]
  error('HCLI_JAC: wrong number of collocation parameters!');
end;

if T<0,
  err=[T],
  error('HCLI_JAC: period became smaller than zero!');
end;

if max(tau)>T,
  err=[max(tau) T]
  error('HCLI_JAC: period became smaller than maximal delay!');
end;

% init J, res:

nml=n*m*l;
ml=m*l;
nml_n_1=nml+n+1;

J=zeros(nml+3*n+(s1+s2)*n+s1+2*s2+1+1,nml_n_1+2*n+(s1+s2)*n+2*s1+s2+nb_par);
res=zeros(nml+3*n+(s1+s2)*n+s1+2*s2+1+1,1);

% initialize gauss points 

if isempty(c)
  c=poly_gau(m);
  gauss_c=c;
  non_gauss=0;
else
  gauss_c=poly_gau(m);
  non_gauss=1;
end;

% determination of gaussian weights:

gauss_abs=ones(1,m);
g=poly_gau(m-1);
for k=1:m
  for j=1:m-1
    gauss_abs(k)=gauss_abs(k)/(gauss_c(k)-g(j));
  end;
  for j=1:m
    if j~=k
      gauss_abs(k)=gauss_abs(k)/(gauss_c(k)-gauss_c(j));
    end;
  end;
end;
gauss_abs=gauss_abs/sum(gauss_abs);

% more initialisation:

for m_i=1:m
  all_dPa(m_i,:)=poly_dla((0:m)/m,c(m_i));
  all_Pa(m_i,:)=poly_lgr((0:m)/m,c(m_i));
end;

% for all collocation points, make equation:

tTT=tau/(T*T);

for l_i=1:l
  % determine index for a in profile:          
  index_a=(l_i-1)*m+1;
  i_index_a=(index_a-1)*n;
  t_start=t(index_a);
  hhh=t(index_a+m)-t_start;
  for m_i=1:m
    i=index_a+m_i-1;
    i_range=(i-1)*n+1:i*n;
    
    % determine c  
    
    col=t_start+c(m_i)*hhh;

    % determine dPa for a:

    dPa=all_dPa(m_i,:)/hhh;

    % add dPa for da in J:
    
    for k=0:m,
      kk=k*n;
      J(i_range,i_index_a+kk+1:i_index_a+kk+n)= ...
        J(i_range,i_index_a+kk+1:i_index_a+kk+n)+eye(n)*dPa(k+1);
    end;

    % add sum a*dPa to res:

    u_prime=profile(:,index_a:index_a+m)*dPa';
    u_prime_previous=previous.profile(:,index_a:index_a+m)*dPa';
    
    res(i_range,1)=u_prime;

    % determine Pa for a:

    Pa=all_Pa(m_i,:); 

    % phase_condition:
    
    if ph & ~non_gauss
      fup=gauss_abs(m_i)*hhh*u_prime_previous';
      i_l_i=(l_i-1)*m*n;
      for q=0:m
        qq=q*n;
        J(nml+2*n+s1*(n+1)+s2*(n+2)+n+1,i_l_i+1+qq:i_l_i+qq+n)= ...
          J(nml+2*n+s1*(n+1)+s2*(n+2)+n+1,i_l_i+1+qq:i_l_i+qq+n) + Pa(q+1)*fup;
      end;
    end;

    % compute x:

    x=profile(:,index_a:index_a+m)*Pa';

    for tau_i=1:nb_tau

      % determine c_tau:
      c_tau_i=col-tau(tau_i)/T;
      c_tau(tau_i)=c_tau_i;  
      if c_tau_i<0 
      initial_function_segment(tau_i)=1;
      else 
      initial_function_segment(tau_i)=0;
      end;
 
      if (~initial_function_segment(tau_i))
        
        % determine index for b in profile:

        index_b_i=length(t)-m;
        while (c_tau_i<t(index_b_i))
          index_b_i=index_b_i-m;
        end;
   
        index_b(tau_i)=index_b_i;

        % c transformed to [0,1] and hhh_tau
     
        hhh_tau(tau_i)=t(index_b_i+m)-t(index_b_i);
        c_tau_trans(tau_i)=(c_tau_i-t(index_b_i))/hhh_tau(tau_i);

        % determine Pb for b:

        Pb(tau_i,:)=poly_elg(m,c_tau_trans(tau_i));

        % compute x_tau:
  
        x_tau(:,tau_i)=profile(:,index_b_i:index_b_i+m)*Pb(tau_i,:)';
      else 
        x_tau(:,tau_i)=x1+epsilon*v*(alpha.*exp(T*lambda_v*c_tau_i));
      end;

    end; % for tau_i=1:d

    % determine A0:
  
    T_A0=T*sys_deri([x x_tau],par,0,[],[]);

    % add -T*A0*Pa for da in J:

    for k=0:m
      kk=k*n; 
      J(i_range,i_index_a+kk+1:i_index_a+kk+n)= ...
        J(i_range,i_index_a+kk+1:i_index_a+kk+n) ...
          -T_A0*Pa(k+1);  
    end;

    % determine parameter derivatives:

    for p_i=1:nb_par
      df=sys_deri([x x_tau],par,[],free_par(p_i),[]);
      J(i_range,nml_n_1+p_i)=-T*df; 
    end;  

    % compute f(x,x_tau):

    f=sys_rhs([x x_tau],par);

    % add -f for dT in J:

    J(i_range,nml_n_1)=J(i_range,nml_n_1)-f;

    % add Tf in res:

    res(i_range,1)=res(i_range,1)-T*f;

    for t_i=1:nb_tau
  
      % determine A1:

      T_A1=T*sys_deri([x x_tau],par,t_i,[],[]);  
      
      if (~initial_function_segment(t_i))
         
        index_b(t_i);
        % add -T*A1*Pb for db in J:

        i_index_b=(index_b(t_i)-1)*n;
        for k=0:m
          kk=n*k;
          J(i_range,i_index_b+kk+1:i_index_b+kk+n)= ...
            J(i_range,i_index_b+kk+1:i_index_b+kk+n) - T_A1*Pb(t_i,k+1);  
        end;

        % determine dPb for b:
 
        dPb=poly_del(m,c_tau_trans(t_i))/hhh_tau(t_i);

        % add -T*A1*sum b*dP*dc_tau for dT in J:
          
        J(i_range,nml_n_1)=J(i_range,nml_n_1)+ ...
          -T_A1*(profile(:,index_b(t_i):index_b(t_i)+m)*dPb')*tTT(t_i);

        % special case where delay is the free parameter
        for p_i=1:length(free_par)
          if free_par(p_i)==n_tau(t_i)  
            J(i_range,nml_n_1+p_i)=J(i_range,nml_n_1+p_i) ...
              +T_A1/T*(profile(:,index_b(t_i):index_b(t_i)+m)*dPb');  
          end;
        end;
 
      else
        % derivatives to alpha_k

        for k=1:s1
          J(i_range,nml_n_1+nb_par+2*n+s1*(n+1)+s2*(n+1)+k)=...
            J(i_range,nml_n_1+nb_par+2*n+s1*(n+1)+s2*(n+1)+k)-...
              T_A1*epsilon*v(:,k)*exp(T*lambda_v(k)*c_tau(t_i));
        end;

        % derivatives to x1

        J(i_range,nml_n_1+nb_par+1:nml_n_1+nb_par+n)=...
        J(i_range,nml_n_1+nb_par+1:nml_n_1+nb_par+n)-T_A1;
      
        % derivatives to v(k)

        for k=1:s1
          J(i_range,nml_n_1+nb_par+2*n+(k-1)*n+1:nml_n_1+nb_par+2*n+k*n)=...
          J(i_range,nml_n_1+nb_par+2*n+(k-1)*n+1:nml_n_1+nb_par+2*n+k*n)-...
            T_A1*alpha(k)*epsilon*exp(T*lambda_v(k)*c_tau(t_i));
        end;

        % derivatives to lambda_v(k)
         
        for k=1:s1
          J(i_range,nml_n_1+nb_par+2*n+s1*n+k)=...
            J(i_range,nml_n_1+nb_par+2*n+s1*n+k)...
            -T_A1*alpha(k)*v(:,k)*epsilon*...
               T*c_tau(t_i)*exp(T*lambda_v(k)*c_tau(t_i)); 
        end;

        % derivative to T
        for k=1:s1
          J(i_range,nml_n_1)=J(i_range,nml_n_1)-...
            T_A1*alpha(k)*v(:,k)*exp(T*lambda_v(k)*c_tau(t_i))*...
             lambda_v(k)*epsilon*(c_tau(t_i)+tau(t_i)/T);
        end;

        % special case where delay is a free parameter
        for p_i=1:length(free_par)
          if free_par(p_i)==n_tau(t_i)  
            for k=1:s1
              J(i_range,nml_n_1+p_i)=J(i_range,nml_n_1+p_i)...
                +T_A1*epsilon*...
                v(:,k)*alpha(k)*exp(T*lambda_v(k)*c_tau(t_i))*lambda_v(k);
            end;
          end;
        end;

      end; % if (index_b(t_i))

    end; % for t_i=1:nb_tau

  end;  % for m_i

end;  % for l_i

% end of collocation equation section

% steady state point condition x1

xx1=x1;
for i=1:nb_tau
  xx1=[xx1 x1];
end;

res(nml+1:nml+n,1)=sys_rhs(xx1,par);

% derivatives to x(t-tau(i))

for i=0:nb_tau
  J(nml+1:nml+n,nml_n_1+nb_par+1:nml_n_1+nb_par+n)=...
      J(nml+1:nml+n,nml_n_1+nb_par+1:nml_n_1+nb_par+n)...
         +sys_deri(xx1,par,i,[],[]);
end;

% derivatives to eta

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


% steady state point condition x2

% derivatives to x(t-tau(i))

xx2=x2;
for i=1:nb_tau
  xx2=[xx2 x2];
end;

res(nml_n_1:nml+2*n)=sys_rhs(xx2,par);

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

% derivatives to eta