FreeFEM Documentation on GitHub

stars - forks

Non-linear static problems

Here we propose to solve the following non-linear academic problem of minimization of a functional:

\[J(u) = \int_\Omega \frac{1}{2} f(|\nabla u|^2) - u*b\]

where \(u\) is function of \(H^1_0(\Omega)\) and \(f\) defined by:

\[f(x) = a*x + x-ln(1+x),\ f'(x) = a+\frac{x}{1+x},\ f''(x) = \frac{1}{(1+x)^2}\]

Newton-Raphson algorithm

Now, we solve the Euler problem \(\nabla J (u) = 0\) with Newton-Raphson algorithm, that is:

\[u^{n+1} = u^n - ( \nabla^2 J (u^{n}))^{-1}*\nabla J(u^n)\]
 1 // Parameters
 2 real a = 0.001;
 3 func b = 1.;
 4 
 5 // Mesh
 6 mesh Th = square(10, 10);
 7 Th = adaptmesh(Th, 0.05, IsMetric=1, splitpbedge=1);
 8 plot(Th, wait=true);
 9 
10 // Fespace
11 fespace Vh(Th, P1);
12 Vh u=0;
13 Vh v, w;
14 
15 fespace Ph(Th, P1dc);
16 Ph alpha; //to store |nabla u|^2
17 Ph dalpha ; //to store 2f''(|nabla u|^2)
18 
19 // Function
20 func real f (real u){
21     return u*a + u - log(1.+u);
22 }
23 func real df (real u){
24     return a +u/(1.+u);
25 }
26 func real ddf (real u){
27     return 1. / ((1.+u)*(1.+u));
28 }
29 
30 // Problem
31 //the variational form of evaluate dJ = nabla J
32 //dJ = f'()*(dx(u)*dx(vh) + dy(u)*dy(vh))
33 varf vdJ (uh, vh)
34     = int2d(Th)(
35           alpha*(dx(u)*dx(vh) + dy(u)*dy(vh))
36         - b*vh
37     )
38     + on(1, 2, 3, 4, uh=0)
39     ;
40 
41 //the variational form of evaluate ddJ = nabla^2 J
42 //hJ(uh,vh) = f'()*(dx(uh)*dx(vh) + dy(uh)*dy(vh))
43 //  + 2*f''()(dx(u)*dx(uh) + dy(u)*dy(uh)) * (dx(u)*dx(vh) + dy(u)*dy(vh))
44 varf vhJ (uh, vh)
45     = int2d(Th)(
46           alpha*(dx(uh)*dx(vh) + dy(uh)*dy(vh))
47         + dalpha*(dx(u)*dx(vh) + dy(u)*dy(vh))*(dx(u)*dx(uh) + dy(u)*dy(uh))
48     )
49     + on(1, 2, 3, 4, uh=0)
50     ;
51 
52 // Newton algorithm
53 for (int i = 0; i < 100; i++){
54     // Compute f' and f''
55     alpha = df(dx(u)*dx(u) + dy(u)*dy(u));
56     dalpha = 2*ddf(dx(u)*dx(u) + dy(u)*dy(u));
57 
58     // nabla J
59     v[]= vdJ(0, Vh);
60 
61     // Residual
62     real res = v[]'*v[];
63     cout << i << " residu^2 = " << res << endl;
64     if( res < 1e-12) break;
65 
66     // HJ
67     matrix H = vhJ(Vh, Vh, factorize=1, solver=LU);
68 
69     // Newton
70     w[] = H^-1*v[];
71     u[] -= w[];
72 }
73 
74 // Plot
75 plot (u, wait=true, cmm="Solution with Newton-Raphson");
Table of content