# Discrete Laplacian operator with finite element discretization

Let a function $$u \in H^1_0(\Omega)$$ defined by its values at the mesh nodes. Can we compute its Laplacian using the matrix resulting from the finite elements discretization of Laplace's equation? I know this may be done using $$\Delta u = \sum_i u_i \Delta \phi_i$$ where $$\phi_i$$ are the shape functions (well defined), and $$u_i$$ are the values on the nodes, but we can't do this for $$\mathbb{P}_1$$ elements. I don't consider finite difference techniques because I'm working on arbitrary high order approximations using the shape functions $$\phi_i$$.

Edit: $$\Delta u$$ here is a weak representation of the Laplacian of $$u$$.

Weak Laplacian

Let $$u,v\in H^2(\Omega)$$ then using the divergence theorem you have the identity

$$\int_{\Omega} v \Delta u = \int_{\partial\Omega} v \partial_{\vec{n}} u - \int_{\Omega} \nabla v \cdot \nabla u.$$

If you specify your boundary conditions properly then the boundary term vanishes. Assume you want to restrict yourself to $$u,v\in H^1(\Omega)$$, then clearly $$\Delta u$$ doesn't have to exist. However if $$\int_{\Omega} v w = -\int_{\Omega} \nabla v\cdot\nabla u$$ you can pretend that $$w\in H^1(\Omega)$$ is a weak analogue of $$\Delta u$$. Thus typically the goal is to find such a $$w$$.

The Case $$u\in\mathcal{V}\subset H^1(\Omega)$$

In your problem you have some finite-dimensional space $$\mathcal{V}\subset H^{1}(\Omega)$$, with a basis $$\{\phi_1,\ldots,\phi_n\}$$. This means that any function $$v\in\mathcal{V}$$ is uniquely represented by its coefficient vector $$a\in\mathbb{R}^n$$, and it is given as $$v(x) = \sum_{j=1}^n a^j \phi_j(x)$$.

Let $$u\in\mathcal{V}$$, we want to find a corresponding function $$\Delta_{\mathcal{V}}u\in\mathcal{V}$$ which satisfies the definition of a weak Laplacian of $$u$$ in $$\mathcal{V}$$ (and thus it is clear why we have named in $$\Delta_{\mathcal{V}}u$$):

$$$$\int_{\Omega} v\Delta_{\mathcal{V}}u = -\int_{\Omega} \nabla v \cdot \nabla u.$$$$

Since $$\Delta_{\mathcal{V}} u\in\mathcal{V}$$ and $$(\phi_i)_{i=1}^n$$ is a basis for $$\mathcal{V}$$ there exist coefficients $$c$$ such that $$\Delta_{\mathcal{V}}u = \sum_j c^j \phi_j$$. Using the identity $$\int_{\Omega} v\Delta_{\mathcal{V}}u = -\int_{\Omega} \nabla v \cdot \nabla u$$, and setting $$v=\phi_i$$ we arrive at:

$$(Mc)_i = \sum_{j=1}^n c^j \int_{\Omega} \phi_i \phi_j = \int_{\Omega} \phi_i\Delta_{\mathcal{V}} u = -\int_{\Omega} \nabla \phi_i \cdot \nabla u = -\sum_{j=1}^n b^j \int_{\Omega}\nabla \phi_i \cdot \nabla \phi_j = (Wb)_i.$$

Here the matrix (or rather the $$(0,2)$$ tensor) $$W_{ij} = -\int_{\Omega}\nabla \phi_i \cdot \nabla \phi_j$$ is the negated stiffness matrix (negated since in the literature one typically considers $$-\Delta$$ since they want a positive semi-definite operator), and $$M_{ij} = \int_{\Omega} \phi_i \phi_j$$ is the mass matrix (or $$(0,2)$$ metric tensor). Since $$(\phi_i)_{i=1}^n$$ is a basis for $$\mathcal{V}$$ then $$M$$ is non-singular, then we have $$Mc = Wb \implies c = M^{-1}Wb$$. From this we arrive at

