Introduction to Linear Algebra in the Finite Element Method

The Finite Element Method (FEM) is a numerical technique that approximates the solution of continuous physical systems governed by partial differential equations (PDEs). In structural mechanics, heat transfer, electromagnetics, and fluid dynamics, the physical behavior of a continuous medium is defined by field equations that apply over an infinite number of points. Because computers cannot solve systems with infinite degrees of freedom, we must discretize the domain.

Discretization subdivides the continuous domain into a finite number of smaller subdomains called elements. These elements are interconnected at discrete points called nodes. By assuming a simple mathematical variation (such as linear or quadratic polynomials) of the unknown field variables within each element, we can approximate the continuous differential equations as a set of algebraic equations.

This discretization process transforms the governing PDEs into a system of linear algebraic equations:

where:

Linear algebra is the mathematical language of FEM. To develop, implement, and analyze finite element formulations, we must establish a rigorous understanding of vector spaces, matrix algebra, coordinate transformations, system solvability, and eigenvalue problems.

Vectors in Finite Element Analysis

1. Geometric and Coordinate Representations

A vector $\mathbf{v}$ in an $n$-dimensional vector space $\mathbb{R}^n$ represents a quantity possessing both magnitude and direction. In engineering mechanics, vectors represent physical fields like displacements, velocities, accelerations, and forces:

Geometrically, a vector is a directed segment in space. Coordinately, it is expressed relative to a reference coordinate system. In a three-dimensional Cartesian space, a displacement vector $\mathbf{u}$ is written as:

where $\mathbf{i}, \mathbf{j}, \mathbf{k}$ are the orthogonal unit vectors along the $X$, $Y$, and $Z$ axes, respectively, and $u_x, u_y, u_z$ are the scalar components.

2. Vector Spaces, Linear Independence, and Bases

A set of vectors $V$ forms a vector space over a field (typically the real numbers $\mathbb{R}$) if it is closed under vector addition and scalar multiplication, satisfying the standard vector space axioms (commutativity, associativity, identity, distributivity).

A set of vectors $\{\mathbf{v}_1, \mathbf{v}_2, \dots, \mathbf{v}_m\}$ is linearly independent if the vector equation:

can only be satisfied by $c_1 = c_2 = \dots = c_m = 0$. If non-zero constants $c_i$ exist that satisfy this relation, the vectors are linearly dependent. In FEA, linear dependence in the displacement interpolations indicates redundant degrees of freedom, which can lead to singular stiffness matrices.

A set of linearly independent vectors that spans the vector space $V$ forms a basis for $V$. The number of vectors in this basis defines the dimension of the space. In a standard 1D linear element, the basis functions (shape functions) $N_1(x) = 1 - x/L$ and $N_2(x) = x/L$ span the space of linear polynomials over the element domain $[0, L]$.

3. Inner Products and Virtual Work

The inner product (or dot product) of two vectors $\mathbf{u}, \mathbf{v} \in \mathbb{R}^n$ is a scalar defined as:

Geometrically, the inner product is related to the angle $\theta$ between the vectors:

In computational mechanics, the inner product represents physical energy quantities. For instance, the virtual work $\delta W$ done by a virtual displacement vector $\delta \mathbf{u}$ under an external force vector $\mathbf{f}$ is the inner product:

Two vectors are orthogonal if their inner product is zero:

Orthogonality represents decoupled physical behaviors. For example, in dynamic modal analysis, the vibration mode shapes are orthogonal with respect to the mass and stiffness matrices, allowing the coupled dynamic equations to be solved as independent single-degree-of-freedom systems.

4. Outer Products

The outer product of two vectors $\mathbf{u} \in \mathbb{R}^m$ and $\mathbf{v} \in \mathbb{R}^n$ is an $m \times n$ matrix defined as:

Outer products are used to construct projection operators and to express stress and strain tensors. In 3D continuum mechanics, stress $\boldsymbol{\sigma}$ and strain $\boldsymbol{\epsilon}$ are second-order tensors, which can be represented as linear combinations of dyadic products (outer products) of basis vectors.

5. Vector Norms

A vector norm is a mapping $\|\cdot\|: \mathbb{R}^n \to \mathbb{R}$ that assigns a positive scalar to a vector, representing its length or magnitude. The norm must satisfy three properties: positive definiteness, scalability, and the triangle inequality. Common norms include:

