# integration of a laplacian

Hi, I solved for a Poisson equation with finite elements, using piecewise linear basis functions on 2d triangles.
Now, I want to evaluate the following expressions:
$$\int_\Omega \Delta u ~dx$$ and $$\int_\Omega (\Delta u)^2 ~dx$$ I want to evaluate these expressions using my approximated solution $u$ which has been computed on piecewise linear basis functions.

For the first one, I thought of using the identity $\int_\Omega div \nabla u = \int_{\partial\Omega} \nabla u . \vec{n}~ds$ and summing this expression over each triangle.However, as expected, the result is strictly 0 since the basis functions are linear.
I also tried to use a kind of jump formula (like $f'(x)=\tilde{f}'(x) + f^+-f^-$ where $\tilde{f}'$ is the derivative of the smooth part of f) but I'm stuck on how to do that for each triangle in 2D (the outer normal is likely to cancel out when computing the same formula for two adjacent triangles sharing an edge) - and I'm wondering if it is supposed to work.

For the second one, I just have no clue.

Am I forced to use higher order elements ? Any idea ? Thanks!

-
Let's say you have approximate solution to $\Delta u = f$. Why don't you integrate $\int_\Omega f$? – Vít Tuček Apr 7 '11 at 18:28
because I want to compute the residual $\int (\Delta u−f)^2$ using the approximated u. – WhitAngl Apr 7 '11 at 18:53

Since your approximate solution is piecewise linear, it is $H^1$ but not $H^2$. Therefore your calculation is impossible. You can do to things to overcome the difficulty:

• Either use higher-order elements,
• or post-process your approximate solution $u$. This means constructing a smoother $\bar u$ using some convolution by an appropriate $\phi(x/h)$, where $h$ is the typical mesh size. In one space dimension, this amounts to using splines. Then $\bar u\in H^2$ and your calculation is meaningful.
-
Thanks, although it's quite disappointing, it perfectly answers my question. – WhitAngl Apr 7 '11 at 20:04
But I thought considering the derivative in the distributional sense would allow to do the integration, like with the "jump formula" (I'm not sure the translation is correct)... – WhitAngl Apr 7 '11 at 20:09

Denis has this exactly right, if your goal is really to calculate these integrals. However, if your real goal (as you say) is to calculate the residual, then this isn't what you want to do at all.

In a weak sense, the Laplacian is a map $\Delta \colon H^1 (\Omega) \to H^{-1} (\Omega)$, so the PDE $\Delta u = f$ makes sense when $u \in H^1 (\Omega)$ and $f \in H^{-1} (\Omega)$. Denoting the approximate FEM solution by $u_h$, the residual is $f - \Delta u_h \in H^{-1} (\Omega)$, so it really makes sense to measure the residual in the $H^{-1} (\Omega)$ norm, not the $L^1 (\Omega)$ or $L^2 (\Omega)$ norm. That is, $$\lVert f - \Delta u_h \rVert _{H^{-1}(\Omega)} = \sup _{ \lVert v \rVert _{H^1 (\Omega)}= 1} \langle f - \Delta u_h , v \rangle _{H^{-1} (\Omega) \times H^1 (\Omega) }.$$

On the other hand, maybe you don't really want to measure the residual itself; you want to estimate the a posteriori error $e _h = u - u_h$. In this case, $e _h \in H^1 (\Omega)$ solves the residual equation $$\Delta e _h = \Delta (u - u_h) = f - \Delta u_h .$$ You can measure $e _h$ a number of ways, e.g., using the energy norm. Typically, of course, you can't actually solve for $e_h$ (since that would mean solving the original PDE exactly!), but you can estimate it by using a more accurate finite-element method for the residual equation (e.g., finer mesh and/or higher-order elements) than you used for $u_h$.

To learn more about these sorts of things, you should look up residual-based a posteriori error estimation (Google returns lots of hits for this phrase).

-
I'm suprised that the distributional derivative shouldn't work : if I want to integrate $\int_{-1}^y |x|'' dx$, it corresponds to $\int_{-1}^y (H(x))'$ where I can use the jump formula to say that $H' = 0 + 1*\delta_0$ (where $1$ is the jump, and $0$ the continuous derivative), where the resulting distribution is a measure which integrates to $H(y)$. This should hold in 2D as well, and exactly corresponds to my case, isn't it ? Indeed the a posteriori error would be even better. However, if I could refine the mesh to compute it again, I would directly use the refined version. Thanks for all – WhitAngl Apr 8 '11 at 12:15
oops, the jump is 2 in the example. I didn't mention that $H$ was the heaviside distribution as well. – WhitAngl Apr 8 '11 at 12:18
also, in practice (ie., on a computer), for the $H^{-1}$ norm, how can I compute the $sup$ over all $v$ with $\|v\|=1$ ? – WhitAngl Apr 8 '11 at 13:46