$$\Delta_{\mathcal{V}}u(x) = \sum_{i=1}^nc^i\phi_i(x) = \sum_{i=1}^n(M^{-1}Wb)^i\phi_i(x).$$

Note that since $$u \in\mathcal{V}$$ then $$\Delta_{\mathcal{V}}u$$ coincides with what you would expect to be the weak Laplacian $$w$$ of $$u$$.

The case $$u\in H^1(\Omega)$$ but $$u\not\in\mathcal{V}$$

If we do not assume that $$u$$ is from $$\mathcal{V}$$ but $$u\in H^1(\Omega)$$ then $$\Delta_{\mathcal{V}}u$$ does not coincide with the weak Laplacian $$w$$ of $$u$$. The derivation also has to be modified as $$u$$ may not be decomposed as a linear combination of $$(\phi_i)_{i=1}^n$$. Instead by analogy to the above we can use $$\beta_i = \int_{\Omega} \nabla \phi_i \cdot \nabla u$$ and compute $$c = M^{-1}\beta$$. Then the "projected" weak Laplacian of $$u$$ on $$\mathcal{V}$$ can be defined as

$$\Delta_{\mathcal{V}}u(x) = \sum_{i=1}^n c^i\phi_i(x) = \sum_{i=1}^n (M^{-1}\beta)^i\phi_i(x).$$

Also note that if $$u\in\mathcal{V}$$ then $$\beta = Wb$$ and you recover the other formula.

The case $$u\in L^2(\Omega)$$ but $$u\not\in H^1(\Omega)$$

Another option if you do not have access to $$\nabla u$$, or it simply doesn't exist, is to project $$u$$ on $$\mathcal{V}$$ and then solve the problem w.r.t. the projection $$u_{\mathcal{V}}$$. Let us define the approximation problem $$$$u_{\mathcal{V}} = \arg\min_{v\in\mathcal{V}}\left\|v-u\right\|^2_2 = \arg\min_{v\in\mathcal{V}}\int_{\Omega}(v(x)-u(x))^2\,dx.$$$$

Using the fact that $$v\in\mathcal{V}$$ we can write $$v = \sum_{i=1}^n b^i\phi_i$$. Then we can rewrite the minimisation problem as follows:

$$\min_{v\in\mathcal{V}}\left\| v-u \right\|^2_2 = \min_{b\in\mathbb{R}^n}\left\|\sum_{i=1}^n b^i\phi_i-u\right\|^2_2.$$

Taking the derivative with respect to $$b$$ and setting to zero yields the normal equations $$$$\begin{bmatrix} \int_{\Omega} \phi_1\phi_1 & \ldots & \int_{\Omega}\phi_1\phi_n \\ \vdots & & \vdots \\ \int_{\Omega}\phi_n\phi_1 & \ldots & \int_{\Omega}\phi_n\phi_n \end{bmatrix} \begin{bmatrix} b^1 \\ \vdots \\ b^n \end{bmatrix} = \begin{bmatrix} \int_{\Omega} \phi_1 u \\ \vdots \\ \int_{\Omega} \phi_n u \end{bmatrix} \implies Mb = d.$$$$ Then you get $$b = M^{-1}d$$ and thus $$u_{\mathcal{V}}(x) = \sum_{i=1}^n b^i\phi_i = \sum_{i=1}^n (M^{-1}d)^i\phi_i$$. After that you can compute the weak Laplacian as before $$c = M^{-1}Wb = M^{-1}WM^{-1}d$$: $$$$\Delta_{\mathcal{V}}u_{\mathcal{V}} = \sum_{i=1}^n c^i\phi_i = \sum_{i=1}^n (M^{-1}WM^{-1}d)^i\phi_i.$$$$

For the case $$u\in\mathcal{V}$$ you of course get $$\Delta_{\mathcal{V}} u_{\mathcal{V}} = \Delta_{\mathcal{V}} u = w$$ (here $$w$$ is the weak Laplacian of $$u$$).

Metric Tensor

Note that $$M$$ here plays the role of a metric tensor, it allows you to convert contravariant vectors to covariant ones (i.e. it lowers indices), and similarly $$M^{-1}$$ allows you to convert covariant vectors to contravariant ones (i.e. it raises indices). The stiffness matrix takes a contravariant vector and spits out a covariant one, since it is also a $$(0,2)$$ tensor, that's why you need the extra $$M^{-1}$$.

