# 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

$$K_{ij}=b(\phi_i,\phi_j) = b(\phi_j,\phi_i) = K_{ji}$$
• @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.