Stiffness matrix

In the finite element method for the numerical solution of elliptic partial differential equations, the stiffness matrix is a matrix that represents the system of linear equations that must be solved in order to ascertain an approximate solution to the differential equation.

The stiffness matrix for the Poisson problem

For simplicity, we will first consider the Poisson problem

$-{\mathrm{\nabla }}^{2}u=f$

on some domain Template:Math, subject to the boundary condition Template:Math on the boundary of Template:Math. To discretize this equation by the finite element method, one chooses a set of basis functions Template:Math defined on Template:Math which also vanish on the boundary. One then approximates

$u\approx u^h=u_1{\phi }_1+\cdots +u_n{\phi }_n.$

The coefficients Template:Math are determined so that the error in the approximation is orthogonal to each basis function Template:Mvar:

$\underset{x\in \mathrm{\Omega }}{\int }{\phi }_i\cdot fdx=-\underset{x\in \mathrm{\Omega }}{\int }{\phi }_i{\mathrm{\nabla }}^2{u}^{h}dx=-\sum _j\left(\underset{x\in \mathrm{\Omega }}{\int }{\phi }_i{\mathrm{\nabla }}^2{\phi }_{j}dx\right)u_j=\sum _j(\underset{x\in \mathrm{\Omega }}{\int }\mathrm{\nabla }{\phi }_i\cdot \mathrm{\nabla }{\phi }_{j}dx)u_j.$

as a consequence of the homogenous Dirichlet boundary conditions. The stiffness matrix is the Template:Mvar-element square matrix Template:Math defined by

${𝐀}_{ij}=\underset{x\in \mathrm{\Omega }}{\int }\mathrm{\nabla }{\phi }_i\cdot \mathrm{\nabla }{\phi }_{j}dx.$

By defining the vector Template:Math with components ${𝐅}_i=\underset{\mathrm{\Omega }}{\int }{\phi }_{i}fdx,$ the coefficients Template:Mvar are determined by the linear system Template:Math. The stiffness matrix is symmetric, i.e. Template:Math, so all its eigenvalues are real. Moreover, it is a strictly positive-definite matrix, so that the system Template:Math always has a unique solution. (For other problems, these nice properties will be lost.)

Note that the stiffness matrix will be different depending on the computational grid used for the domain and what type of finite element is used. For example, the stiffness matrix when piecewise quadratic finite elements are used will have more degrees of freedom than piecewise linear elements.

The stiffness matrix for other problems

Determining the stiffness matrix for other PDEs follows essentially the same procedure, but it can be complicated by the choice of boundary conditions. As a more complex example, consider the elliptic equation

$-\sum _{k,l}\frac{\partial }{\partial x_k}\left(a^{kl}\frac{\partial u}{\partial x_l}\right)=f$

where $𝐀(x)=a^{kl}(x)$ is a positive-definite matrix defined for each point Template:Mvar in the domain. We impose the Robin boundary condition

$-\sum _{k,l}{\nu }_k{a}^{kl}\frac{\partial u}{\partial x_l}=c(u-g),$

where Template:Mvar is the component of the unit outward normal vector Template:Mvar in the Template:Mvar-th direction. The system to be solved is

$\sum _j(\sum _{k,l}\underset{\mathrm{\Omega }}{\int }{a}^{kl}\frac{\partial {\phi }_i}{\partial x_k}\frac{\partial {\phi }_j}{\partial x_l}dx+\underset{\partial \mathrm{\Omega }}{\int }c{\phi }_i{\phi }_{j}ds)u_j=\underset{\mathrm{\Omega }}{\int }{\phi }_{i}fdx+\underset{\partial \mathrm{\Omega }}{\int }c{\phi }_{i}gds,$

as can be shown using an analogue of Green's identity. The coefficients Template:Mvar are still found by solving a system of linear equations, but the matrix representing the system is markedly different from that for the ordinary Poisson problem.

In general, to each scalar elliptic operator Template:Mvar of order Template:Math, there is associated a bilinear form Template:Mvar on the Sobolev space Template:Mvar, so that the weak formulation of the equation Template:Math is

$B[u,v]=(f,v)$

for all functions Template:Mvar in Template:Mvar. Then the stiffness matrix for this problem is

${𝐀}_{ij}=B[{\phi }_j,{\phi }_i].$

Practical assembly of the stiffness matrix

In order to implement the finite element method on a computer, one must first choose a set of basis functions and then compute the integrals defining the stiffness matrix. Usually, the domain Template:Math is discretized by some form of mesh generation, wherein it is divided into non-overlapping triangles or quadrilaterals, which are generally referred to as elements. The basis functions are then chosen to be polynomials of some order within each element, and continuous across element boundaries. The simplest choices are piecewise linear for triangular elements and piecewise bilinear for rectangular elements.

The element stiffness matrix Template:Math for element Template:Mvar is the matrix

${𝐀}_{ij}^{[k]}=\underset{T_k}{\int }\mathrm{\nabla }{\phi }_i\cdot \mathrm{\nabla }{\phi }_{j}dx.$

The element stiffness matrix is zero for most values of Template:Mvar and Template:Mvar, for which the corresponding basis functions are zero within Template:Mvar. The full stiffness matrix Template:Math is the sum of the element stiffness matrices. In particular, for basis functions that are only supported locally, the stiffness matrix is sparse.

For many standard choices of basis functions, i.e. piecewise linear basis functions on triangles, there are simple formulas for the element stiffness matrices. For example, for piecewise linear elements, consider a triangle with vertices Template:Math, Template:Math, Template:Math, and define the 2×3 matrix

$𝐃=\left[\begin{array}{ccc}{x}_3-x_2& x_1-x_3& x_2-x_1\\ y_3-y_2& y_1-y_3& y_2-y_1\end{array}\right].$

Then the element stiffness matrix is

${𝐀}^{[k]}=\frac{{𝐃}^{𝖳}𝐃}{4\mathrm{area}(T)}.$

When the differential equation is more complicated, say by having an inhomogeneous diffusion coefficient, the integral defining the element stiffness matrix can be evaluated by Gaussian quadrature.

The condition number of the stiffness matrix depends strongly on the quality of the numerical grid. In particular, triangles with small angles in the finite element mesh induce large eigenvalues of the stiffness matrix, degrading the solution quality.

References

• Template:Citation
• Template:Citation
• Template:Citation
• Template:Citation
• Template:Citation