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
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
|
// cleanning version 07/2008 FH in Sevilla.
int nerr =0;
int debugJ =0;
int debugdJ =0;
real umax=0;
{
func bool stop(int iter,real[int] u,real[int] g)
{
cout << " stop = " << iter << " " << u.linfty << " " << g.linfty << endl;
return g.linfty < 1e-5 || iter > 15;;
}
// minimisation of $J(u) = \frac12\sum (i+1) u_i^2 - b_i $
// work array
real[int] b(10),u(10);
func real J(real[int] & u)
{
real s=0;
for (int i=0;i<u.n;i++)
s +=(i+1)*u[i]*u[i]*0.5 - b[i]*u[i];
if(debugJ) cout << "J ="<< s << " u =" << u[0] << " " << u[1] << "...\n" ;
return s;
}
// the grad of J (this is a affine version (the RHS is in )
func real[int] DJ(real[int] &u)
{
for (int i=0;i<u.n;i++)
u[i]=(i+1)*u[i];
if(debugdJ) cout << "dJ0 ="<< " u =" << u[0] << " " << u[1] << " " << u[2] << "...\n" ;
u -= b;
if(debugdJ) cout << "dJ-b ="<< " u =" << u[0] << " " << u[1] << " " << u[2] << "...\n" ;
return u; // return of global variable ok
};
// the grad of the bilinear part of J (the RHS in remove)
func real[int] DJ0(real[int] &u)
{
for (int i=0;i<u.n;i++)
u[i]=(i+1)*u[i];
if(debugdJ) cout << "dJ0 ="<< " u =" << u[0] << " " << u[1] << " " << u[2] << "...\n" ;
return u; // return of global variable ok
};
func real error(real[int] & u,real[int] & b)
{
real s=0;
for (int i=0;i<u.n;i++)
s += abs((i+1)*u[i] - b[i]);
return s;
}
func real[int] matId(real[int] &u) { return u;};
int verb=5; // verbosity of algo ..
b=1. ; u=0.; // set right hand side and initial gest
LinearCG(DJ,u,eps=1.e-6,nbiter=20,precon=matId,verbosity=verb);
cout << "LinearGC (Affine) : J(u) = " << J(u) << " err=" << error(u,b) << endl;
nerr += !(error(u,b) < 1e-5);
if(nerr) cout << " sol: "<< " u =" << u[0] << " " << u[1] << " " << u[2] << "...\n" ;
b=1; u=0;
LinearCG(DJ,u,eps=1.e-15,nbiter=20,precon=matId,verbosity=50,stop=stop);
cout << "LinearGC (Affine with stop) : J(u) = " << J(u) << " err=" << error(u,b) << endl;
nerr += !(error(u,b) < 1e-5);
if(nerr) cout << " sol: "<< " u =" << u[0] << " " << u[1] << " " << u[2] << "...\n" ;
b=1; u=0; // set right hand side and initial gest
LinearCG(DJ0,u,b,eps=1.e-6,nbiter=20,precon=matId,verbosity=verb);
cout << "LinearGC (Linear) : J(u) = " << J(u) << " err=" << error(u,b) << endl;
nerr += !(error(u,b) < 1e-5);
if(nerr) cout << " sol: "<< " u =" << u[0] << " " << u[1] << " " << u[2] << "...\n" ;
b=1; u=0; // set right hand side and initial gest
AffineGMRES(DJ,u,eps=1.e-6,nbiter=20,precon=matId,verbosity=verb); // correct in version 3.11
cout << "LinearGMRES (Affine) : J(u) = " << J(u) << " err=" << error(u,b) << endl;
nerr += !(error(u,b) < 1e-5);
if(nerr) cout << " sol: "<< " u =" << u[0] << " " << u[1] << " " << u[2] << "...\n" ;
b=1; u=0; // set right hand side and initial gest
LinearGMRES(DJ0,u,b,eps=1.e-6,nbiter=20,precon=matId,verbosity=verb);
cout << "LinearGMRES (Linear) : J(u) = " << J(u) << " err=" << error(u,b) << endl;
nerr += !(error(u,b) < 1e-5);
if(nerr) cout << " sol: "<< " u =" << u[0] << " " << u[1] << " " << u[2] << "...\n" ;
b=1; u=0; // set right hand side and initial gest
NLCG(DJ,u,eps=1.e-6,nbiter=20,precon=matId,verbosity=verb);
cout << "NLCG: J(u) = " << J(u) << " err=" << error(u,b) << endl;
nerr += !(error(u,b) < 1e-5);
if(nerr) cout << " sol: "<< " u =" << u[0] << " " << u[1] << " " << u[2] << "...\n" ;
// warning BFGS use a full matrix of size nxn (where n=u.n)
b=1; u=2; // set right hand side and initial gest
BFGS(J,DJ,u,eps=1.e-6,nbiter=20,nbiterline=20);
cout << "BFGS: J(u) = " << J(u) << " err=" << error(u,b) << endl;
assert(error(u,b) < 1e-5);
if(nerr) cout << " sol: "<< " u =" << u[0] << " " << u[1] << " " << u[2] << "...\n" ;
assert(nerr==0);
};
{ // --- a real non linear test ---
mesh Th=square(10,10); // mesh definition of $\Omega$
fespace Vh(Th,P1); // finite element space
fespace Ph(Th,P0); // make optimization
Vh b=1; // to defined b
// $ J(u) = 1/2\int_\Omega f(|\nabla u|^2) - \int\Omega u b $
// $ f(u) = a*u + u-ln(1+u), \quad f'(u) = a+\frac{u}{1+u}, \quad f''(u) = \frac{1}{(1+u)^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));}
// the functionnal J
func real J(real[int] & u)
{
Vh w;w[]=u;
real r=int2d(Th)(0.5*f( dx(w)*dx(w) + dy(w)*dy(w) ) - b*w) ;
cout << "J(u) =" << r << " " << u.min << " " << u.max << endl;
return r;
}
// -----------------------
Vh u=0; // the current value of the solution
Ph alpha; // of store $df(|\nabla u|^2)$
int iter=0;
alpha=df( dx(u)*dx(u) + dy(u)*dy(u) ); // optimization
func real[int] dJ(real[int] & u)
{
Vh w;w[]=u;
alpha=df( dx(w)*dx(w) + dy(w)*dy(w) ); // optimization
varf au(uh,vh) = int2d(Th)( alpha*( dx(w)*dx(vh) + dy(w)*dy(vh) ) - b*vh)
+ on(1,2,3,4,uh=0);
u= au(0,Vh);
return u; // warning no return of local array
}
varf alap(uh,vh)=
int2d(Th)( alpha *( dx(uh)*dx(vh) + dy(uh)*dy(vh) )) + on(1,2,3,4,uh=0);
varf amass(uh,vh)= int2d(Th)( uh*vh) + on(1,2,3,4,uh=0);
matrix Amass = amass(Vh,Vh,solver=CG);
matrix Alap= alap(Vh,Vh,solver=Cholesky,factorize=1);
func real[int] C(real[int] & u)
{
real[int] w=Amass*u;
u = Alap^-1*w;
return u; // no return of local array variable
}
int conv=0;
real eps=1e-6;
for(int i=0;i<20;i++)
{
conv=NLCG(dJ,u[],nbiter=10,precon=C,veps=eps,verbosity=5);
if (conv) break;
alpha=df( dx(u)*dx(u) + dy(u)*dy(u) ); // optimization
Alap = alap(Vh,Vh,solver=Cholesky,factorize=1);
cout << " restart with new preconditionner " << conv << " eps =" << eps << endl;
}
plot (u,wait=1,cmm="solution with NLCG");
umax = u[].max;
Vh sss= df( dx(u)*dx(u) + dy(u)*dy(u) ) ;
plot (sss,wait=0,fill=1,value=1);
// the method of Newton Ralphson to solve dJ(u)=0;
// see Newton.edp example
}
assert(nerr==0);
|
