I’m trying to estimate the residual error when solving a pde with solve. More precisely freefem is solving a system which looks like : find U st AU = B where A is a matrix, B a vector. I would like to compute |AU-B|/|B|. I come up with the following code, but the results are far from the one i expected. Could someone explain to me what I’ve been doing wrong ? I suspect that i misunderstood how to use varf…

Thanks a lot,

Simon

real h = 0.01;
int Nm=floor(1./h); //Number of FE nodes on a unit segment
mesh Th=square(Nm,Nm,[x,y]);
fespace Vh(Th,P1);
Vh u,v,f;
f=1;
{
solve VA(u,v, solver=sparsesolver)=
int2d(Th)(dx(u)*dx(v) + dy(u)*dy(v))
- int2d(Th)(f*v)
+ on(1, 2, 3, 4, u=0);
}
varf bilinear(u,v) = int2d(Th)(dx(u)*dx(v) + dy(u)*dy(v))
+ on(1, 2, 3, 4, u=0);
varf linear(unused,v) = int2d(Th)(f*v);
matrix A = bilinear(Vh, Vh);
Vh U, B;
U[] = A*u[];
B[] = linear(0, Vh);
U[] -= B[];
cout << "residuals (expected small): " << ", " << sqrt(u[]'*u[])/sqrt(B[]'*B[]) << endl;

Indeed thanks for the help !
The issue remains though : for this code example, the relative error is about 10% (which i find quite big, but i may be wrong).

I see but i’m looking for the error stemming from the resolution of the matricial system AU=B, where U is the unknown and represents the coefficients in the decomposition of the solution of the PDE. I see U as a big vector of R^n which contains the residual terms a(u,\phi_i)-b(\phi_i) for all i (vertex of the mesh) and want to sum (in a L2 sense) these local errors to get a “macroscopic” error.
This error is obviously linked to the one you are mentioning, but i guess this is not the one I’m looking for.