Direct Versus Iterative Solution Methods

This section describes the direct and iterative solution methods used to solve the linear system of equations obtained after spatial and temporal discretization of the governing equations.

Most of the methods of discretization discussed (FDM, FVM and FEM) yield a linear system of equations that need to be solved.

The resultant expression is of the form: (1)

A x = b

that is (2)

[\begin{matrix} \begin{matrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{matrix} & \dots & \begin{matrix} a_{1 n} \\ a_{2 n} \end{matrix} \\ ⋮ & ⋱ & ⋮ \\ \begin{matrix} a_{n 1} & a_{n 2} \end{matrix} & \dots & a_{n n} \end{matrix}] [\begin{matrix} \begin{matrix} x_{1} \\ x_{2} \end{matrix} \\ ⋮ \\ x_{3} \end{matrix}] = [\begin{matrix} \begin{matrix} b_{1} \\ b_{2} \end{matrix} \\ ⋮ \\ b_{3} \end{matrix}]

where

$A$ is a matrix of known coefficients.
$b$ is a vector of given coefficients.
$x$ is a vector of unknowns such as density, velocity and temperature.

The methods used to solve the linear system of equations can be classified into two categories:

Direct Methods

Direct methods obtain an exact solution of the linear system of equations. These methods generally consist of elementary matrix operations to simplify the equations that will allow an exact computation. Some of these methods are:

Cramer’s Rule: The solution to $A x = b$ is computed directly by $x = A^{- 1} b$ and using Cramer’s rule to compute the inverse of the coefficient matrix by calculating the co factor and adjoint of matrix $A$ .
LU Decomposition: This approach is based on the assumption that matrix $A$ is composed of an upper triangular and a lower triangular matrix, that is $A = L U$ . This results in $A x = (L U) x = L (U x) = L y$ after which the systems $L y = b$ and $U x = y$ are solved separately.
Gaussian Elimination: The given linear system of equations is converted to an augmented form given by $[\begin{matrix} \begin{matrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{matrix} & \dots & \begin{matrix} a_{1 n} \\ a_{2 n} \end{matrix} \\ ⋮ & ⋱ & ⋮ \\ \begin{matrix} a_{n 1} & a_{n 2} \end{matrix} & \dots & a_{n n} \end{matrix} | \begin{matrix} \begin{matrix} b_{1} \\ b_{2} \end{matrix} \\ ⋮ \\ b_{3} \end{matrix}] [\begin{matrix} \begin{matrix} x_{1} \\ x_{2} \end{matrix} \\ ⋮ \\ x_{3} \end{matrix}]$ which is then converted to an upper triangular form after elementary operations. The solution to the original system of equations is found by back substitution from the last row of the augmented matrix.

Direct methods work well for simple problems which generate lower order matrices but they are computationally expensive and inefficient for large sparse matrices. Moreover the coefficient matrix $A$ depends on the properties of the fluid, such as viscosity or density which evolves as the solution is refined for the primary variables. Therefore it is not efficient to seek an exact solution to $A x = b$ for every global iteration because the variables used to build $A$ and $b$ are evolving.

Iterative Methods

Iterative methods obtain an approximate solution of the linear system of equations that is based on iterations. An approximate solution of the primary variables is calculated which is then fed back into the primary equations to refine the variables. For large scale CFD problems it is more efficient to slightly improve over each variable at one turn and then cycle over to all the coupled variables in the subsequent iteration to reduce the error.

Iterative methods generate an approximate solution to the system that tends to converge to an exact solution. After

m

iterations the approximation is described as: (3)

A x^{(m)} = b - r^{(m)}

where

r^{(m)}

is the residual after

m

iterations. The error in the approximation is defined as: (4)

ϵ^{(m)} = x - x^{(m)}

which results in the expression (5)

A ϵ^{(m)} = r^{(m)}

The purpose of the iterative methods is to drive this residual below the convergence criteria set.

The solution of the linear system of equations can be expressed as: (6)

x^{(m + 1)} = C x^{(m)} + d

Based on the way that the solution at every iteration is updated, iterative methods can be classified into two types:

Stationary Iterative Methods: The terms $C$ and $d$ are not dependent on the iteration count $m$ . Some examples of these methods include Jacobi method, Gauss-Seidel method, Successive Overrelaxation method (SOR) and Symmetric Successive Overrelaxation method (SSOR).
Non-Stationary Iterative Methods: The terms $C$ and $d$ are updated at every iteration. Some examples of these methods include Conjugate Gradient method (CG), Generalized Minimal Residual method (GMRES) and Conjugate Gradient Squared Method (CGS).