In non-linear FEA solvers (such as the Newton-Raphson method), convergence is determined by checking if the norm of the residual force vector $\mathbf{R} = \mathbf{f} - \mathbf{K}(\mathbf{u})\mathbf{u}$ is below a specified tolerance $\epsilon$:

Matrices as Linear Transformations

A matrix $\mathbf{A} \in \mathbb{R}^{m \times n}$ represents a linear transformation that maps a vector $\mathbf{x} \in \mathbb{R}^n$ to a vector $\mathbf{y} \in \mathbb{R}^m$:

In finite element formulations, matrices map nodal displacements to nodal forces, transform coordinates between reference frames, and relate strains to stresses through constitutive laws.

1. Matrix Multiplication and Algebraic Properties

For two matrices $\mathbf{A} \in \mathbb{R}^{m \times p}$ and $\mathbf{B} \in \mathbb{R}^{p \times n}$, their product $\mathbf{C} = \mathbf{A} \mathbf{B}$ is a matrix in $\mathbb{R}^{m \times n}$ where each element is defined as:

Matrix multiplication is associative: $(\mathbf{AB})\mathbf{C} = \mathbf{A}(\mathbf{BC})$, and distributive: $\mathbf{A}(\mathbf{B} + \mathbf{C}) = \mathbf{AB} + \mathbf{AC}$, but non-commutative: $\mathbf{AB} \neq \mathbf{BA}$ in general. The identity matrix $\mathbf{I}$ satisfies $\mathbf{AI} = \mathbf{IA} = \mathbf{A}$.

2. Matrix Transpose and Symmetry

The transpose of a matrix $\mathbf{A}$, denoted $\mathbf{A}^T$, is formed by swapping its rows and columns: $(A^T)_{ij} = A_{ji}$. The transpose operation satisfies:

A square matrix $\mathbf{K} \in \mathbb{R}^{n \times n}$ is symmetric if:

In structural mechanics, the global stiffness matrix $\mathbf{K}$ is symmetric. This symmetry is a direct mathematical consequence of Betti's Reciprocal Theorem (or the Maxwell-Betti reciprocal theorem), which states that the work done by a first set of forces acting through the displacements produced by a second set is equal to the work done by the second set of forces acting through the displacements produced by the first set.

Symmetry is computationally valuable. For a system of size $n$, a symmetric matrix requires storing only the upper or lower triangular part ($n(n+1)/2$ entries instead of $n^2$), and allows the use of specialized solvers like symmetric Cholesky decomposition ($LDL^T$).

3. Positive Definiteness

A symmetric matrix $\mathbf{K} \in \mathbb{R}^{n \times n}$ is positive-definite if, for every non-zero vector $\mathbf{x} \in \mathbb{R}^n$, the quadratic form is strictly positive:

If $\mathbf{x}^T \mathbf{K} \mathbf{x} \ge 0$, the matrix is positive-semidefinite.

In structural engineering, the quadratic form $\frac{1}{2} \mathbf{u}^T \mathbf{K} \mathbf{u}$ represents the strain energy ($U$) stored in the deformed structure:

Since a physical system must store positive energy when deformed from its equilibrium position, the stiffness matrix of a properly constrained, stable structure is symmetric positive-definite (SPD). If a structure is under-constrained (i.e., it can move as a rigid body without storing internal strain energy), there exist non-zero displacement vectors $\mathbf{u}_{rb}$ (rigid body modes) such that:

In this case, the stiffness matrix is positive-semidefinite and singular.

4. Orthogonal Matrices

A square matrix $\mathbf{Q}$ is orthogonal if its transpose is equal to its inverse:

Orthogonal matrices represent rotations in space. They preserve vector lengths (Euclidean norms) and angles:

Coordinate transformations in FEM rely on orthogonal rotation matrices to transform local element quantities to the global coordinate system.

5. Sparsity and Bandedness

In a finite element mesh, a node is only physically connected to its immediate neighboring nodes within adjacent elements. Therefore, the degree of freedom at a particular node does not interact directly with distant nodes.

