p_correc.m

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

function [point,success]=p_correc(point,free_par,step_cnd,method,p_nr,...
                                  previous,d_nr,tz)

% function [point,success]=p_correc(point0,free_par,step_cnd,method,adapt,
%                                 previous)
% INPUT:
%   point0 initial point guess
%   free_par free parameter numbers in N^d 
%       step_cnd steplength condition(s) as point(s)
%   method method parameters
%   adapt if zero or absent, do not adapt mesh; if one, always adapt 
%       previous (optional) previously computed branch point (used, in case 
%            of periodic solutions or connecting orbits, to
%            minimize phase shift)
% OUTPUT:
%       point final solution
%       success nonzero for successfull computation
%
% OPTIONAL EXTRA INPUT:
%       d_nr number of delay crossing zero
%       tz (for psol only) time point for which tau(tz)=0 and dtau(tz)/dt=0

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

% optional parameters:

if ~exist('p_nr'),
  p_nr=0;
end;
if ~exist('previous'),
  previous=[];
end;

if ~exist('d_nr'),
  d_nr=0;
else
  d_tau=sys_ntau;
end;

% some method parameters:

max_iter=method.newton_max_iterations; % max number of newton iterations
nmon_iter=method.newton_nmon_iterations; % max number of nonmonotone iterations
conv_r=method.halting_accuracy; % required accuracy for convergence
print_r=method.print_residual_info; % print residual evolution

% remove old stability info if present:

if isfield(point,'stability')
  point.stability=[];
end;

% initialize:

switch point.kind,
  case 'stst',
    n=size(point.x,1);
    p_start=n;
  case 'fold',
    n=size(point.x,1);
    p_start=2*n;
    c=point.v'/(point.v'*point.v);
  case 'hopf',
    n=size(point.x,1);
    p_start=3*n+1;
    vr=real(point.v);
    vi=imag(point.v);
    vn=vr'*vr+vi'*vi;
    c=point.v'/vn; % complex conjugate transpose
  case {'psol','hcli'},
    n=size(point.profile,1);
    col=method.collocation_parameters;
    mm=point.degree;
    if ~isempty(col) & length(col)~=mm
      err=[point.degree mm]
      error('P_CORREC: number of collocation parameters differs from degree!');
    end;
    ll=floor((size(point.profile,2)-1)/mm);
    if ll~=(size(point.profile,2)-1)/mm,
      error('P_CORREC: point mesh does not contain l intervals of m points!');
    end;
    p_start=n*mm*ll+1+n;
    if d_nr==0
      ph=method.phase_condition;
    else
      ph=0;
    end;
    if isempty(point.mesh)
      mesh=0:1/(ll*mm):1;
    else
      mesh=point.mesh;
      ma=method.adapt_mesh_before_correct;
      if p_nr>1 & mod(p_nr,ma)==0 
        % do not adapt mesh when p_nr=0
        % do not (yet) adapt mesh when p_nr=1
        % adapt mesh when p_nr>1 & mod(p_nr,ma)=0
        new_mesh=psol_msh(mesh,mm,point.profile,ll,point.degree);
        switch point.kind,
          case 'psol',
            point.profile=psol_eva(point.profile,mesh,new_mesh,mm);
          case 'hcli',
            point.profile=hcli_eva(point.profile,mesh,new_mesh,mm);
        end;
        point.mesh=new_mesh;
        mesh=new_mesh;
      end;
    end;
  otherwise,
    err=point.kind
    error('P_CORREC: point kind not recognized.');
  end;

% store parameters:

par=point.parameter;

% perform Newton-Raphson iterations:

