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

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?

• 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? Commented Mar 29, 2019 at 15:46
• 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. Commented Mar 29, 2019 at 15:59
• 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? Commented Mar 29, 2019 at 16:07
• 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? Commented Mar 29, 2019 at 16:09

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.