bifdetec

 Check for change of sign(det(A)). If yes, locate bif-point and save. 
 Inputs: 
 p.u, p.lam, p.deta=old point data;   
 u1, lam1, tau1, Gu, Glam at new point 

 For p.bifchecksw=1 check for consistency with EVals of Gu; 
 For p.bifchecksw>1 only det(A) used

0001 function p=bifdetec(p,u1,lam1,tau1,Gu,Glam,ineg0)
0002 % Check for change of sign(det(A)). If yes, locate bif-point and save.
0003 % Inputs:
0004 % p.u, p.lam, p.deta=old point data;
0005 % u1, lam1, tau1, Gu, Glam at new point
0006 %
0007 % For p.bifchecksw=1 check for consistency with EVals of Gu;
0008 % For p.bifchecksw>1 only det(A) used
0009 ds=p.ds; neq=p.neq; np=p.np; xi=0.5; % use balanced xi for det(A)
0010 amat=[[Gu Glam]; [xi*tau1(1:neq*np)' (1-xi)*tau1(neq*np+1)]]; % at new point
0011 % use sign(det(A))=sign(product of neig EVals closest to zero)
0012 % (assumption: negative Evals are all contained in muv)
0013 if p.eigsstart==1; vs=size(amat,1); p.evopts.v0=ones(vs,1)/vs; end 
0014 [V,mu]=eigs(amat,p.neig,p.eigref,p.evopts); muv=mu*ones(1,p.neig)'; 
0015 deta0=p.deta; p.deta=prod(sign(real(muv))); 
0016 if(p.restart>0 || abs(p.deta-deta0)<1.1) % restart, or no bif, do nothing
0017     return; end 
0018 % Thus now bif detected, possibly cross-check with ineg;
0019 ineg=-1; % only occurs for spcalcsw=0
0020 if(p.spcalcsw>0 && p.bifchecksw<2) % catch spurious detA sign change
0021    ineg=p.ineg;
0022    if(ineg0==ineg)
0023       fprintf('Warning: ignored detA_old-detA=%i with ineg_old=ineg=%i.\n',...
0024       deta0-p.deta,ineg);
0025       fprintf('\tSet bifchecksw>1 to force detection.\n'); 
0026       return; end
0027 end 
0028 if(p.vsw>0); fprintf('bif between %g and %g, ds=%g, detA_old=%i detA=%i \n',...
0029            p.lam,lam1, ds, deta0, p.deta);  end 
0030 % Now localize!
0031 u0=p.u; lam0=p.lam; tau0=p.tau; ds=ds/2; bisecc=0; % start bisection
0032 while (bisecc<p.bisecmax && abs(ds)>p.dsmin)     
0033    lamd=tau0(neq*np+1); 
0034    if (p.biflocsw~=1 || abs(lamd)<p.lamdtol/100) % tangent predictor
0035      u1=u0+ds*tau0(1:neq*np); lam1=lam0+ds*lamd; % predictor
0036      if(p.vsw>1) fprintf('checking lam=%g ...',lam1); end 
0037      if(p.parasw==0 || (p.parasw==1 && abs(lamd)>p.lamdtol)) % fixed lam corrector
0038         [u1,res,iter,Gu,Glam]=nlooppde(p,u1,lam1);
0039      else
0040         [u1,lam1,res,iter,Gu,Glam]=nloopext(p,u1,lam1,ds); % arclength corrector
0041      end
0042    end
0043    if p.biflocsw==1 % secant predictor with fixed lam-correction
0044      un=u0+ds*(u1-u0); lamn=lam0+ds*(lam1-lam0); % predictor
0045      if(p.vsw>1) fprintf('checking lam=%g ...',lamn); end 
0046      [un,res,iter,Gu,Glam]=nlooppde(p,un,lamn);
0047    end
0048    if(res<2*p.tol) % If step was (rather) OK, then
0049      amat=[[Gu Glam]; [p.xi*tau0(1:neq*np)' (1-p.xi)*tau0(neq*np+1)]];
0050      tau1=amat\[zeros(neq*np,1);1]; tau1=tau1/xinorm(tau1,p.xi); % tau at new point
0051      amat=[[Gu Glam]; [xi*tau1(1:neq*np)' (1-xi)*tau1(neq*np+1)]]; % amat at new point
0052      if p.eigsstart==1; vs=size(amat,1); p.evopts.v0=ones(vs,1)/vs; end
0053      [V,mu]=eigs(amat,p.neig,p.eigref,p.evopts); 
0054      muv=mu*ones(1,p.neig)'; detab=sign(prod(real(muv))); 
0055      if(p.vsw>1) fprintf('ok, deta=%i\n',detab); end 
0056      if p.biflocsw~=1 % tangent bisection
0057        u0=u1; lam0=lam1; tau0=tau1; % for the next step in bisec
0058        if(p.vsw>1) fprintf('ok, deta=%i\n',detab); end 
0059        if(detab==deta0) ds=ds/2; % still before bif
0060        else ds=-ds/2; deta0=detab; end; % 'midpoint' is behind bif., flip direction
0061      end
0062      if p.biflocsw==1 % secant bisection
0063        if(detab==deta0) u0=un; lam0=lamn; ds=0.5; % still before bif
0064        else u1=u0; lam1=lam0; u0=un; lam0=lamn; deta0=detab; ds=0.5; end; 
0065        % 'midpoint' is behind bif., flip direction
0066      end
0067    else 
0068      if(p.vsw>1) fprintf('no convergence, loc. might be bad ...\n'); end
0069      ds=3*ds/4; % might have hit sing point, try a smaller step in predictor
0070    end % if res<2*p.tol
0071    bisecc=bisecc+1;
0072 end % localization done
0073 
0074 % postprocessing!
0075 us=p.u;lams=p.lam;taus=p.tau; % remember values behind bif
0076 p.u=u0; p.lam=lam0; p.tau=tau0; % and store current values for save
0077 brout=p.outfu(p,u0,lam0); % append to bifvals
0078 p.bifvals=[p.bifvals [p.count; p.bcount; ineg; p.lam; p.err; brout]]; 
0079 if(p.errchecksw>0); p.err=errcheck(p); end 
0080 p.branch=[p.branch [p.count; ineg; p.lam; p.err; brout]]; % put on branch
0081 fname=[p.bpname,sprintf('%i',p.bcount),'.mat'];save(fname,'p');    
0082 fprintf('%4i  %7.5e (BP, saved to %s)\n',p.count,lam0,fname);
0083 if(p.biflocsw~=1) figure(p.brfig);plot(lam1,brout(p.bpcmp),'o'); end
0084 if(p.biflocsw==1) figure(p.brfig);plot(lam0,brout(p.bpcmp),'o'); end
0085 if((p.count>0) && (mod(p.count,p.smod)==0)) % save to file
0086    fname=[p.pname,sprintf('%i',p.count),'.mat']; save(fname,'p'); 
0087    fprintf('also saved to %s\n',fname);
0088 end 
0089 p.u=us;p.lam=lams; p.tau=taus; % restore values behind bif
0090 p.bcount=p.bcount+1; p.count=p.count+1; 
0091 
0092