swibra

 load data from file pre/fname and attempt branch switch, 
 new pre and stepsize ds are auxiliary arguments

0001 function q=swibra(pre,fname,varargin)
0002 % load data from file pre/fname and attempt branch switch,
0003 % new pre and stepsize ds are auxiliary arguments
0004 ffname=[pre '/' fname]; s=load(ffname,'p'); p=s.p; u=p.u; 
0005 lam=p.lam; tau=p.tau; xi=p.xi; 
0006 del=1e-3; % don't choose del too small!
0007 if (length(varargin)>1) p.ds=varargin{2}; gotds=1; 
0008 else p.ds=p.dsmax/10; gotds=0; 
0009 end 
0010 r=resi(p,u,lam); [Gu,Glam]=getder(p,u,lam,r); % residual and jacobians
0011 % now calc al0, al1, phi0
0012 u0d=tau(1:p.neq*p.np); al0=tau(p.neq*p.np+1); 
0013 j=1; % use j=1 eigenval nearest to zero
0014 if p.eigsstart==1; vs=size(Gu,1); p.evopts.v0=ones(vs,1)/vs; end % set eigs startvector
0015 [phi1v,mu]=eigs(Gu,j,0,p.evopts);
0016 phi1=phi1v(:,j);
0017 if al0==0 fprintf('Degenerate branch point with d/ds lambda=0');return; end
0018 fprintf('zero eigenvalue is %g\n',mu(j,j));
0019 [psi1v,mu]=eigs(Gu',j,0,p.evopts); psi1=psi1v(:,j);psi1=psi1/(phi1'*psi1); 
0020 al1=psi1'*u0d;  phi0=(u0d-al1*phi1)/al0; 
0021 rp=resi(p,u+del*phi1,lam); % residual and jacobian at u+del*phi1
0022 [Gup,Glamp]=getder(p,u+del*phi1,lam,rp); Gud=Gup-Gu; 
0023 if any(Gud) 
0024    Glamd=Glamp-Glam; a1=psi1'*(Gud*phi1)/del; 
0025    b1=psi1'*(Gud*phi0+Glamd)/del; al1b=-(a1*al1/al0+2*b1); 
0026    fprintf('al0=%g, a1=%g, b1=%g, al1b=%g\n',al0,a1,b1,al1b); 
0027    if al1b==0 fprintf('No distinct branch to switch to.');return; end
0028    tau1=[al1b*phi1+a1*phi0; a1];  
0029 else % trivial branch
0030     tau1=[phi1; 0];fprintf('trivial swibra\n');
0031 end
0032 tau1=tau1/xinorm(tau1,xi);
0033 figure(p.ifig); clf; n0=(p.pcmp-1)*p.np+1; n1=(p.pcmp)*p.np; 
0034 pdeplot(p.points,p.edges,p.tria,'xydata',tau1(n0:n1)); axis tight; box on; 
0035 title(['\tau_' mat2str(p.pcmp) ' at ' fname]); colormap(p.cm); 
0036 xi=1/p.np; % the "natural choice"
0037 if(p.isw>1); p.xi=asknu('xi',xi); 
0038     if ~gotds; p.ds=asknu('ds',p.ds); end 
0039 end
0040 p.tau=tau1; p.bcount=1; p.count=1; % set counters and filenames
0041 bs=p.branch(:,size(p.branch,2)); bvs=p.bifvals(:,size(p.bifvals,2)); 
0042 bs(1)=1; % put last of p on new barnch
0043 if(~isempty(varargin)); [p,ok]=setfn(p,varargin{1}); if ok~=1; q=p; return; end
0044 else fprintf('warning: file name prefixes unchanged.\n'); end
0045 p.branch=bs; p.bifvals=bvs; p.deta=0; 
0046 fname=[p.pname,sprintf('%i',p.count),'.mat'];save(fname,'p'); 
0047 p.count=2;p.bcount=1;q=p;