# FEM: Symmetry of stiffness matrix

After my first course in the Finite element method I understand that elemental and global stiffness matrices are symmetric. What is the reason behind symmetry of stiffness matrices? Is it because symmetric property of the bilinear term in weak form?

• Your last comment is already the correct answer (although the global stiffness matrix is not necessarily symmetric if the matrix is modified after assembly to incorporate Dirichlet boundary conditions). – Christian Clason Jul 16 '15 at 13:14
• The problem $\partial_t u + \vec{v}\cdot \nabla u - \kappa \Delta u = f$ for example is not symmetric. As pointed out in the answers, a non-symmetric operator yields a non symmetric stiffness matrix. – Pablo Jeken Apr 11 '20 at 9:42

## 2 Answers

If the bilinear form is symmetric (including however you choose to impose boundary conditions), and you choose to use the same kinds of trial and test functions, then the matrix will be symmetric. You can see that by writing out the form of your problem:

$$K_{ij}=b(\phi_i,\phi_j) = b(\phi_j,\phi_i) = K_{ji}$$

• @ChristianClason, if the imposition of the boundary conditions preserves the symmetry of the bilinear form, then you're good. I did not use the word "however" to imply that any imposition of the boundary conditions would be fine, thus the word "including". – Bill Barth Jul 16 '15 at 13:24
• @ChristianClason, I used a shorthand here implying that the choice of how to implement the boundary conditions is equivalent to some modifications of the bilinear form (yes, at the discrete level). Since we don't even have a definition here of $b$ take it to be $b(\cdot,\cdot) = a(\cdot,\cdot) + {\rm BCs}$ where $a$ is the original form. I wasn't trying to be rigorous, since that could take many pages to get right in general. – Bill Barth Jul 16 '15 at 13:38

The stiffness matrix is symmetric if the operator $L$ of the PDE is self-adjoint, i.e. if you have $\langle Lf,g\rangle = \langle f,Lg\rangle$ for any pair of functions $(f,g)$ in the suitable function space, where $\langle u,v\rangle$ denotes the inner product between two functions $u,v$, for instance $\int_\Omega uv dx$ ($L_2$ inner product).

The stiffness matrix is obtained by projecting the PDE $Lf = g$ onto the finite element function basis $(\Phi_i)$, in the following form: $\forall i, \langle Lf,\Phi_i\rangle = \langle g, \Phi_i\rangle$ (Galerkin formulation), then looking for the solutions in the form $f = \sum \alpha_j \Phi_j$, by injecting the second equation into the first one, you obtain the discretized equation, in the form $K \alpha = B \beta$, where $\alpha$ is the vector of unknowns (coordinates of the solution in the function basis $(\Phi_i)$), $\beta$ is the RHS projected onto the function basis ($\beta_i = \langle g, \Phi_i\rangle$), $B$ is the mass matrix ($b_{i,j} = \langle \Phi_i,\Phi_j\rangle$) and the coefficients of the matrix $k_{i,j} = \langle L \Phi_i, \Phi_j\rangle$. Therefore, if the operator $L$ is self-adjoint, you have $\langle L \Phi_i, \Phi_j\rangle = \langle L \Phi_j, \Phi_i\rangle$ and the matrix $K$ is symmetric.

This is the case of most commonly encountered PDEs, but there are exceptions when the physics gets more complicated...