// p-Laplace problem, p>2, Rayleigh-Ritz int nint = 100; real p = 4.0; // p > 2.0 otherwise NLCG input "too" exact. Do not forget .0 real r = 1.0; border Gamma (t=0,2*pi) { x = r * cos(t); y = r * sin(t); label = 1; } mesh Th = buildmesh(Gamma(nint)); fespace Vh(Th,P1); Vh u0,u,v; func f = 1.0; func uexact = (p-1)/p*(0.5)^(1/(p-1))*(1-(x^2+y^2)^(p/(2*p-2))); real Juexact = 2*pi*(0.5^(p/(p-1.0))*(p-1)/(3*p-2.0)/p- (p-1)/p*0.5^(1.0/(p-1))*(0.5-(p-1)/(3*p-2))); func real J(real[int] & ucoeff) { Vh uloc; uloc[] = ucoeff; return int2d(Th)((dx(uloc)^2+dy(uloc)^2)^(p/2)/p-f*uloc); } func real[int] deltaJR(real[int] & ucoeff) // gradient of J^R applied to each phi function { Vh uloc; uloc[] = ucoeff; varf F(unused,v) = int2d(Th)((dx(uloc)^2+dy(uloc)^2)^((p-2)/2)* (dx(uloc)*dx(v)+dy(uloc)*dy(v))-f*v) +on(1,unused=0); // it takes value 0 at points on the boundaries /* In NLCG, x_m = x_{m-1} + \alpha_{m-1} d_{m-1} d_{m-1} is a composition of nablaJR directions Therefore, if x_0 has the right boundary conditions, we should not touch them */ real[int] r = F(0,Vh); return r; } // initial solution given by Poisson (p=2) solve Poisson(u0,v) = int2d(Th)((dx(u0)*dx(v)+dy(u0)*dy(v))) // bilinear form -int2d(Th)(f*v) // linear functional +on(1,u0=0); // boundary conditions /* The default solver is sparsesolver (multi-frontal LU) It is possible to specify solver=CG, and possibly eps */ // u0 is the solution of Poisson's problem // simple solution, no preconditioner u = u0; NLCG(deltaJR,u[],nbiter=200,eps=-1e-4); // a simple minimizer cout << "L2 error " << sqrt(int2d(Th)((uexact-u)^2)) << endl; cout << "J error " << abs(J(u[])-Juexact) << endl; plot(u);