nloopext

 newton loop (with damping parameter alpha) for extended system

0001 function [u,lam,res,iter,Gu,Glam]=nloopext(p,u1,lam1,ds) 
0002 % newton loop (with damping parameter alpha) for extended system
0003 almin=0.01; % minimal damping
0004 u0=p.u; lam0=p.lam; % last point on branch, don't change!
0005 iter=0; u=u1; lam=lam1; alpha=1; 
0006 r=resi(p,u,lam);res0=norm(r,p.normsw); res=res0; % the residual
0007 [Gu,Glam]=getder(p,u,lam,r); % the derivatives
0008 if(res<p.tol); return; end; % initial res. small, do nothing!
0009 
0010 % now start the actual loop
0011 stepok=1; % stepok=1 indicates that step of size alpha decreased residual!
0012 while(abs(res0)>p.tol && iter<p.imax && stepok) 
0013   amat=[[Gu Glam];[p.xi*p.tau(1:p.np*p.neq)',(1-p.xi)*p.tau(p.np*p.neq+1)]];  
0014   p1=p.tau'*[p.xi*(u-u0);(1-p.xi)*(lam-lam0)]-ds; 
0015   upd=p.blss(amat,[r;p1],p,lam); stepok=0; 
0016   while(stepok==0 && alpha>almin) 
0017     u1=u-alpha*upd(1:p.np*p.neq); lam1=lam-alpha*upd(p.np*p.neq+1); 
0018     r=resi(p,u1,lam1); res=norm(r,p.normsw); 
0019     if(res<res0); % good step
0020        if(res<res0/2 && alpha<1) alpha=alpha*2; end % very good step, possibly increase alpha
0021        stepok=1; u=u1; lam=lam1; res0=res;
0022        if(p.nsw==0) [Gu,Glam]=getder(p,u,lam,r); end % full Newton, get new derivatives
0023     else alpha=alpha/2; % bad step, try smaller alpha
0024     end % res<res0
0025   end % while stepok==0
0026   iter=iter+1;
0027 end % while res<p.tol && stepok
0028 
0029 % postprocessing
0030 if(p.nsw>0);  % if chord, then now get derivatives at new point!
0031     [Gu,Glam]=getder(p,u,lam,r); end
0032 if(p.vsw>1); if(alpha<1) fprintf('\nnloopext: damp alpha=%g, res=%g, ds=%g\n',...
0033             alpha,res,ds); end; end; % inform user if damping was used
0034 res=res0;