As long as the functions $$(\phi_i)_{i=1}^n$$ are linearly independent, then $$M$$ is non-singular, and in fact positive definite since $$$$a^TMa = a^T\begin{bmatrix} \int_{\Omega} \phi_1 \sum_{j=1}^n a^j \phi_j \\ \ldots \\ \int_{\Omega}\phi_n \sum_{j=1}^n a^j \phi_j \end{bmatrix} = \int_{\Omega}\left(\sum_{i=1}^na^i\phi_i(x)\right)^2\,dx = \left\|\sum_{i=1}^n a^i\phi_i\right\|^2\geq 0.$$$$ Since the $$2$$-norm is positive definite, then $$\left\|\sum_{i=1}^n a^i\phi_i\right\|^2=0 \iff \sum_{i=1}^n a^i\phi_i=0.$$ Since $$(\phi_i)_{i=1}^n$$ are linearly independent, then $$\sum_{i=1}^n a^i\phi_i=0 \iff a = 0.$$ This means that $$M$$ is positive definite, and thus non-singular and invertible.

Regarding why this is called the metric tensor, it's because you can use it compute the inner product (which can give you angles and lengths) between two functions $$v,u\in\mathcal{V}$$ using just their representations $$a,b$$ as follows $$$$\langle v, u\rangle = \int_{\Omega} v(x)u(x)\,dx = \sum_{i=1}^n\sum_{j=1}^n a^ib^j\int_{\Omega}\phi_i\phi_j = a^TMb.$$$$

A Contravariant Vector is the Coordinate Representation of a Function

A contravariant vector here is the representation of a function from your mesh, i.e. $$b = (b^1, \ldots, b^n)$$ is the representation of $$u$$ if $$u = \sum_i b^i \phi_i$$. This means that $$b$$ are the coordinates of $$u$$ in the basis $$\phi_i$$.

A Covariant Vector is the Coordinate Representation of a Functional

A covariant vector $$\tilde{b} =(b_1, \ldots, b_n)$$ is the representation of a functional on your mesh (i.e. an object that eats a function and spits out a number), the intuitive interpretation of a functional would be a measuring device (and its evaluation for a specific function is a measurement). One possible basis for the space of functionals $$\mathcal{V}'$$ (the dual space) is $$\phi^i(f) = \int_{\Omega}\phi_i f$$. Then the covariant vector $$\tilde{b}$$ defines some functional $$\tilde{u}$$:

$$\tilde{u}(f) = \sum_j b_j \phi^j(f) = \sum_j b_j \int_{\Omega} \phi_j f.$$

Edit: Edited to address Bangerth's comment regarding $$\Delta u$$ being misleading.

• This is an elegant derivation, but it is not what the OP actually looked for. You are making the assumption that you can represent $\Delta u$ as a finite element function that can be expanded in terms of basis functions, but for a given finite element function $u$ as requested, its Laplacian can of course not be written as a sum of basis functions. What you are computing here is some kind of projection of the Laplacian of the function $u$ into the finite element space. Feb 13, 2023 at 16:24
• @WolfgangBangerth I now see how this can be misunderstood, I'll rectify this. Yes $\Delta u$ is not the Laplacian applied to u, it's a function from $\mathcal{V}$. Maybe it would have been clearer if I had written $\Delta_{\mathcal{V}} u$ so there is no confusion with the usual Laplacian. I named it $\Delta u$ simply because it is the analogue of the Laplacian of $u$ on the mesh. As far as I understood gbmreda was asking precisely for the discrete analogue of the Laplacian. Feb 13, 2023 at 16:30
• @WolfgangBangerth I updated it. I had taken "Let a function $u\in H^1(\Omega)$ defined by its values at the mesh nodes" to mean that $u\in\mathcal{V}$, since gbremda noted that $u$ is defined by the values at the mesh nodes. I believe I have covered the other cases too with the rewrite however. Feb 13, 2023 at 18:09