/usr/share/doc/freefem++/examples/examples++-tutorial/Newton.edp is in freefem++-doc 3.19.1-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 | { // --- a real non linear test ---
mesh Th=square(10,10); // mesh definition of $\Omega$
Th = adaptmesh(Th,0.05,IsMetric=1,splitpbedge=1);
//plot(Th,wait=1);
//Th = adaptmesh(Th,0.1,IsMetric=1,splitpbedge=1);
plot(Th,wait=0);
fespace Vh(Th,P1); // finite element space
fespace Ph(Th,P1dc); // make optimization
Vh b=1; // to defined b
// $ J(u) = 1/2 \int_\Omega f(|\nabla u|^2) - \int\Omega u b $
// $ f(x) = a*x + x-ln(1+x), \quad f'(x) = a+\frac{x}{1+x}, \quad f''(x) = \frac{1}{(1+x)^2}$
real a=0.001;
func real f(real u) { return u*a+u-log(1+u); }
func real df(real u) { return a+u/(1+u);}
func real ddf(real u) { return 1/((1+u)*(1+u));}
Vh u=0; // the current value of the solution
Ph alpha; // to store $ (|\nabla u|^2)$
Ph dfalpha; // to store $f' (|\nabla u|^2)$
Ph ddfalpha ; //to store = $2 f''( |\nabla u|^2) $ optimisation
int iter=0;
// methode of Newton Ralphson to solve dJ(u)=0;
// $$ u^{n+1} = u^n - (\frac{\partial dJ}{\partial u_i})^{-1}*dJ(u^n) $$
// ---------------------------------------------
// the variationnal form of evaluate dJ
// --------------------------------------
// dJ = f'()*( dx(u)*dx(vh) + dy(u)*dy(vh)
varf vdJ(uh,vh) = int2d(Th)( dfalpha*( dx(u)*dx(vh) + dy(u)*dy(vh) ) - b*vh)
+ on(1,2,3,4, uh=0);
// the variationnal form of evaluate ddJ
// hJ(uh,vh) = f'()*( dx(uh)*dx(vh) + dy(uh)*dy(vh)
// + f''()( dx(u)*dx(uh) + dy(u)*dy(uh) ) * (dx(u)*dx(vh) + dy(u)*dy(vh))
varf vhJ(uh,vh) = int2d(Th)( dfalpha*( dx(uh)*dx(vh) + dy(uh)*dy(vh) )
+ ddfalpha*( dx(u)*dx(vh) + dy(u)*dy(vh) )*( dx(u)*dx(uh) + dy(u)*dy(uh) ) )
+ on(1,2,3,4, uh=0);
// the newton algorithm
Vh v,w;
u=0;
for (int i=0;i<100;i++)
{
alpha = dx(u)*dx(u) + dy(u)*dy(u);// optimization
dfalpha = df( alpha ) ; // optimization
ddfalpha = 2*ddf(alpha ) ; // optimization
v[]= vdJ(0,Vh);
real res= v[]'*v[];
cout << i << " residu^2 = " << res << endl;
matrix H;
if( res< 1e-12) break;
H= vhJ(Vh,Vh,factorize=1,solver=LU);
w[]=H^-1*v[];
u[] -= w[];
plot (u,wait=0,cmm="solution with Newton Raphson");
}
savemesh(Th,"N",[x,y,u*1.5]);
{ ofstream file("N.bb");
file << "2 1 1 "<< u[].n << " 2 \n";
int j;
for (j=0;j<u[].n ; j++)
file << u[][j] << endl;
}
exec("ffmedit N");
exec("rm N.*");
}
|