FEniCS demo: demo_singular-poisson.py

Question

I am looking at demo_singular-poisson.py (#15 from the Documented DOLFIN demos), which solves a Neumann boundary-value problem for the Laplace operator over the unit square. I am writing with two questions:

1) I understand that, once one (builds a mesh on the domain and) specifies a function space $V$, one goes on to consider the finite-dimensional approximation (on $V$) to the exact variational problem. This finite-dimensional approximation is equivalent to a system of linear equations of the form $Ax = b$, where $A$ is a symmetric matrix.

It is stated in demo_singular.py that, for this Neumann boundary-value problem, the matrix $A$ has a non-trivial null space and, furthermore, that this null space is one-dimensional - may I ask why (or under what conditions) that is the case?

2) Can I also ask for an elaboration of what this block of code from demo_singular.py does:

# Create vector that spans the null space
null_vec = Vector(u.vector())
V.dofmap().set(null_vec, 1.0)
null_vec *= 1.0/null_vec.norm("l2")

Answer 1

Hi,

1) In short, the Poisson problem with Neumann boundary conditions
[ -\Delta u = f \quad \text{in } \Omega , \quad \nabla u \cdot n = g \quad \text{on } \partial\Omega ]

• admits an unique solution (up to an additive constant) if the following compatibility condition is satisfied:

[ \int_\Omega f d\Omega + \int_\Partial\Omega g d\Sigma = 0 ]

• admits no solution otherwise.

In fact it is easy to verify that if u is a solution, then also u + constant is a solution.
You can find more details in any PDE book.

2) That block of code is a little odd, but it basically computes the finite element interpolant of a constant function (which is the null space of the Poisson operator).
The following code will have the same effect:

null_vec = interpolate( Constant(1), V).vector()
null_vec *= 1./null_vec.norm("l2")

Comments

@umberto: Thanks for the feedback on both points.

I just have a follow-up question about 1): I understand the need for the compatibility condition on $f$ and $g$ when one considers the exact variational problem, but it's not clear to me how this compatibility condition implies that the matrix $A$ (which arises when one focuses on a finite-dimensional approximation to the exact variational problem) has a 1-dimensional null space. Are you able to help me understand this implication?

---

For P1 elements on an uniform 1-D mesh, the Poisson finite element matrix is given by

A = spdiags([-ones(n,1), 2*ones(n,1), -ones(n,1)], [-1,0,1], n, n);
A(1,1) = A(n,n) = 1

You can easily check in Matlab that the null space of A has dimension 1, being ones(n,1) a basis of the null space.

The same hold in the general 3D case and unstructured grid.

---

Thanks for your reply.

Can you suggest a reference which discusses these properties of $A$ in detail? (In addition, I would like to learn about properties of the stiffness matrix associated with a $P_1$-finite-element approximation to the weak formulation for more general second-order elliptic operators.)

Any book suggestions will be appreciated.