# Large residual error when solving a 2D BVP with NDSolve and Finite Element Method

I'm trying to solve a simple PDE (an inhomogeneous Helmholtz equation) in a disk domain with center at $(1,0)$ and radius 0.5 using the following commands

NDSolve[{Laplacian[u[x,y],{x,y}]== -0.010 u[x,y] + 0.01 + 0.05 x^2, DirichletCondition[u[x, y] == 0, True]},u, {x, y}\[Element] Disk[{1, 0}, 0.5], Method ->{"PDEDiscretization" -> {"FiniteElement","MeshOptions" -> {"MaxCellMeasure" -> 0.00001},"IntegrationOrder" -> 5}}, "InterpolationOrder" -> 5];

p[x_, y_] = %[[1, 1, 2]][x, y];

errp[x_, y_] = Laplacian[p[x,y],{x,y}]- (-0.01 p[x, y] + 0.01 + 0.05 x^2);

Plot[{errp[x, 0], p[x, 0]}, {x, 0.5, 1.5}]


It seems that the residual error is great compared to the computed solution. Is there any improvement in order to overcome this problem?

• The Laplacian of the soution of FEM discretization can only converge weakly (in $H^{-1}$to the right hand side. High oscillations are allowed for weak convergence. Jul 25 '18 at 12:13
• The list of arguments in errp[r_, z_]=... seems to be wrong? Jul 25 '18 at 12:30
• @UlrichNeumann yes right! I corrected the arguments, it was just a typo...
– DK13
Jul 25 '18 at 12:34
• @DK13: I (MMA 11.0.1.0 Windows )cannot reproduce your result. Which Mathematica version you are using? Jul 25 '18 at 13:04
• version 11.1 Linux...
– DK13
Jul 25 '18 at 13:25

This is expected. The mesh NDSolve uses is of second order and you will not get a good approximation if you apply the Laplacian to that. The solution would be bit better for for third order meshes, which we do not have in version V11.3.

Now, to verify PDE solutions one use a process called manufactured solutions. For that you take a PDE and plug in an arbitrary function. The result of that will be a right hand side. Next you compute the solution to the PDE with that right hand side just computed. Now you can check the quality of the solution found by looking at the difference between the solution found and the function you plugged into the PDE. In code:

pde = Laplacian[
u[x, y], {x, y}] - (-0.010 u[x, y] + 0.01 + 0.05 x^2);
uExact = Function[{x, y}, x^2 + y^2];
f = pde /. u -> uExact
(* 3.99 - 0.05 x^2 + 0.01 (x^2 + y^2) *)
r = Disk[{1, 0}, 0.5];
ufun = NDSolveValue[{pde == f,
DirichletCondition[u[x, y] == uExact[x, y], True]},
u, {x, y} \[Element] r];
Plot3D[{ufun[x, y] - uExact[x, y]}, Element[{x, y}, r]]


Proof for the claim above:

ufun["ElementMesh"]["MeshOrder"]
2
`
• well what does a small error resulting from the method of manufactured solutions mean since the error of the original problem remains large?
– DK13
Jul 25 '18 at 13:51
• It means that the FEM works. Recovering a derivative from an interpolation function is hard already, recovering a second derivative (like a Laplacian) from that is even harder. So is the error you see from FEM or from a hard recovery of the derivatives? Jul 25 '18 at 14:04
• Oh I think that I understand know, this means that the problem comes after the solution when taking the Laplacian of the interpolating function. Now I have an additional question. What if I use such an interpolating function into the NDSolve as a source term (not its derivative) for a secondary PDE? Will this give a reliable solution or similar problems will emerge?
– DK13
Jul 25 '18 at 14:11
• The problem is two fold. The accuracy of the FEM plus the hard derivative recovery. And from the above test we know that FEM performs OK for the specific exact solution. Using interpolating functions:The higher the derivative of that interpolating function you take the less accurate it will be. Jul 25 '18 at 14:18
• in your specific example I added linearly an interpolating function in the source term of the PDE and I reproduced your procedure. The error of the manufactured solution method remains small. Can I trust this in order to conclude that the final solution is acceptable? Why near the boundary the error oscillates so much?
– DK13
Jul 25 '18 at 14:51