FEM: singularity of the stiffness matrix

Question

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?

Answer 1

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.

Answer 2

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.