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
|
// FFCS: making a numerical value visible for regression tests
real regtest;
{ // --- 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");
// FFCS: regression tests
regtest=u[]'*u[]; //'
}
load "medit" load "msh3"
mesh3 Th3= movemesh23(Th,transfo=[x,y,u*1.5]);
medit("N",Th3);
}
|
