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
|
// file optimcontrol.edp
border aa(t=0, 2*pi) { x = 5*cos(t); y = 5*sin(t); };
border bb(t=0, 2*pi) { x = cos(t); y = sin(t); };
border cc(t=0, 2*pi) { x = -3+cos(t); y = sin(t); };
border dd(t=0, 2*pi) { x = cos(t); y = -3+sin(t); };
mesh th = buildmesh(aa(70)+bb(35)+cc(35)+dd(35));
fespace Vh(th,P1);
Vh Ib=((x^2+y^2)<1.0001),
Ic=(((x+3)^2+ y^2)<1.0001),
Id=((x^2+(y+3)^2)<1.0001),
Ie=(((x-1)^2+ y^2)<=4),
ud,u,uh,du;
real[int] z(3);
problem A(u,uh) =int2d(th)((1+z[0]*Ib+z[1]*Ic+z[2]*Id)*(dx(u)*dx(uh)
+dy(u)*dy(uh))) + on(aa,u=x^3-y^3);
z[0]=2; z[1]=3; z[2]=4;
A; ud=u;
ofstream f("J.txt");
func real J(real[int] & Z)
{
for (int i=0;i<z.n;i++)z[i]=Z[i];
A; real s= int2d(th)(Ie*(u-ud)^2);
f<<s<<" "; return s;
}
real[int] dz(3), dJdz(3);
problem B(du,uh)
=int2d(th)((1+z[0]*Ib+z[1]*Ic+z[2]*Id)*(dx(du)*dx(uh)+dy(du)*dy(uh)))
+int2d(th)((dz[0]*Ib+dz[1]*Ic+dz[2]*Id)*(dx(u)*dx(uh)+dy(u)*dy(uh)))
+on(aa,du=0);
func real[int] DJ(real[int] &Z)
{
for(int i=0;i<z.n;i++)
{ for(int j=0;j<dz.n;j++) dz[j]=0;
dz[i]=1; B;
dJdz[i]= 2*int2d(th)(Ie*(u-ud)*du);
}
return dJdz;
}
real[int] Z(3);
for(int j=0;j<z.n;j++) Z[j]=1;
BFGS(J,DJ,Z,eps=1.e-6,nbiter=15,nbiterline=20);
cout << "BFGS: J(z) = " << J(Z) << endl;
for(int j=0;j<z.n;j++) cout<<z[j]<<endl;
plot(ud,value=1,ps="u.eps");
|
