init_gpd_lp2m.m

计算动力学系统的分岔图

function [x0,v0]= init_DPDm_LP2(mapfile,eps,x,p,ap,n,varargin)
% 
% [x0,v0]= init_DPDm_LP2(mapfile,eps,x,p,ap,n)
% Initializes a fold continuation from a degenerate flip point
%
global cds lpmds
% check input
if size(ap,2)~=2
  errordlg('Two active parameter are needed for a Limitpoint bifurcation curve continuation');
end
% initialize lpmds
lpmds.mapfile = mapfile;
func_handles = feval(lpmds.mapfile);
lpmds.func = func_handles{2};
lpmds.Jacobian  = func_handles{3};
lpmds.JacobianP = func_handles{4};
lpmds.Hessians  = func_handles{5};
lpmds.HessiansP = func_handles{6};
lpmds.Der3      = func_handles{7};
lpmds.Der4      = func_handles{8};
lpmds.Der5      = func_handles{9};
lpmds.Niterations=2*n;
siz = size(func_handles,2);
if siz > 9
  j=1;
  for i=10:siz
    lpmds.user{j}= func_handles{i};
    j=j+1;
  end
end
lpmds.nphase = size(x,1);
lpmds.ActiveParams = ap;
lpmds.P0 = p;
if size(varargin,1)>0,lpmds.BranchParams=varargin{1};else lpmds.BranchParams=[];end 
cds.curve = @limitpointmap;
cds.ndim = length(x)+2;
%-----Defining Symbolic derivatives-----
  symjac  = ~isempty(lpmds.Jacobian);
  symhes  = ~isempty(lpmds.Hessians);
  symDer3 = ~isempty(lpmds.Der3);
  symDer4 = ~isempty(lpmds.Der4);
  symDer5 = ~isempty(lpmds.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
  cds.options.SymDerivative = symord;
  symjacp  = ~isempty(lpmds.JacobianP); 
  symhessp = ~isempty(lpmds.HessiansP); 
  symordp = 0;
  if symjacp,  symordp = 1; end
  if symhessp, symordp = 2; end
  cds.options.SymDerivativeP = symordp;
%----------------------
p0 = p;p = n2c(p);
global T1global T2global T3global T4global T5global
if (cds.options.SymDerivative >=5)
  T1global=tens1(lpmds.func,lpmds.Jacobian,x,p,n);
  T2global=tens2(lpmds.func,lpmds.Hessians,x,p,n);
  T3global=tens3(lpmds.func,lpmds.Der3,x,p,n);
  T4global=tens4(lpmds.func,lpmds.Der4,x,p,n);
  T5global=tens5(lpmds.func,lpmds.Der5,x,p,n);   
end
  hessIncrement =(cds.options.Increment)^(3.0/4.0);
  ten3Increment =(cds.options.Increment)^(3.0/5.0);
  ten4Increment =(cds.options.Increment)^(3.0/6.0);
  ten5Increment =(cds.options.Increment)^(3.0/7.0);
%---Branch Switching Algorithm----
  nphase = size(x,1);
  A = lpmjac(x,p,n);
  [X,D] = eig(A+eye(nphase));
  [Y,i] = min(abs(diag(D)));
  vext = X(:,i)/norm(X(:,i));
  [X,D] = eig(A'+eye(nphase));
  [Y,i] = min(abs(diag(D)));
  wext = X(:,i)/(X(:,i)'*vext);
  c2 = nf_DPDm(lpmds.func,lpmds.Jacobian,lpmds.Hessians,lpmds.Der3,lpmds.Der4,lpmds.Der5,A,vext,wext,nphase,x,p,n);
  A1 = lpmjacp(x,p,n);   							%jacobianp
  temp = (eye(nphase)-A)\A1;							%temp=(I-A)^{INV}*J1
  s1=[1;0];s2=[0;1];	
  AA=zeros(nphase,nphase,n);
  wwt(:,:,n)=eye(nphase);
  x1=x;
  xit=zeros(nphase,n);xit(:,1)=x1;
  AA(:,:,1)=lpmjac(x1,p,1);
   for m=2:n
    x1=feval(lpmds.func,0,x1,p{:});
    xit(:,m)=x1;
    AA(:,:,m)=lpmjac(x1,p,1);
   end%define standard vectors  
  
  test1 = lphesspvect(xit,p,vext,AA,n)*s1; 						% A1(q,s1) 
  test1 = test1 +  multilinear2(lpmds.func,vext,temp*s1,x,p,n,hessIncrement); 	% +B(q,temp*s1)
  gamma1= wext'*test1;
  test2 = lphesspvect(xit,p,vext,AA,n)*s2; 						% A1(q,s2) 
  test2 = test2 +  multilinear2(lpmds.func,vext,temp*s2,x,p,n,hessIncrement); 	% +B(q,temp*s2)
  gamma2= wext'*test2;
  s1 = -[gamma1;gamma2]/(gamma1^2 + gamma2^2);					% new orthogonal basis
  s2 = [-gamma2;gamma1];
  h200 = (eye(nphase)-A)\(multilinear2(lpmds.func,vext,vext,x,p,n,hessIncrement));
  Abor = [A+eye(nphase) vext ; wext' 0];
  hh110 = Abor\[(gamma1*test1 + gamma2*test2)/(gamma1^2 + gamma2^2)-vext ; 0];
  hh101 = Abor\[(gamma2*test1 - gamma1*test2) ; 0];
  hh110 = hh110(1:nphase);hh101 = hh101(1:nphase); 		
%h110 = hh110 +delta1*hh101; h101 = delta2*hh101;
%
%Computation of B_1 and C_1 vectors, symbolic derivatives have to be redefined
%since they are fixed with parameters. When numerically differentiating we scale
%s1 and s2 to norm 1.
%temp1 = B_1(q,q,s1); temp3 = C_1(q,q,q,s1) + 3 B_1(q,h200,s1);
%temp2 = B_1(q,q,s2); temp4 = C_1(q,q,q,s2) + 3 B_1(q,h200,s2);
%---------------------------------------------
%wrt to s1
  p1 =p0;p1(ap) = p1(ap) + cds.options.Increment*s1/norm(s1);p1=n2c(p1);
  if (cds.options.SymDerivative >=3)
    T1global=tens1(lpmds.func,lpmds.Jacobian,x,p1,n);
    T2global=tens2(lpmds.func,lpmds.Hessians,x,p1,n);
    T3global=tens3(lpmds.func,lpmds.Der3,x,p1,n);
  end
  temp1 = multilinear2(lpmds.func,vext,vext,x,p1,n,hessIncrement);
  temp3 = multilinear3(lpmds.func,vext,vext,vext,x,p1,n,ten3Increment);
  temp3 = temp3 + 3.0*multilinear2(lpmds.func,vext,h200,x,p1,n,hessIncrement);
  p1 =p0;p1(ap) = p1(ap) - cds.options.Increment*s1/norm(s1);p1=n2c(p1);
  if (cds.options.SymDerivative >=3)
    T1global=tens1(lpmds.func,lpmds.Jacobian,x,p1,n);
    T2global=tens2(lpmds.func,lpmds.Hessians,x,p1,n);
    T3global=tens3(lpmds.func,lpmds.Der3,x,p1,n);
  end
  temp1 = temp1 - multilinear2(lpmds.func,vext,vext,x,p1,n,hessIncrement);
  temp3 = temp3 - multilinear3(lpmds.func,vext,vext,vext,x,p1,n,ten3Increment);
  temp3 = temp3 - 3.0*multilinear2(lpmds.func,vext,h200,x,p1,n,hessIncrement);
  temp1 = temp1*norm(s1)/(2.0*cds.options.Increment);
  temp3 = temp3*norm(s1)/(2.0*cds.options.Increment);
%wrt to s2
  p1 =p0;p1(ap) = p1(ap) + cds.options.Increment*s2/norm(s2);p1=n2c(p1);
  if (cds.options.SymDerivative >=3)
    T1global=tens1(lpmds.func,lpmds.Jacobian,x,p1,n);
    T2global=tens2(lpmds.func,lpmds.Hessians,x,p1,n);
    T3global=tens3(lpmds.func,lpmds.Der3,x,p1,n);
  end
  temp2 = multilinear2(lpmds.func,vext,vext,x,p1,n,hessIncrement);
  temp4 = multilinear3(lpmds.func,vext,vext,vext,x,p1,n,ten3Increment);
  temp4 = temp4 + 3.0*multilinear2(lpmds.func,vext,h200,x,p1,n,hessIncrement);
  p1 =p0;p1(ap) = p1(ap) - cds.options.Increment*s2/norm(s2);p1=n2c(p1);
  if (cds.options.SymDerivative >=3)
    T1global=tens1(lpmds.func,lpmds.Jacobian,x,p1,n);
    T2global=tens2(lpmds.func,lpmds.Hessians,x,p1,n);
    T3global=tens3(lpmds.func,lpmds.Der3,x,p1,n);
  end
  temp2 = temp2 - multilinear2(lpmds.func,vext,vext,x,p1,n,hessIncrement);
  temp4 = temp4 - multilinear3(lpmds.func,vext,vext,vext,x,p1,n,ten3Increment);
  temp4 = temp4 - 3.0*multilinear2(lpmds.func,vext,h200,x,p1,n,hessIncrement);
  temp2 = temp2*norm(s2)/(2.0*cds.options.Increment);
  temp4 = temp4*norm(s2)/(2.0*cds.options.Increment);;

%wrt to original parameter
  if (cds.options.SymDerivative >=3)
    T1global=tens1(lpmds.func,lpmds.Jacobian,x,p1,n);
    T2global=tens2(lpmds.func,lpmds.Hessians,x,p1,n);
    T3global=tens3(lpmds.func,lpmds.Der3,x,p1,n);
  end
%--------------------------------------------
%Continue
  test1 = lphesspvect(xit,p,h200,AA,n)*s1; 							%  A1(h200,s1)
  test1 = test1 + 2*multilinear2(lpmds.func,vext,hh110,x,p,n,hessIncrement);		%+2B(q,hh110)
  test1 = test1 +   multilinear2(lpmds.func,h200,temp*s1,x,p,n,hessIncrement);		%+ B(h200,h010)
  test1 = test1 +   multilinear3(lpmds.func,vext,vext,temp*s1,x,p,n,ten3Increment);	%+ C(q,q,h010)
  test1 = test1 +   temp1;								% see above for temp1
  hh210 = (A-eye(nphase))\(2*h200-test1);
  
  test2 = lphesspvect(xit,p,h200,AA,n)*s2; 							%  A1(h200,s2)
  test2 = test2 + 2*multilinear2(lpmds.func,vext,hh101,x,p,n,hessIncrement);		%+2B(q,hh101)
  test2 = test2 +   multilinear2(lpmds.func,h200,temp*s2,x,p,n,hessIncrement);		%+ B(h200,h001)
  test2 = test2 +   multilinear3(lpmds.func,vext,temp*s2,vext,x,p,n,ten3Increment);	%+ C(q,q,h001)
  test2 = test2 +   temp2;								% see above for temp2
  hh201 = (eye(nphase)-A)\test2;
%h210 = (A-I)\(2h200-test1-delta1*test2) = hh210 + delta1*hh201 
%h201 = (A-I)\(-delta2*test2) = delta2*hh201
  RHS3 = multilinear3(lpmds.func,vext,vext,vext,x,p,n,ten3Increment);			%   C(q,q,q)
  RHS3 = RHS3 + 3.0*multilinear2(lpmds.func,vext,h200,x,p,n,hessIncrement);   		% +3B(q,h200)
  a = wext'*RHS3/6.0;
  h300 =  Abor\[6.0*a*vext - RHS3; 0];h300 = h300(1:nphase);
%----
  test3 = lphesspvect(xit,p,h300,AA,n)*s1; 							%  A1(h300,s1)
  test3 = test3 + multilinear4(lpmds.func,vext,vext,vext,temp*s1,x,p,n,ten4Increment);	%+ D(q,q,q,h010)
  test3 = test3 + 3.0*multilinear3(lpmds.func,vext,vext,hh110,x,p,n,ten3Increment);	%+3C(q,q,hh110)