# FEM: singularity of the stiffness matrix

I'm solving the differential equation $$\left( \sigma^{2}(x) u ''(x) \right)'' = f(x), \;\;\; 0 \leqslant x \leqslant 1$$ with initial conditions $u(0) = u(1) = 0$, $u''(0) = u''(1) = 0$. Here $\sigma(x) \geqslant \sigma_{0} > 0$ is parameter. In operator form we can rewrite the differential equation as $Au = f$, where operator $A$ is positive definite.

Following FEM scheme, I reduce my problem to an optimisation problem $$J(u) = (Au,u) - 2(f,u) \to \min_{u}$$ I introduce finite elements $h_{k}(x)$ as $$v_{k}(x) = \left\{ \begin{array}{rl} 1 - \left( \frac{x-x_{k}}{h} \right)^2, & x \in [x_{k-1},x_{k+1}] \\ 0, & \text{otherwise} \end{array} \right.$$ for any $k = 1,\ldots,n-1$, where $x_{k} = hk$, $h = \frac{1}{n}$. Finite elements $v_{0}(x)$ and $v_{n}(x)$ are introduced similarly.

I try to find numericaly the vector $\alpha$ such that $u(x) = \sum_{k=0}^{n} \alpha_{k} v_{k}(x)$ solves the optimisation problem. We have $$J(u) = \sum\limits_{i=0}^{n} \sum\limits_{j=0}^{n} \alpha_{i} \alpha_{j} (Av_{i},v_{j}) - \sum\limits_{i=0}^{n} 2\alpha_{i} (v_{i},f) = \alpha^{T} V \alpha - 2\alpha^{T} b \to \min\limits_{\alpha},$$ where $b_{i} = (f,v_{i})$ and $V_{i,j} = (Av_{i},v_{j})$. After differentiation with respect to $\alpha$ I receive $$V\alpha = b,$$ but here the stiffness matrix $V$ is singular. So what I have to do? Maybe I have to choose other finite elements?

-
 Hi, Nimza, do you have a test problem that you know the exact solution? If yes, try solving $V^T V \alpha = V^T b$ first to test if your basis is correct inside the domain, if everything looks correct, then maybe it is the incorrectly posed BC makes the matrix singular. But the BC seems OK to me though. – Shuhao Cao Apr 12 '12 at 14:49

## 2 Answers

In decreasing order of likelihood

1. Incorrect basis. From you description, it appears that you have exactly two quadratic functions with support on each element. That space is not a partition of unity and is not $C^1$ (continuous first derivatives). To discretize your fourth order problem directly (instead of reducing it to a system of second order equations, for example), you will need a $C^1$ basis. Note that the $C^1$ basis should be able to exactly reproduce all linear functions.

2. Insufficient boundary conditions. This will be blatantly obvious if you compute and plot the null space.

3. Incorrect assembly. Check the map from elements to assembled ordering to confirm that it is what you expected, for example that it isn't reversing the orientation of elements.

4. Incorrect local assembly. In 1D, you can analytically compute what the element stiffness matrix looks like (perhaps for a simplified case) and check that the code reproduces it.

-
Thank you. 1. I think that I will need a $C^2$ basis because $(Au,v) = \int_{0}^{1} \sigma^{2}(x) u''(x) v''(x) dx$. Then, if I consider only functions that satisfy boundary conditions then $\ker A = \{ 0 \}$. – Nimza Apr 11 '12 at 21:11
A $C^1$ basis is enough, the integrand need not be continuous. Note that the boundary conditions on second derivatives will become a boundary integral. You can use a $C^0$ basis for direct discretization of a fourth order problem, but you will need to integrate jump terms as with discontinuous Galerkin methods for first and second order systems. It's not a bad method, but it's unnecessarily complicated in 1D because it's so easy to construct bases with any order of continuity (e.g. splines). This paper is an example of "$C^0$ DG". – Jed Brown Apr 11 '12 at 23:25
Ok. I corrected my basis: now $v_{k}(x) = \cos^2 \left(\frac{\pi}{2h}(x-x_{i}) \right)$ on $[x_{i-1},x_{i+1}]$ and $i = 1,\ldots,n-1$. Now it is $C^1$. But method still doesn't work. – Nimza Apr 12 '12 at 8:58
The $C^1$ basis should be able to reproduce linear functions, but this cannot. Once you fix that, check that integrals are being performed correctly, then check boundary conditions. – Jed Brown Apr 12 '12 at 13:30

Clearly the problem has an ODD order derivative. More specifically for larger Péclet numbers, the stiffness matrix might not maintain 'fine' shape, which creates zeros during assembly and hence gets singular or sometimes very small determinant that are noticeable by the oscillations in solution plot.

The solution to these kind of problem is the use of penalty, among other methods. More specifically this is called Petrov-Galerkin method.

Sorry for my bad English comprehension.

-