This local connectivity means that the global stiffness matrix $\mathbf{K}$ is highly sparse (most of its entries are zero). If the nodes are numbered sequentially, the non-zero entries are clustered around the main diagonal, forming a banded matrix.

The bandwidth $b$ of a matrix is defined as the maximum distance between a non-zero element and the diagonal:

The total width of the band is $2b + 1$. Minimizing the bandwidth is critical for computational efficiency, as the cost of solving equations is proportional to $n \cdot b^2$ for banded direct solvers. Bandwidth minimization algorithms, such as the Cuthill-McKee and Reverse Cuthill-McKee (RCM) algorithms, reorder the node numbers to cluster non-zero terms as close to the diagonal as possible.

Commercial solvers utilize advanced sparse storage schemes to avoid storing zero values, such as:

Coordinate Transformations

In finite element analysis, elements are oriented in arbitrary directions. For example, the members of a space truss or frame structure point in different directions in 3D space.

It is easiest to formulate the stiffness equations of an individual element in a local coordinate system ($x', y', z'$) aligned with the element's primary axis. However, to assemble these individual element equations into a single global system, we must express all element matrices in a common global coordinate system ($X, Y, Z$).

1. Direction Cosines and Rotation Matrices

Let the local coordinate axis $x'$ be defined by a unit vector $\hat{\mathbf{e}}_1'$ and the global Cartesian coordinate axes $X, Y$ be defined by $\mathbf{i}, \mathbf{j}$. The direction cosines are the cosines of the angles between the local and global axes:

For a 2D vector, the transformation from global components $\mathbf{u} = [u, v]^T$ to local components $\mathbf{u}' = [u', v']^T$ is written as:

where $\mathbf{T}$ is the coordinate transformation (rotation) matrix:

Since $\mathbf{T}$ is orthogonal, $\mathbf{T}^{-1} = \mathbf{T}^T$, which means the inverse transformation from local to global coordinates is:

2. Vector Transformation

For an element with multiple nodes, we construct a block-diagonal transformation matrix. For a 2D truss element with two nodes (each node having 2 DOFs: $u$ and $v$), the relation between the local displacement vector $\mathbf{u}' = [u_1', v_1', u_2', v_2']^T$ and the global displacement vector $\mathbf{u} = [u_1, v_1, u_2, v_2]^T$ is:

where $\mathbf{T}_e$ is the $4 \times 4$ element transformation matrix.

3. Matrix Transformation (Similarity Transformation)

To find how the stiffness matrix transforms, we use the principle of energy conservation. The strain energy $U$ stored in the element must be the same regardless of which coordinate system is used to evaluate it:

where $\mathbf{k}'$ is the local element stiffness matrix and $\mathbf{k}$ is the global element stiffness matrix.

Substituting the displacement transformation $\mathbf{u}' = \mathbf{T}_e \mathbf{u}$ into the energy expression:

Equating this to the global energy expression $\frac{1}{2} \mathbf{u}^T \mathbf{k} \mathbf{u}$ yields the matrix transformation formula:

This similarity transformation maps the element stiffness from local coordinates to global coordinates, allowing for global system assembly.

Systems of Linear Equations

Once the element matrices are transformed and assembled, we obtain the global system equations $\mathbf{K} \mathbf{u} = \mathbf{f}$. To solve this system, we must evaluate its algebraic properties.

1. Solvability, Determinants, and Rank

According to the Rouche-Capelli theorem, a system of linear equations $\mathbf{K} \mathbf{u} = \mathbf{f}$ with $n$ unknowns has a unique solution if and only if the rank of the stiffness matrix $\mathbf{K}$ is equal to the rank of the augmented matrix $[\mathbf{K} | \mathbf{f}]$ and equals $n$:

For a square stiffness matrix, a unique solution exists if the determinant of $\mathbf{K}$ is non-zero:

If $\det(\mathbf{K}) = 0$, the matrix is singular, and its rank is less than $n$. In this case, the system has either no solution or infinitely many solutions.

2. Physical Singularity (Rigid-Body Modes)

In structural mechanics, a singular stiffness matrix indicates that the structure is not sufficiently constrained in space. If we apply an external force to an unconstrained structure (such as a beam floating in space), the structure will translate or rotate indefinitely without producing internal strain energy.

Mathematically, these unconstrained rigid body motions correspond to the null space (or kernel) of the stiffness matrix:

Any displacement vector $\mathbf{u}_{rb} \in \text{null}(\mathbf{K})$ represents a rigid body mode. Because $\mathbf{K} \mathbf{u}_{rb} = \mathbf{0}$, the determinant of $\mathbf{K}$ is zero, making the matrix singular and impossible to invert.

To resolve this singularity and solve the system, we must apply boundary conditions (constraints) that eliminate all rigid-body modes. Applying boundary conditions reduces the system size by removing degrees of freedom with known displacements, rendering the remaining active stiffness matrix non-singular and invertible.

3. Ill-Conditioned Systems and the Condition Number

In practice, even if a system is theoretically invertible, it can be difficult to solve accurately on a computer due to floating-point precision limits. Such systems are called ill-conditioned.

The sensitivity of the solution $\mathbf{u}$ to small perturbations in the load vector $\mathbf{f}$ or the stiffness matrix $\mathbf{K}$ is measured by the condition number $\kappa(\mathbf{K})$:

If we use the $L_2$ norm, the condition number of a symmetric positive-definite stiffness matrix is the ratio of its maximum eigenvalue to its minimum eigenvalue:

In finite element models, ill-conditioning typically occurs when:

Eigenvalues and Eigenvectors

The eigenvalue problem is central to dynamic and stability analyses in engineering mechanics.

1. Mathematical Definitions and Spectral Decomposition

For a square matrix $\mathbf{A} \in \mathbb{R}^{n \times n}$, a scalar $\lambda$ and a non-zero vector $\mathbf{\phi}$ are an eigenvalue and its corresponding eigenvector if they satisfy:

To find the eigenvalues, we rewrite the equation as:

Since $\mathbf{\phi} \neq \mathbf{0}$, the matrix $(\mathbf{A} - \lambda \mathbf{I})$ must be singular, meaning its determinant is zero:

This determinant equation yields an $n$-th degree polynomial in $\lambda$, known as the characteristic equation. The $n$ roots of this polynomial are the eigenvalues $\lambda_i$.

If $\mathbf{A}$ is a symmetric matrix, its eigenvalues are real numbers, and its eigenvectors are orthogonal. We can compile the eigenvectors as columns in a modal matrix $\mathbf{\Phi} = [\mathbf{\phi}_1, \mathbf{\phi}_2, \dots, \mathbf{\phi}_n]$ and the eigenvalues on the diagonal of a matrix $\mathbf{\Lambda} = \text{diag}(\lambda_1, \lambda_2, \dots, \lambda_n)$. This allows us to decompose $\mathbf{A}$ (known as spectral decomposition or diagonalization):

Since $\mathbf{\Phi}$ is orthogonal, $\mathbf{\Phi}^{-1} = \mathbf{\Phi}^T$, which simplifies the decomposition to:

2. Physical Applications in FEM

Eigenvalue problems are used in three primary areas of finite element analysis:

A. Structural Dynamics (Vibration Analysis)

For a system under free vibration without damping, the equations of motion are derived from Newton's second law:

Assuming a harmonic solution of the form $\mathbf{u}(t) = \mathbf{\phi} \sin(\omega t)$, we substitute it into the governing equation:

Since $\sin(\omega t) \neq 0$, this reduces to the generalized eigenvalue problem:

Here, the eigenvalues represent the squares of the natural frequencies ($\lambda_i = \omega_i^2$), and the eigenvectors $\mathbf{\phi}_i$ represent the mode shapes (the relative displacement geometries of the structure at those natural frequencies).

B. Structural Stability (Buckling Analysis)

When a structure is subjected to compressive loads, it can buckle before reaching its material yield strength. The linearized buckling problem is formulated as:

where:

C. Principal Stress and Strain Transformations

At any point in a deformed solid body, the stress state is represented by a $3 \times 3$ symmetric tensor $\boldsymbol{\sigma}$:

To find the planes on which shear stresses are zero, we solve the eigenvalue problem:

The eigenvalues $\sigma_p$ are the principal stresses (maximum, middle, and minimum normal stresses), and the eigenvectors $\mathbf{n}$ are the principal directions (the unit normal vectors of the planes on which these stresses act).