neimarksackermapr.m

计算动力学系统的分岔图

function out = neimarksackermap
%
% Neimark_sacker curve definition file for a problem in mapfile
% 

global cds nsmds
    out{1}  = @curve_func; 
    out{2}  = @defaultprocessor;
    out{3}  = @options; 
    out{4}  = @jacobian;
    out{5}  = @hessians;
    out{6}  = @testf;
    out{7}  = @userf;
    out{8}  = @process;
    out{9}  = @singmat;
    out{10} = @locate;
    out{11} = @init;
    out{12} = @done;
    out{13} = @adapt;
    
return
%-------------------------------------------------------
function func = curve_func(arg)
global cds nsmds
  [x,p,k] = rearr(arg); p=n2c(p);n=nsmds.Niterations;
  jac = nsmjac(x,p,n);
  RED=jac*jac-2*k*jac+eye(nsmds.nphase);
  Bord=[RED nsmds.borders.w;nsmds.borders.v' zeros(2)];
  bunit=[zeros(nsmds.nphase,2);eye(2)];
  vext=Bord\bunit;
  vext1=vext(nsmds.nphase+nsmds.index1(1),nsmds.index1(2));
  vext2=vext(nsmds.nphase+nsmds.index2(1),nsmds.index2(2));
  x1=x;
  for i=1:n
  x1=feval(nsmds.func,0,x1,p{:});
  end
  func = [x1-x; vext1;vext2];
   
%--------------------------------------------------------  
function jac = jacobian(varargin)
 global cds nsmds


    n=nsmds.Niterations;
    nap = length(nsmds.ActiveParams); 
    x0 = varargin{1}; [x,p,k] = rearr(x0); p = n2c(p);     
    J=nsmjac(x,p,n);
    RED=J*J-2*k*J+eye(nsmds.nphase);
    Bord=[RED nsmds.borders.w;nsmds.borders.v' zeros(2)];
    bunit=[zeros(nsmds.nphase,2);eye(2)];
    vext=Bord\bunit;
    wext=Bord'\bunit;
    jac=[J-eye(nsmds.nphase)  nsmjacp(x,p,n) zeros(nsmds.nphase,1)];
    hess=nsmhess(x,p,n);
     for i=1:nsmds.nphase
        jac(nsmds.nphase+1,i)=-wext(1:nsmds.nphase,nsmds.index1(1))'*(J*hess(:,:,i)+hess(:,:,i)*J-2*k*hess(:,:,i))*vext(1:nsmds.nphase,nsmds.index1(2));
        jac(nsmds.nphase+2,i)=-wext(1:nsmds.nphase,nsmds.index2(1))'*(J*hess(:,:,i)+hess(:,:,i)*J-2*k*hess(:,:,i))*vext(1:nsmds.nphase,nsmds.index2(2));
    end
    jac(nsmds.nphase+1,nsmds.nphase+nap+1)=2*wext(1:nsmds.nphase,nsmds.index1(1))'*J*vext(1:nsmds.nphase,nsmds.index1(2));
    jac(nsmds.nphase+2,nsmds.nphase+nap+1)=2*wext(1:nsmds.nphase,nsmds.index2(1))'*J*vext(1:nsmds.nphase,nsmds.index2(2));
    hessp=nsmhessp(x,p,n);
    for i=1:nap
        jac(nsmds.nphase+1,nsmds.nphase+i)=-wext(1:nsmds.nphase,nsmds.index1(1))'*(J*hessp(:,:,i)+hessp(:,:,i)*J-2*k*hessp(:,:,i))*vext(1:nsmds.nphase,nsmds.index1(2));
        jac(nsmds.nphase+2,nsmds.nphase+i)=-wext(1:nsmds.nphase,nsmds.index2(1))'*(J*hessp(:,:,i)+hessp(:,:,i)*J-2*k*hessp(:,:,i))*vext(1:nsmds.nphase,nsmds.index2(2));
    end 
%x,jac,
%------------------------------------------------------
function hess = hessians(varargin)  
hess =[];
%------------------------------------------------------
function varargout = defaultprocessor(varargin)
global cds nsmds
%elseif strcmp(arg,'defaultprocessor')
  if nargin > 2
    s = varargin{3};
    varargout{3} = s;
  end
 % compute eigenvalues?
  if (cds.options.Multipliers==1)
%       [tspan,y0,odeopt] = feval(hds.mapfile, [], [], 'init', []);
      n=nsmds.Niterations;
      x0 = varargin{1}; [x,p] = rearr(x0); p = n2c(p);
      jac =nsmjac(x,p,n);
      varargout{2} = eig(jac);
  else
      varargout{2}= nan;
  end  

  % all done succesfully
  varargout{1} = 0;
%-------------------------------------------------------  
function option = options
global cds nsmds
  option = contset;
  % Check for symbolic derivatives in mapfile
  
  symjac  = ~isempty(nsmds.Jacobian);
  symhes  = ~isempty(nsmds.Hessians);
  symDer3 = ~isempty(nsmds.Der3);
  symDer4 = ~isempty(nsmds.Der4);
  symDer5 = ~isempty(nsmds.Der5);
    
  symord = 0; 
  if symjac, symord =  1; end
  if symhes, symord =  2; end
  if symDer3, symord = 3; end
  if symDer4, symord = 4; end
  if symDer5, symord = 5; end

  option = contset(option, 'SymDerivative', symord);
  option = contset(option, 'Workspace', 1);
  option = contset(option, 'Locators', [0 0 0]);

  symjacp = ~isempty(nsmds.JacobianP); 
  symhessp = ~isempty(nsmds.HessiansP); 
  symordp = 0;
  if symjacp,  symordp = 1; end
  if symhessp, symordp = 2;end
  option=contset(option,'SymDerivativeP',symordp);
  
  cds.symjac  = 1;%1
  cds.symhess = 1;

% ---------------------------------------------------------------

function [out, failed] = testf(id, x, v)
global  cds nsmds 
n=nsmds.Niterations;nphase = nsmds.nphase;
[x0,p,k] = rearr(x); p1 = n2c(p);
jac=nsmjac(x0,p1,n);%x,eig(jac),pause
RED=jac*jac-2*k*jac+eye(nphase);
Bord=[RED nsmds.borders.w;nsmds.borders.v' zeros(2)];
bunit=[zeros(nphase,2);eye(2)];
vext=Bord\bunit;
wext=Bord'\bunit;
failed = [];
for i=id
  lastwarn('');
switch i
  case 1 % CH (Chenciner) 
     [V,D] = eig(jac);
     % find pair of complex eigenvalues
     d = diag(D);idx1=0;idx2=0;
     for i=1:nphase
        for j=i+1:nphase
            if (d(i)== conj(d(j))) & (imag(d(i))~=0)&(abs(real(d(i)-k))<0.000001)&(norm(d(i)*d(j)-1)<0.000001)
                idx1=i;idx2=j;norm(d(i)*d(j)-1);
            end
        end
     end
    if (idx1==0)||(imag(d(idx1))==0);
       out(1)=111;
         
    else
       %temp=idx1;
       if imag(d(idx1))<0
       idx1=idx2;
       end  
       q=V(:,idx1);
    end
  if (idx1==0)||(imag(d(idx1))==0);
       out(1)=111;         
  else
   [V,D] = eig(jac');  
   d1  = diag(D);
   [Y,i]=min(abs(d1-conj(d(idx1))));    
   p=V(:,i);
   q=q/norm(q);
   p=p/(q'*p);
   out(1)=nf_NSm(nsmds.func,nsmds.Jacobian,nsmds.Hessians,nsmds.Der3,jac,q,p,nphase,x0,p1,d(idx1),n);
 end
  
  case 2 % PDNS  (Flip+Neimark-Sacker)
      out(2) =det(jac+eye(nsmds.nphase));
  case 3 % LPNS  (Fold+Neimark-Sacker) 
      out(3)=det(jac-eye(nsmds.nphase));
  case 4   %R1 (Resonance 1:1)
      out(4)= k-1;
  case 5 % NSNS  (Double Neimark-Sacker) 
      
     [Q,R]=qr(wext(1:nsmds.nphase,1:2));
      Q1=Q(1:nsmds.nphase,3:nsmds.nphase); 
      A = [Q1'*jac*Q1 zeros(nsmds.nphase-2,1)];
      [bialt_M1,bialt_M2,bialt_M3,bialt_M4]=bialtaa(nsmds.nphase-2); 
      B=A(bialt_M1).*A(bialt_M2)-A(bialt_M3).*A(bialt_M4);
      B=B-eye(size(B,1));
      out(5) = det(B);
      
  case 6  % R2  (Resonance1:2)
      out(6)=k+1;
  case 7  % R3  (Resonance1:3)
      out(7)=k+1/2.0;
  case 8  % R4  (Resonance1:4)
      out(8)=k;
  otherwise
    error('No such testfunction');
  end  
end

%-------------------------------------------------------
function [out, failed] = userf(userinf, id, x, v)
global cds nsmds
dim =size(id,2);
failed = [];
for i=1:dim
  lastwarn('');
  [x0,p] = rearr(x); p = n2c(p);
  if (userinf(i).state==1)
      out(i)=feval(nsmds.user{id(i)},0,x0,p{:});
      else
      out(i)=0;
  end
  if ~isempty(lastwarn)
    msg = sprintf('Could not evaluate userfunction %s\n', id(i).name);
    failed = [failed i];
  end
end
% ---------------------------------------------------------------

function [failed,s] = process(id, x, v, s)
global cds nsmds opt
ndim = cds.ndim;
n=nsmds.Niterations;
[x0,p,k] = rearr(x); p1 = n2c(p);
nphase=size(x0,1);
% WM: Removed SL array
fprintf('label = %s, x = ', s.label); printv(x);
jac=nsmjac(x0,p1,n);
RED=jac*jac-2*k*jac+eye(nsmds.nphase);
Bord=[RED nsmds.borders.w;nsmds.borders.v' zeros(2)];
bunit=[zeros(nsmds.nphase,2);eye(2)];
vext=Bord\bunit;
wext=Bord'\bunit;
d1=k+sqrt(-1.0)*abs(sqrt(1-k*k));
d2=conj(d1);
alpha=vext(1:nphase,1)'*(jac*vext(1:nphase,2)-d1*vext(1:nphase,2));
beta=-vext(1:nphase,1)'*(jac*vext(1:nphase,1)-d1*vext(1:nphase,1));
q=alpha*vext(1:nphase,1)+beta*vext(1:nphase,2);
alpha=wext(1:nphase,1)'*(jac'*wext(1:nphase,2)-d2*wext(1:nphase,2));
beta=-wext(1:nphase,1)'*(jac'*wext(1:nphase,1)-d2*wext(1:nphase,1));
p=alpha*wext(1:nphase,1)+beta*wext(1:nphase,2);
q=q/norm(q);
p=p/(q'*p);
switch id
  case 1 % CH
    s.data.c=nf_CHm(nsmds.func,nsmds.Jacobian,nsmds.Hessians,nsmds.Der3,nsmds.Der4,nsmds.Der5,jac,q,p,nphase,x0,p1,n);
    fprintf('Normal form coefficient of CH = %d\n', s.data.c); 
    s.msg  = sprintf('Chenciner bifurcation');  
  case 2 % PDNS
    [V,D]=eig(jac+eye(nphase));
    [Y,i]=min(abs(diag(D)));
    q2=real(V(:,i));
    [V,D]=eig(jac'+eye(nphase));
    [Y,i]=min(abs(diag(D)));
    p2=real(V(:,i));
    q2=q2/norm(q2);
    p2=p2/(q2'*p2);
    s.data.c=nf_PDNSm(nsmds.func,nsmds.Jacobian,nsmds.Hessians,nsmds.Der3,jac,q2,p2,q,p,nphase,x0,p1,n);
    fprintf('Normal form coefficient for PDNS :[a , b, c, d]= %d, %d, %d, %d\n', s.data.c);
    s.msg  = sprintf('Flip+Neimark-Sacker');
  case 3 % LPNS
      
    [V,D]=eig(jac-eye(nphase));
    [Y,i]=min(abs(diag(D)));
    q2=real(V(:,i));
    [V,D]=eig(jac'-eye(nphase));
    [Y,i]=min(abs(diag(D)));
    p2=real(V(:,i));
    q2=q2/norm(q2);
    p2=p2/(q2'*p2);