for i=1:max_iter

  % compute & apply corrections:

  switch point.kind  

    case 'stst',
      % produce jacobian
      [J,res]=stst_jac(point.x,par,free_par);
      % add (linear) steplength conditions
      for j=1:length(step_cnd)
        J(n+j,1:n)=step_cnd(j).x';
        for k=1:length(free_par)
          J(n+j,n+k)=step_cnd(j).parameter(free_par(k));
        end;
        res(n+j,1)=0;
      end;
      if d_nr~=0
        % add condition on zero delay
        kk=size(J,1);
        for jj=0:d_tau
          xx(:,jj+1)=point.x;
        end;
        J(kk+1,1:n)=0;
        for j=0:d_tau
          J(kk+1,1:n)=J(kk+1,1:n)+sys_dtau(d_nr,xx,par,j,[]);
        end;
        for k=1:length(free_par)
          J(kk+1,n+k)=sys_dtau(d_nr,xx,par,[],free_par(k));
        end;
        res(kk+1,1)=sys_tau(d_nr,xx,par);
      end;
      % add extra conditions
      if method.extra_condition
        [resi,condi]=sys_cond(point);
        for j=1:length(condi)
          k=size(J,1);
          res(k+1)=resi(j);
          J(k+1,1:n)=condi(j).x';
          for o=1:length(free_par)
            J(k+1,n+o)=condi(j).parameter(free_par(o));
          end;
        end;
      end;
      % solve linear system
      if size(J,1)~=size(J,2)
        disp('P_CORREC warning: use of nonsquare Jacobian.');
      end;

      dx=J\res;
      % apply non-parameter corrections
      point.x=point.x-dx(1:n);

    case 'fold',
      % produce jacobian
      [J,res]=fold_jac(point.x,point.v,par,free_par,c);
      % add (linear) steplength conditions
      for j=1:length(step_cnd)
        J(2*n+1+j,1:n)=step_cnd(j).x';
        J(2*n+1+j,n+1:2*n)=step_cnd(j).v';
        for k=1:length(free_par)
          J(2*n+1+j,2*n+k)=step_cnd(j).parameter(free_par(k));
        end;
        res(2*n+1+j,1)=0;
      end;
      if d_nr~=0
        % add condition on zero delay
        kk=size(J,1);
        for jj=0:d_tau
          xx(:,jj+1)=point.x;
        end;
        J(kk+1,1:2*n+2)=0;
        for j=0:d_tau
          J(kk+1,1:n)=J(kk+1,1:n)+sys_dtau(d_nr,xx,par,j,[]);
        end;
        for k=1:length(free_par)
          J(kk+1,2*n+k)=sys_dtau(d_nr,xx,par,[],free_par(k));
        end;
        res(kk+1,1)=sys_tau(d_nr,xx,par);
      end;
      % add extra conditions
      if method.extra_condition
        [resi,condi]=sys_cond(point);
        for j=1:length(condi)
          k=size(J,1);
          res(k+1)=resi(j);
          J(k+1,1:n)=condi(j).x';
          J(k+1,n+1:2*n)=condi(j).v';
          for o=1:length(free_par)
            J(k+1,2*n+o)=condi(j).parameter(free_par(o));
          end;
        end; 
      end;
      % solve linear system
      if size(J,1)~=size(J,2)
        disp('P_CORREC warning: use of nonsquare Jacobian.');
      end;
      dx=J\res;
      % apply non-parameter corrections
      point.x=point.x-dx(1:n);
      point.v=point.v-dx(n+1:2*n);

    case 'hopf'  
      % produce jacobian
      [J,res]=hopf_jac(point.x,point.omega,point.v,par,free_par,c);
        % add (linear) steplength conditions
        for j=1:length(step_cnd)
          J(3*n+2+j,1:n)=step_cnd(j).x';
          J(3*n+2+j,n+1:2*n)=real(step_cnd(j).v)';
          J(3*n+2+j,2*n+1:3*n)=imag(step_cnd(j).v)';
          J(3*n+2+j,3*n+1)=step_cnd(j).omega;
          for k=1:length(free_par)
            J(3*n+2+j,3*n+1+k)=step_cnd(j).parameter(free_par(k));
          end;
          res(3*n+2+j)=0;
        end;
      if d_nr~=0
        % add condition on zero delay
        kk=size(J,1);
        for jj=0:d_tau
          xx(:,jj+1)=point.x;
        end;
        J(kk+1,1:3*n+3)=0;
        for j=0:d_tau
          J(kk+1,1:n)=J(kk+1,1:n)+sys_dtau(d_nr,xx,par,j,[]);
        end;
        for k=1:length(free_par)
          J(kk+1,3*n+1+k)=sys_dtau(d_nr,xx,par,[],free_par(k));
        end;
        res(kk+1,1)=sys_tau(d_nr,xx,par);
      end;
      % add extra conditions
      if method.extra_condition
        [resi,condi]=sys_cond(point);
        for j=1:length(condi)
          k=size(J,1);
          res(k+1)=resi(j);
          J(k+1,1:n)=condi(j).x';
          J(k+1,n+1:2*n)=real(condi(j).v');
          J(k+1,2*n+1:3*n)=imag(condi(j).v');
          J(k+1,3*n+1)=condi(j).omega;
          for o=1:length(free_par)
            J(k+1,3*n+1+o)=condi(j).parameter(free_par(o));
          end;
        end;
      end;
      % solve linear system
      if size(J,1)~=size(J,2)
        disp('P_CORREC warning: use of nonsquare Jacobian.');
      end;
      dx=J\res;
      % apply non-parameter corrections
      point.x=point.x-dx(1:n);
      point.v=point.v-dx(n+1:2*n)-sqrt(-1)*dx(2*n+1:3*n);
      point.omega=point.omega-dx(3*n+1);

    case 'psol',
      % produce jacobian
      [J,res]=psol_jac(col,point.period,point.profile,mesh,mm,par,...
                         free_par,ph);
      % add (linear) steplength conditions
      for j=1:length(step_cnd)
        for k=1:ll*mm+1
          J(n*(mm*ll+1)+ph+j,(k-1)*n+1:k*n)=step_cnd(j).profile(:,k)';
        end;
        J(n*(mm*ll+1)+ph+j,n*(mm*ll+1)+1)=step_cnd(j).period;
        for k=1:length(free_par)
          J(n*(mm*ll+1)+ph+j,n*(mm*ll+1)+1+k)=...
                                  step_cnd(j).parameter(free_par(k));
        end;
        res(n*(mm*ll+1)+ph+j)=0;
      end;
      if d_nr~=0
        %%%% add 2 conditions: on zero delay (1) and 
        %%%% on zero time derivative of delay (2)

        T=point.period;

        % compute x_tau, dx_tau, d2x_tau 

        for tau_i=1:d_tau+1
          if tau_i==1
            tau(tau_i)=0;
          else