Conjugate Gradient for singular 2D poisson finite element with Neumann Boundary Conditions

Question

Heavily edited question after I realised partly what the problem was

I have programmed a simple 2D square finite element solution to the Poisson equation

$-\Delta u = f$

The source function integrates to zero and I'm implementing as natural boundary conditions. This leads to a tridiagonal block matrix structure. As follows: $$A = \begin{bmatrix} D_e & T_1 & & \\\ T_1 & D & T_1 & \\\ & \ddots & \ddots & \ddots \\\ & & T_1 & D_e \end{bmatrix} $$ $A$ is symmetric, positive semi-definite, but singular (invariant under transformations of $u+c$ Since $\int_\Omega f dx = 0$ I don't do anything on the right hand side to implement $\int_{\partial \Omega} g dx = 0$.

Of course the system is not uniquely solvable (up to a constant) so I included a row of 1s on the bottom of the sparse matrix to represent $\int_\Omega u dx = 0$. I also add a column of ones on the right side to represent the lagrangian to solve for (?). The bottom right corner is a zero.

This answer has some information about how my system is setup.

Anyway it is no problem to solve the sparse system using scipy's library. The result is as follows:

Now see the result when I try to solve with my own CG method or that of scipy, similar result:

The blocks clearly have a drop in value at their centres.

Anybody have experience similar to this? Any ideas what I've done wrong? I'm using np.linalg.cond and have found that the condition number is 626.

The error sometimes jumps upwards which seems incorrect.

Thanks in advance

Answer 1

I took a look at the eigenvalues produced by scipy.sparse.linalg.eigs and np.linalg.eig and noticed that in the process of modifying the stiffness matrix, the 0 eigenvalue is removed but two new eigenvalues appear at $\pm \sqrt n$ where n is the size of $A \in \mathbb{R}^{n\times n}$.

So the matrix now has a negative eigenvalue. Solution to this was to use a scipy.sparse.linalg.bicgstab gradient method which produced the following and is better but not really great:

Solution was to code my own Biconjugate Gradient Stabilised which led to this:

for a different (randomly generated) source: