I am solving the Dirichlet problem $$ \begin{cases} \Delta u = 0, \\ u|_{\partial D} = f, \end{cases} $$ in a $2d$ domain $D$ using the finite element method. What I want to get is the normal derivative of the solution $u$ on the boundary: $\tfrac{\partial u}{\partial \nu}$. In other words, I'm interested in FEM implementation of the Dirichlet-to-Neumann (Poincare-Steklov) map $\Lambda$.

Hovewer, I could not find an explicit answer to this problem in Google. I'm programming in Matlab and I would extremely appreciate an optimized solution.

I tried to solve this problem by formulating it in a weak form: find $h \in L^2(\partial D)$ such that for any $w \in H^{\frac 1 2}(\partial D)$, $$ \int_{\partial D} hw \,dl = \int_D \nabla w \nabla u dS. $$ As usually, I represent $h = \sum' h_j \psi_j$, $w = \sum w_i \psi_i$, $u = \sum u_l \psi_l$, where $\psi_i$ are the basis elements and $\sum'$ denotes summation over the boundary elements only, and I get $$ \sum_i' \sum_j' h_j w_i \int_{\partial D} \psi_i \psi_j dl = \sum_k \sum_l w_k u_l \int_D \nabla \psi_k \nabla \psi_l dS. $$ Since $\{w_j\}$ are arbitrary, we get the system $$ \sum_j' \int_{\partial D} \psi_i \psi_j dl \cdot h_j = \sum_l\int_D \nabla \psi_k \nabla \psi_l dS \, \cdot u_l. $$ Introducing the matrices $$ B_{ij} = \int_{\partial D} \psi_i \psi_j \, dl, \quad A_{ij} = \int_D \nabla \psi_i \nabla \psi_j dS, $$ where in the first matrix $i$, $j$ correspond to boundary elements and in the second matrix $i$ corresponds to the boundary elements and $j$ corresponds to all elements, I get the matrix equation $$ B h = Au, $$ for finding $h$.

However, if I implement this approach, I get some nonsense. For example, consider the following circular domain and boundary values $u_0(x,y) = -y$

enter image description here

Then the flux along the boundary must coincide with $u_0(x,y)$

I am trying to compute the flux of with the equations used above. However I was wondering what functions to use for $h$? If I use piecewise constant on elements, I get that the $B$ matrix has rank 1 too small (because its a circle), but how can one know this beforehand? (I tried to adjoin this with that the total flux must be $0$ but it gives wrong results). If I try to use a continous at vertices and linear on elements function, I get a flux which misses a factor of $1.4$. If I try to use discontinous but linear elements the $B$ matrix has more columns then rows. For $w$ I used continous piecewise linear elements which are $1$ only at nodes on the boundary, is this right?

What is the precise elements to use for all the functions, and can someone write a pseudo code for this problem?

  • $\begingroup$ You'd get more answers if you inlined the problem description, rather than link to some other site. $\endgroup$ Commented Dec 6, 2017 at 14:05
  • $\begingroup$ This problem was also considered here scicomp.stackexchange.com/questions/27002/… $\endgroup$
    – badmf
    Commented Dec 6, 2017 at 14:16
  • $\begingroup$ I'm not sure your formulation is correct. Can you explain how you ended up with the weak formulation for $h$? Secondly, the test function $w$ lives everywhere, so it can't be in $H^{1/2}(\partial\Omega)$ but in fact needs to be in $H^1(\Omega)$. (Though this then implies that its restriction to the boundary is in $H^{1/2}(\partial\Omega)$.) $\endgroup$ Commented Dec 6, 2017 at 16:38
  • $\begingroup$ The weak form is from 2nd answer here scicomp.stackexchange.com/questions/26774/… $\endgroup$ Commented Dec 6, 2017 at 17:51
  • $\begingroup$ But at least the issue with the space of the test functions is wrong there already. $\endgroup$ Commented Dec 6, 2017 at 18:21

1 Answer 1


The formulation of this problem is tricky. Here is what you have in your original post:

Find $h \in L^2(\partial D)$ such that for any $w \in H^{\frac 1 2}(\partial D)$, $$ \int_{\partial D} hw \,dl = \int_D \nabla w \nabla u dS. $$ This already can't be quite right because you have $w$ on the right in a domain integral, so it cannot be a function $w \in H^{\frac 1 2}(\partial D)$ that only lives on the boundary. Rather, you need to consider all test functions $w\in H^1(D)$ for this.

