In FEM classes, it's usually taken for granted that the stiffness matrix is positive definite, but I just can't understand why. Could anyone give some explanation?

For instance, we can consider the Poisson problem: $$ -\nabla^2 u = f,$$ whose stiffness matrix is: $$K_{ij} = \int_\Omega\nabla\varphi_i\cdot\nabla\varphi_j\, d\Omega,$$ which is symmetric and positive definite. Symmetry is an obvious property, but the positive definiteness is not so explicite to me.

  • 1
    $\begingroup$ This actually depends on the partial differential equation you are trying to solve. Can you add the one you are interested in? $\endgroup$ – Christian Clason Nov 28 '15 at 9:09
  • $\begingroup$ Hi, @ChristianClason, thank you for your comment. I have added a concrete example of this problem. $\endgroup$ – user123 Nov 28 '15 at 9:33
  • 4
    $\begingroup$ Caveat: Without boundary conditions, the complete system stiffness matrix, as assembled from element matrices, does not have full rank, as it has to map the equivalent of rigid body motions to zero forces. Thus the complete stiffness matrix can at best be positive semidefinite. With proper boundary conditions however, rigid body motions are disabled, and the constrained system is then nonsingular. (Otherwise one could not solve it). Therefore, to find actual positive definiteness, you have to look at the condensed matrix resulting from the application of boundary conditions. $\endgroup$ – ccorn Nov 28 '15 at 12:30
  • $\begingroup$ This isn't always the case. For rigid body modal analysis with pre-stress negative eigenvalues occur (wrong ofc). $\endgroup$ – Guido Jul 21 '20 at 16:10

The property follows from the property of the corresponding (weak form of the) partial differential equation; this is one of the advantages of finite element methods compared to, e.g., finite difference methods.

To see that, first recall that the finite element method starts from the weak form of the Poisson equation (I'm assuming Dirichlet boundary conditions here): Find $u\in H^1_0(\Omega)$ such that $$ a(u,v):= \int_\Omega \nabla u\cdot \nabla v \,dx = \int_\Omega fv\,dx \qquad\text{for all }v\in H^1_0(\Omega).$$ The important property here is that $$ a(v,v) = \|\nabla v\|_{L^2}^2 \geq c \|v\|_{H^1}^2 \qquad\text{for all }v\in H^1_0(\Omega). \tag{1}$$ (This follows from Poincaré's inequality.)

Now the classical finite element approach is to replace the infinite-dimensional space $H^1_0(\Omega)$ by a finite-dimensional subspace $V_h\subset H^1_0(\Omega)$ and find $u_h\in V_h$ such that $$ a(u_h,v_h):= \int_\Omega \nabla u_h\cdot \nabla v_h \,dx = \int_\Omega fv_h\,dx \qquad\text{for all }v_h\in V_h.\tag{2}$$ The important property here is that you are using the same $a$ and a subspace $V_h\subset H^1_0(\Omega)$ (a conforming discretization); that means that you still have $$ a(v_h,v_h) \geq c \|v_h\|_{H^1}^2 >0 \qquad\text{for all }v_h\in V_h. \tag{3}$$

Now for the last step: To transform the variational form to a system of linear equations, you pick a basis $\{\varphi_1,\dots,\varphi_N\}$ of $V_h$, write $u_h =\sum_{i=1}^N u_i\varphi_i$ and insert $v_h=\varphi_j$, $1\leq j\leq N$ into $(2)$. The stiffness matrix $K$ then has the entries $K_{ij}=a(\varphi_i,\varphi_j)$ (which coincides with what you wrote).

Now take an arbitrary vector $\vec v=(v_1,\dots,v_N)^T\in \mathbb{R}^N$ and set $v_h:=\sum_{i=1}^Nv_i \varphi_i\in V_h$. Then we have by $(3)$ and the bilinearity of $a$ (i.e., you can move scalars and sums into both arguments) $$ \vec v^T K \vec v = \sum_{i=1}^N\sum_{j=1}^N v_iK_{ij} v_j = \sum_{i=1}^N\sum_{j=1}^N a(v_i\varphi_i,v_j\varphi_j) = a(v_h,v_h) >0.$$ Since $\vec v$ was arbitrary, this implies that $K$ is positive definite.

TL;DR: The stiffness matrix is positive definite because it comes from a conforming discretization of a (self-adjoint) elliptic partial differential equation.


If stiffness of element is not positive, then the system is not stable. So the model is most likely not correct. Look at the most basic equation of harmonic oscillator

$$m x''(t) + k x(t) = f(t)$$

The solution is unstable if $k$ is negative (look at the roots of the characteristic equation). It means the solution will blow up. The stiffness has to be a restoring force. At least for a physical spring. The stiffness matrix extends this to large number of elements (global stiffness matrix). That is all. But it is the same basic idea. FEM basis is in the stiffness matrix method for structural analysis where each element has a stiffness associated with it.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.