Boundary condition in Newton method

General Discussion
1 
         Hi everybody,

I search how to put the boundary condition in the Newton method , for Dirichlet , Neumann or Robin condition.

For example, solving the toy equation d2 u(x) /dx2 + u(x) -1 =0, u(0)=1, u(1)=10. It’s just to try the program. Here limits are added to each loop, so it’s an error.

fespace Vh(Th,P1);
Vh u, v, du, dv;

varf right(u, v) = int2d(Th)( -dx(u) * dx(v) +v * u ) -int2d(Th)( v*1.0 )

  • on(2,u=10 ) + on(4,u=1);

varf left(du , dv ) = int2d(Th)( -dx(du) *dx(dv) +dv *du )

  • on(2,du=0 ) + on(4,du=0) ;

for (int n = 0; n < 6; n++){

real[int] b= right(0, Vh) ;
matrix a=left(Vh, Vh );

du[]= a^-1 * b ;
u[] += du[];

}

Thank you
2

Some remark on Newton method and Boundary condition of every kind

You solve F(u) = 0
let call DF (u)v = the derivative in direction v
The Newton method is

let call w^n solution of $ DF(u^n) w = F(u^n)$ ;
when $ u^{n=1} = u^n - w^n $
and the Boundary condition on w^n must be homogenous
now remember $DF(u^n) (w)$ is linear in so we can commute directly $u^{n+1} u^n - w^n $
so

DF(u^n) u^{n+1} = DF(u^n) u^{n} - F(u^n)
en the B.C on $u^{n+1}$ are the given B.C and the unitialisation is more simple (without B.C) in
general:

exemple
mesh Th=square(10,10); // mesh definition of $\Omega$
fespace Vh(Th,P1); // finite element space
macro grad(u) [dx(u),dy(u)]//
macro F(u,v) (int2d(Th)( grad(u)'grad(v) - u^3v) + int2d(Th)((x+y)v))// + on(1,2,3,4,u=1);
// the probleme is Find u : F(u,v) = 0 + on(1,2,3,4,u=1);
macro DF(u,v,w) int2d(Th)( grad(w)'grad(v) - 3u^2w*v) // + on(1,2,3,4,u=1);
Vh u=0, up,v ; // the current value of the solution
for (int iter=1; iter < 100; ++iter)
{
up[]=u[];
solve PbNewton(u,v) = DF(up,v,u) - DF(up,v,up) + F(up,v) + on(1,2,3,4,u=1);
up[] -=u[]; // to compute the err
real err = up[].linfty ;
cout << iter << " " << err << endl;
plot(u);
if(err<1e-6) break;
}

3 
           Hi,

Thanks you very much Professor Hecht.
If I understand correctly, we transform the equation into a sequence , u^n → u^(n+1) , w^n is no more used.

I would like to share an another toy equation in 1d :
y’(x) -y(x)^2 -1 =0 , which have a analytical solution, y(x) = tan(x + constant)

PbNewton(u,v) = int2d(Th)( dx(u)v - 2uupv ) -int2d(Th)( dx(up)v - 2upupv )
+int2d(Th)( dx(up)v - up^2v -v )
+on(4, u= 0.2 ) ; // the B.C

But I have a difficulty for Neumann condition.
If I take y(x)=tan(x), cst=0 for simplicity, y’(x)=1+tan² (x) , y’(0)=1
and change the last line by :
-int1d(Th,4) ( 1.0v ) + int2d(Th)( 1E-6u*v)

It don’t work, with or without the penalty term. Someone could remind me how to put cleanly the B.C ?

             Sincerely
4

You want to solve u’(x) -u(x)^2 -1 =0 on ]0,1[, u’(0)=1 but you have not good variational form
The problem not a pb of Newton but a Pb of formulation in FreeFem++``

Sorry no idea your problem is not a good problem for variational formulation.