Interpolation of function onto mesh gives different results, depending on mesh density

Question

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?

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.