Startpage >> Main >> Newton2

Newton 2

Newton method for \(F:\mathbb{R}\to \mathbb{R}\) is given by:

Want \(x_\star\) such that \(F(x_\star)=0\). Given \(x \approx x_\star\), by Taylor expansion we have \(F(x_\star)=0\approx F(x)+(x_\star-x)F'(x)\) . Then, we improve \(x\) by taking \(x:=x+h\) where \(h=x_\star-x\) satisfies \begin{array}{rcl} F'(x)h & =& -F(x). \end{array}

Take now \(F:U\to V\) where \(U\) and \(V\) are Banach spaces. The Gateaux derivative of \(F\) at a point \(x_\star\in U\) is the linear mapping \(dF(x_\star):U\to V\) that satisfies (denote the value of \(dF(x_\star) \) at \(h\in U\) by \(<dF(x_\star),h>\in V\) )

$$ \begin{array}{rcl} <dF(x_\star),h> & :=&\lim_{s\to 0}\frac{1}{s}(F(x_\star+sh)-F(x_\star)),\quad \forall h\in U. \end{array} $$

Newton method for \(F:U\to V\) is given by:

Want \(x_\star\in U\) such that \(F(x_\star)=0\). Given \(x \approx x_\star\) improve \(x\) by taking \(x:=x+h\) where \(h\) satisfies \begin{array}{rcl} <dF(x),h> & =& -F(x). \end{array}

Example 1. For the nonlinear advection equation on \(\Omega\) with \(\partial\Omega=\Gamma_D\cup\Gamma_N\)

\begin{array}{rcl} \alpha\vec{u}-\mu \Delta \vec{u}+ (\vec{u}\cdot \vec{\nabla})\vec{u} & =& \vec{f},\quad \Omega\\ \vec{u} & =& 0,\quad \Gamma_D\\ \frac{\partial{\vec{u}}}{\partial\vec{n}} & =& 0,\quad \Gamma_N \end{array}

include Dirichlet boundary conditions on \(U\) and take

$$F(\vec{u}):=\alpha\vec{u}-\mu \Delta \vec{u}+ (\vec{u}\cdot \vec{\nabla})\vec{u} -\vec{f}$$

a simple computation gives

$$<dF(\vec{u}),\vec{h}>=\alpha\vec{h}-\mu \Delta \vec{h}+ (\vec{u}\cdot \vec{\nabla})\vec{h} +(\vec{h}\cdot \vec{\nabla})\vec{u} $$

and Newton method is

Want \(\vec{u}_\star\in U\) such that \(F(\vec{u}_\star)=0\). Given \(\vec{u} \approx \vec{u}_\star\) improve \(\vec{u}\) by taking \(\vec{u}:=\vec{u}+\vec{h}\) where \(\vec{h}\) satisfies the linear problem \begin{array}{rcl} \alpha\vec{h}-\mu \Delta \vec{h}+ (\vec{u}\cdot \vec{\nabla})\vec{h} +(\vec{h}\cdot \vec{\nabla})\vec{u} & =& -(\alpha\vec{u}-\mu \Delta \vec{u}+ (\vec{u}\cdot \vec{\nabla})\vec{u} -\vec{f}). \end{array}

Example 2. For the nonlinear Navier-Stokes equation on \(\Omega\) with \(\partial\Omega=\Gamma_D\cup\Gamma_N\)

\begin{array}{rcl} -\mu \Delta \vec{u}+ (\vec{u}\cdot \vec{\nabla})\vec{u} +\vec{\nabla p}& =& \vec{f},\quad \Omega\\ div(\vec{u})& =& 0,\quad \Omega\\ \vec{u} & =& 0,\quad \Gamma_D\\ \frac{\partial{\vec{u}}}{\partial\vec{n}} & =& 0,\quad \Gamma_N \end{array}

include Dirichlet boundary conditions on \(U\) space for the pair \((\vec{u},p)\) and maintain the restriction \(div( \vec{u})=0\). Then take

$$F(\vec{u},p):=(-\mu \Delta \vec{u}+ (\vec{u}\cdot \vec{\nabla})\vec{u}+\vec{\nabla p} -\vec{f},\; div(\vec{u}))$$

again, a simple computation gives

$$<dF((\vec{u},p)),(\vec{h},q)>=(-\mu \Delta \vec{h}+ (\vec{u}\cdot \vec{\nabla})\vec{h} +(\vec{h}\cdot \vec{\nabla})\vec{u}+\vec{\nabla q},\; div(\vec{h})) $$

and Newton method is

Want \((\vec{u}_\star,p_\star)\in U\) such that \(F(\vec{u}_\star,p_\star)=0\). Given \((\vec{u},p) \approx (\vec{u}_\star,p_\star)\) improve \((\vec{u},p)\) by taking \(\vec{u}:=\vec{u}+\vec{h}\) and \(p:=p+q\) where \((\vec{h},q)\) satisfies the linear problem \begin{array}{rcl} (-\mu \Delta \vec{h}+ (\vec{u}\cdot \vec{\nabla})\vec{h} +(\vec{h}\cdot \vec{\nabla})\vec{u}+\vec{\nabla q}, \; div(\vec{h}) ) & =& (-\mu \Delta \vec{u}+ (\vec{u}\cdot \vec{\nabla})\vec{u}+\vec{\nabla p} -\vec{f},\; div(\vec{u})). \end{array} Solved by a variational formulation using a pair of test function \((\vec{v} ,s)\).

Incompressible Navier Stokes around a 3d Cylinder. Continuation on viscosity

Adapted on April 2014 from code by Author: F. Hecht, jan 2012.

TriangulationViscosity 0.02Viscosity 0.002
Close up viscosity 0.02Close up viscosity 0.002

download example: NSNewtonv1.edp or return to main page

// 
// adapted on April 2014 from code by Author: F. Hecht  
// jan 2012 
//
//Stationnary imcompressible Navier Stokes Equation with Newton method.
//  a round a 3d Cylinder 
// build the  Mesh
//

real R = 5,L=15;
border cc(t=0,2*pi){ x=cos(t)/2;y=sin(t)/2;label=1;}
border ce(t=pi/2,3*pi/2) { x=cos(t)*R;y=sin(t)*R;label=1;}
//cluster triangulation points on boundary y=-R near x=0
border beb(tt=0,1) { real t=tt^1.2; x= t*L; y= -R; label = 1;}
//cluster triangulation points on boundary y=R near x=0
border beu(tt=1,0) { real t=tt^1.2; x= t*L; y= R; label = 1;}
border beo(t=-R,R) {  x= L; y= t; label = 0;}
/*
//
// to see the boundary curves one by one
//
plot(cc(-50),wait=1);
plot(cc(-50)+ce(30),wait=1);
plot(cc(-50),ce(30),beb(20),wait=1);
plot(cc(-50),ce(30),beb(20),beu(20),wait=1);
plot(cc(-50),ce(30),beb(20),beu(20),beo(10),wait=1);
*/
int mm=30;
mesh Th=buildmesh(cc(-2*mm)+ce(3*mm)+beb(3*mm)+beu(3*mm)+beo(2*mm));
plot(cmm=" This is the mesh...Wait for computations...",Th);
//
// for a close view
//
func bb=[[-1,-2],[4,2]];
//
// Macro operators
//
macro Grad(u1,u2) [ dx(u1),dy(u1), dx(u2),dy(u2)]// 
macro UgradV(u1,u2,v1,v2) [ [u1,u2]'*[dx(v1),dy(v1)] , [u1,u2]'*[dx(v2),dy(v2)] ]// 
macro div(u1,u2)  (dx(u1)+dy(u2))//

//  FE Space 
fespace Xh(Th,P2);// Use P1b for a less expensive computation
fespace Mh(Th,P1);

Xh u1,u2;// solutions for velocity variables
Xh du1,du2;// increments for velocity variables
Xh uc1,uc2;// computed solutions for velocity variables at previous viscosity
Mh p;// solution for the pressure variable
Mh dp;// increment for the pressure variable
Xh v1,v2; // test functions for the velocity variables
Mh q;// test function for pressure variable
Mh pc;//  computed solutions for pressure variable at previous viscosity

// intial guess with B.C. 

u1 = ( x^2+y^2) > 2;
u2=0;
p=0;
//  Physical parameter 
real nu= 1./50; // starting viscosity
real nufinal=1/500.; // wanted final viscosity
//
// cnu reducing factor for viscosity
// We reduce viscosity by dividing by cnu
//
real cnu=2; 
//
// Number of reductions on viscosity
//
int countermax=10;
int counter=0;
// stop test for Newton 
real eps=1e-6;

//
verbosity=0;
real err=0;
real nuc=nu;// current viscosity
cout <<  "  ----------------------------  " << endl;
cout <<  "  --- START COMPUTATIONS ---  " << endl;
cout <<  "  ----------------------------  " << endl;
//
//  Loop on vicosity until reach nufinal, countermax computations or stop reducing 
//
while((nuc>nufinal) && (counter<countermax)&&(cnu>1.05))  
{ // start loop on viscosity 
       counter =counter+1;
	int n; // need value of n out of the for loop
	cout <<  " Compute with viscosity  " << nuc <<" (counter "<<counter<<")"<< endl;
	for( n=0;n< 15;n++) // Newton Loop for current value of cnu
	  { // --->star for loop
	    //   compute increments for the three variables u1,u2,p
	     solve Oseen([du1,du2,dp],[v1,v2,q],init=n) =
		   int2d(Th) (  nuc*(Grad(du1,du2)'*Grad(v1,v2) )
		              + UgradV(du1,du2, u1, u2)'*[v1,v2]
		              + UgradV( u1, u2,du1,du2)'*[v1,v2]
		              - div(du1,du2)*q - div(v1,v2)*dp 
		              - 1e-8*dp*q // stabilization term 
		             )
		+ int2d(Th) (  nuc*(Grad(u1,u2)'*Grad(v1,v2) )
		              + UgradV(u1,u2, u1, u2)'*[v1,v2]
		              - div(u1,u2)*q - div(v1,v2)*p 
		              - 1e-8*p*q 
		             )
	    + on(1,du1=0,du2=0)
	    ;

//        
// update values u1, u2,p to compute again
//
	    u1[] += du1[]; //u1[]=u1[]+du1[]
	    u2[] += du2[]; //u2[]=u2[]+du2[]
	     p[]  += dp[];  //p[]=p[]+dp[]

//
// computation of norms for the increments
//
	    real Lu1=u1[].linfty,  Lu2 = u2[].linfty , Lp = p[].linfty;
//
// normalization of the norms
//
    err= du1[].linfty/Lu1 + du2[].linfty/Lu2 + dp[].linfty/Lp;
    cout<<" it = "+ n+ ", err = "+err <<endl;
//
// plot iterates 
//   
 plot([u1,u2],p,wait=0,cmm="Newton Method for viscosity = " + nuc +", it = "+ n+ ", err = "+err, coef=0.3);
//
//--> In case of no convergence...
//      (n>3 AND err>2) is blowup 
//      (n>12 AND err>eps) is lack of convergence or very slow
//
   if( (n>3 && err > 2.) || (n>12 && err > 1e4*eps) )  // (n>3 AND err>2) or (n>12 AND err>eps)       
   {// start ---> if
    cout<<" NO convergence for nu = "+nuc+ ", it = "+ n+ ", err = "+err <<endl;
    nuc=nuc*cnu; // get back to previous nu
    cout << " LAST viscosity TO CONVERGE = " << nuc<< endl;
    //
    // take computed values for previous nuc to start again
    //
    u1=uc1; u2=uc2; p=pc;
    //
    // New reduction for cnu
    //
    cnu= cnu^(0.6); // no conv. => change lower (cnu^(0.6) makes cnu smaller)
    nuc = nuc/cnu;  // new vicosity
    cout << " RESTART with NEW viscosity  = " << nuc << endl;
    cout << " Using as Inital Data the solution for PREVIOUS viscosity "  << endl;
    break; //  Blowup ???? --->out of  for loop  go to while loop 
   }// end ---> if( n>3 && err > 2.) 
 //
//-->In case of Convergence then...	
// 
   if(err < eps) 
   {// start ---> if(err < eps)
  cout << " CONVERGENCE  for viscosity = " << nuc << " on iters = "<< n <<", err = "+err << endl;
    uc1=u1; uc2=u2; pc=p; // save computed values of this convergent solution
    plot([u1,u2],p,wait=1,cmm="Solution for viscosity = " + nuc+". Press Enter for a close up ",coef=0.3);
    plot([u1,u2],p,wait=1,coef=0.3,bb=bb);
     if( n < 4)  cnu=cnu^1.5; // fast converge => change faster ( cnu^ 1.5 makes bigger cnu)   
      nuc = max(nufinal, nuc/cnu); // reduce to a new vicosity 
      cout << " Continuation to new viscosity  = " << nuc <<" (dividing factor "<<cnu<<")"<< endl;
      cout << " Initial data IS the FINAL computation  for previous viscosity " << endl;
   break; // converge ---> out of for loop. Start again with  u1,u2,p and new nuc
   }// end ---> if(err < eps)
 }//end ---> for loop
} // end  while loop on viscosity 
plot([uc1,uc2],pc,cmm="Last convergence of Newton Method for nu = " + nuc+". Press Enter for a close up ",coef=0.3,wait=1);
plot([uc1,uc2],pc,cmm="Close up nu = " + nuc,coef=0.3,bb=bb);
cout << " LAST viscosity TO CONVERGE nu = " << nuc << endl;

return to main page

Page last modified on April 15, 2014, at 06:04 PM