Now, this means that on the left you have that the functions $h,w$ are from different spaces: you have "many more" test functions $w$ than solution functions $h$. If you discretized this, you'd get a rectangular matrix with more rows than columns.

Of course, in practice the "vast majority" of $H^1$ functions $w$ are in fact zero on the boundary. Indeed, the set of test functions $w\in H^1$ whose traces $w|_{\partial D}$ are linear independent (i.e.,for which the corresponding rows of the matrix on the left will be linearly independent) can all be written as the extension of $H^{1/2}(\partial D)$ functions into the interior. How exactly one extends into the interior does not actually matter; let us just say that $Ew \in H^1(D)$ is an extension of $w\in H^{1/2}(\partial D)$ so that $Ew|_{\partial D}=w$.

Then, the now correct problem formulation is as follows: Find $h \in H^{-1/2}(\partial D)$ such that for any $w \in H^{\frac 1 2}(\partial D)$, $$ \int_{\partial D} hw \,dl = \int_D \nabla (Ew) \cdot \nabla u dS. $$

This is now the point where you can discretize. To this end, let $V_h \subset H^{1/2}(\partial D)$ be a finite dimensional subspace, for example consisting of piecewise constant functions on the boundary relative to the faces of finite element mesh of $D$. Then, if $\{\psi_j\}_{j=1}^N$ is a basis of $V_h$, you would get the following linear system for the finite element approximation $h_h = \sum_j H_j \psi_j$: $$ \sum_j \left\{ \int_{\partial D} \psi_i \psi_j dl \right\} H_j = \int_D \nabla (E\psi_i) \nabla u_h dS. $$ But, importantly, the finite element approximation $u_h$ is defined with respect to a different set of shape functions, say $\varphi_k \in W_h \subset H^1(D)$ so that $u_h = \sum_k U_k \varphi_k$. It is important to realize that the $\varphi_k$ are not just extensions of the $\psi_k$ -- their number is entirely different, and they form an unrelated space. With this, you can then write the problem as you wanted to: $$ \sum_j \left\{ \int_{\partial D} \psi_i \psi_j dl \right\} H_j = \sum_k \left\{\int_D \nabla (E\psi_k) \nabla \varphi_k dS\right\} U_k. $$ If you want to write this as a linear system with matrices, you'd get this: $$ B_{ij} = \int_{\partial D} \psi_i \psi_j \, dl, \quad A_{ij} = \int_D \nabla (E\psi_i) \nabla \varphi_j dS, $$ and the linear system to find the coefficients $H_j$ of $h_h$ is given by $$ B H = A U. $$ Here, $B$ is simply a mass matrix on the boundary, and consequently symmetric and positive definite. $A$, on the other hand, is a rectangular matrix.

Finally, we had not chosen an extension operator. For practical reasons, it would of course be nice if $Ew_i$ is a finite element function, for example one in $W_h$ (or at least a polynomial on each cell) because then you can just apply quadrature as always. This is easily achieved by making sure that $Ew_i$ is zero at all interior nodes. If you do this, then $A_{ij}$ is zero for all $j$ corresponding to nodes that are neither at the boundary or on cells that are adjacent to the boundary. (This would suggest more efficient ways to store the matrix $A$, but that's besides the point for the current post.)

  • $\begingroup$ What is this $H^{-1/2}$ space? $\endgroup$
    – badmf
    Commented Dec 11, 2017 at 12:00
  • $\begingroup$ @WolfgangBangerth that is the great answer. Correct me if I am wrong, you restricting yourself to piecewise continuous in the domain and piecewise discontinuous functions on the boundary. You are taking advantage of that because $V_h$ has dual pace $V^*_h$, and space $V_h$ is pivot space, so those are the same and thus $V^*_h \subset H^{-1/2}$. As a result you have simple mass matrix. $\endgroup$
    – likask
    Commented Dec 11, 2017 at 17:02
  • 1
    $\begingroup$ @JustMe: $H^{-1/2}$ is the space of all functions that are Neumann values $\partial u/\partial n$ of $H^1$ functions $u$ restricted to the boundary. $\endgroup$ Commented Dec 11, 2017 at 18:30
  • $\begingroup$ @likask: That's one way to see it. But also because $H^{1/2} \subset L_2 \subset H^{-1/2}$, you automatically have that $V_h \subset H^{-1/2}$ for any of the common finite element spaces you may choose. $\endgroup$ Commented Dec 11, 2017 at 18:32

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.