I wanted to test the numerical accuracy of my program. For that I wanted to interpolate the function $$f=I_0\exp\left(-100x^2\right)\exp(-100y^2)$$ onto a grid, defined on $$\Omega=[0,1]^2$$ by using the equation $$u=f$$ and the continuous galerkin method.
In addition I compared the result with the expected (interpolated) function. Boundary conditions are $$\vec{n}\cdot u=0$$
Now I would expect the lowest value to be $0$ ($\exp(-100)\exp(-100)\approx1.3\cdot10^{-87}$), regardless of grid size, but instead I get (for $I_0=1$)

| cell number   |  Lowest value |
| 16            |  -0.003398    |
| 64            |  -1.434e-5    |
| 256           |  -5.822e-10   |
| 1024          |  -4.559e-32   |
| 4096          |  0            |

Those values are not problematic for small values for $I_0$, but as soon as $I_0$ increases I have to use a finer mesh, else my values will go into negative values. It also does not help by refining the mesh around the origin, the values at $(1,1)$ will still be negative, unless I use the critical mesh density everywhere.
Are those numerical errors, or errors resulting from the method itself? And is there a way to have the density of the grid below the critical density at places with a low gradient (as for example at $(1,1)$), without getting "incorrect" results?

  • $\begingroup$ When you say you refine the mesh around the origin, do you mean around $(0,0)$ or $(1,1)$? Also, what space of functions are you using for your solution? If you’re just using a linear approximation, you should definitely expect to require pretty hefty refinement throughout your grid to capture things accurately. If you sample your approximate solution at a grid of $(x,y)$ values and compare to the your exact function $f$, does $(1,1)$ have the worst error? Or are they similar orders of magnitude? $\endgroup$
    – spektr
    Commented Mar 29, 2019 at 15:46
  • $\begingroup$ I am using gaussian base elements, with degree 3 for calculation. If I refine the mesh adaptively, I refine it around $(0, 0)$, to catch the bigger gradients. Usually the biggest error is at $(0, 0)$, but the area which concerns me more is towards $(1,1)$, where I should get $0$, but get negative values instead. $\endgroup$
    – arc_lupus
    Commented Mar 29, 2019 at 15:59
  • $\begingroup$ When I think about how the surface of $f$ looks, it makes sense to me that you’d have negative values at $(1,1)$. Your solution field is going to have a hard time matching the curvature of the exponential, so I would expect the solution at $(1,1)$ to undershoot the exact value. What are your thoughts? What happens if you refine around $(1,1)$ too? $\endgroup$
    – spektr
    Commented Mar 29, 2019 at 16:07
  • $\begingroup$ If I refine around $(1,1)$, the values stay below zero as long as I stay below the critical refinement level (when refining at $(1,1)$, I refine the mesh globally, else it would not make sense). Or would you suggest to refine the mesh locally at $(0,0)$ and $(1,1)$ for testing? $\endgroup$
    – arc_lupus
    Commented Mar 29, 2019 at 16:09

1 Answer 1


You assume that a property that holds for the original function $f$ is also true for the $L_2$ projection $u_h=P_h f$. In your case, the property is that if $f$ is non-negative, then $u_h$ should also be non-negative.

But this assumption is not in general correct. You already know this from taking a few terms of a Fourier or Taylor series: Just because a function is non-negative, there is really no guarantee that its (truncated) Fourier series or (truncated) Taylor series is non-negative. In fact, if you think of the Fourier series of the step ("Heaviside") function, then we know that this yields in Gibbs phenomenon which means that the maximum and minimum of the truncated Fourier series exceeds the maximum and minimum of the original function. The $L_2$ project onto a finite dimensional finite element space is really of the same kind: it has over- and undershots, and that's what you are observing.

I will add that if you interpolate the function onto the mesh, the minima and maxima of the interpolated function do not exceed the extrema of the original function. That is a property of the interpolation operation that is not shared by the projection operation.

  • $\begingroup$ So, basically, the same reason as scicomp.stackexchange.com/q/29660/20545 and (I assume) scicomp.stackexchange.com/q/30270/20545, and the only way to circumvent that problem in a proper way is to increase the amount of points (similar to resolution problems in fourier transformations). Is that correct? $\endgroup$
    – arc_lupus
    Commented Mar 30, 2019 at 10:53
  • $\begingroup$ Yes, if you increase the number of points, the problem will go away asymptotically. But for each choice of the number of points, you may still get a negative value. You can't avoid that with a linear projection, but there are nonlinear projection operations you could use. $\endgroup$ Commented Apr 1, 2019 at 18:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.