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: Solution using scipy'S spsolve

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

CG approximation of solution 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

  • $\begingroup$ Before I accept my own answer, I will wait to see if someone else can shine light on this... $\endgroup$ Oct 4, 2019 at 20:41

1 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:

enter image description here

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

Soln with self programmed BiCGstab

for a different (randomly generated) source:

Source map for previous soln


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.