Solution of Linear Algebraic Equations – Direct & Iterative Methods

Introduction

Linear algebraic equations are the backbone of numerical simulation, engineering design, and data science. A system of $n$ linear equations can be written in matrix form

$$A\mathbf{x}=\mathbf{b},$$

where $A\in {\mathbb{R}}^{n\times n}$ is the coefficient matrix, $\mathbf{x}$ the unknown vector, and $\mathbf{b}$ the right‑hand side.
The goal of this unit is to equip you with systematic procedures to obtain $\mathbf{x}$ efficiently and accurately, and to understand when each method is appropriate.

Definitions

• Coefficient matrix $A$: square matrix of constants.
• Triangular matrix: all entries below (lower) or above (upper) the main diagonal are zero.
• LU decomposition: $A=LU$ where $L$ is lower triangular and $U$ is upper triangular.
• Pivoting: swapping rows (partial) or columns (full) to avoid division by small numbers.
• Iterative method: starts with an initial guess ${\mathbf{x}}^{(0)}$ and refines it using a recurrence.
• Convergence: $\Vert {\mathbf{x}}^{(k+1)}-{\mathbf{x}}^{(k)}\Vert \to 0$.
• Condition number $\kappa (A)=\Vert A\Vert \Vert A^{-1}\Vert$: measures sensitivity of the solution to perturbations.

Direct Methods

Gaussian Elimination (with Partial Pivoting)

The classic algorithm transforms $A$ into an upper triangular matrix $U$ by successive elimination of sub‑diagonal entries. The solution is then obtained by back‑substitution.

flowchart TD
    A["Start with \(A\) and \(\mathbf{b}\)"] --> B["Forward elimination"]
    B --> C["Partial pivoting"]
    C --> D["Upper triangular matrix \(U\)"]
    D --> E["Back substitution"]
    E --> F["Solution \(\mathbf{x}\)"]

Figure 1 – Gaussian elimination process

Worked Example
Solve

$$\{\begin{array}{l}3x+2y-z=1\\ 2x-5y+3z=-2\\ -1x+4y+2z=3\end{array}$$

Forward elimination (no pivot needed):

$$\begin{array}{rl}{R}_2& \leftarrow R_2-\frac{2}{3}{R}_1\\ R_3& \leftarrow R_3+\frac{1}{3}{R}_1\end{array}$$

Resulting upper triangular system:

$$\begin{array}{rl}3x+2y-z& =1\\ 0x-\frac{19}{3}y+\frac{11}{3}z& =-\frac{8}{3}\\ 0x+\frac{13}{3}y+\frac{5}{3}z& =\frac{10}{3}\end{array}$$

Back substitution:

$$z=1,y=2,x=-1$$

Thus $\mathbf{x}=(-1,2,1{)}^T$.

LU Decomposition – Crout’s Method

Crout’s algorithm factors $A$ as $A=LU$ where $L$ has non‑zero diagonal entries and $U$ has unit diagonal.

flowchart TD
    A["Start with \(A\)"] --> B["Compute \(L\) and \(U\) (Crout)"]
    B --> C["Solve \(L\mathbf{y}=\mathbf{b}\) (forward)"]
    C --> D["Solve \(U\mathbf{x}=\mathbf{y}\) (backward)"]
    D --> E["Solution \(\mathbf{x}\)"]

Figure 2 – Crout’s LU decomposition

Worked Example
Factor

$$A=\left[\begin{array}{ccc}4& 3& 2\\ 2& 5& 1\\ 1& 0& 3\end{array}\right]$$

Crout’s algorithm yields

$$L=\left[\begin{array}{ccc}4& 0& 0\\ 2& \frac{9}{2}& 0\\ 1& -\frac{3}{4}& \frac{11}{4}\end{array}\right],U=\left[\begin{array}{ccc}1& \frac{3}{4}& \frac{1}{2}\\ 0& 1& -\frac{1}{9}\\ 0& 0& 1\end{array}\right]$$

Given $\mathbf{b}=[7,8,5{]}^T$, forward solve $L\mathbf{y}=\mathbf{b}$ to get $\mathbf{y}=[7,\frac{7}{2},\frac{5}{4}{]}^T$, then back‑solve $U\mathbf{x}=\mathbf{y}$ to obtain $\mathbf{x}=[1,2,3{]}^T$.

Cholesky Decomposition

For symmetric positive‑definite $A$, $A=L{L}^T$. It requires only half the operations of LU.

Figure 3 – Cholesky decomposition

Matrix Inversion – Gauss‑Jordan

To find $A^{-1}$, augment $A$ with the identity matrix and perform row operations until the left side becomes $I$.

flowchart TD
    A["Start with \([A|I]\)"] --> B["Row operations"]
    B --> C["Left side becomes \(I\)"]
    C --> D["Right side is \(A^{-1}\)"]

Figure 4 – Gauss‑Jordan inversion

Iterative Methods

Jacobi Method

Each component of ${\mathbf{x}}^{(k+1)}$ is computed from the previous iterate:

$$x_i^{(k+1)}=\frac{1}{a_{ii}}(b_i-\sum _{j\ne i}{a}_{ij}{x}_j^{(k)}).$$

flowchart TD
    X0["Initial guess \(\mathbf{x}^{(0)}\)"] --> X1["Compute \(\mathbf{x}^{(1)}\) using Jacobi"]
    X1 --> X2["Check convergence"]
    X2 -->|"Not converged"| X1
    X2 -->|"Converged"| X3["Solution \(\mathbf{x}\)"]

Gauss‑Seidel Method

Improves Jacobi by using newly computed components immediately:

$$x_i^{(k+1)}=\frac{1}{a_{ii}}(b_i-\sum _{j<i}{a}_{ij}{x}_j^{(k+1)}-\sum _{j>i}{a}_{ij}{x}_j^{(k)}).$$

Successive Over‑Relaxation (SOR)

Adds a relaxation factor $\omega$ (1 < $\omega$ < 2) to accelerate convergence:

$$x_i^{(k+1)}=(1-\omega )x_i^{(k)}+\frac{\omega }{a_{ii}}(b_i-\sum _{j<i}{a}_{ij}{x}_j^{(k+1)}-\sum _{j>i}{a}_{ij}{x}_j^{(k)}).$$

Conjugate Gradient (CG)

Specialized for symmetric positive‑definite $A$. Uses residuals and search directions to converge in at most $n$ steps.

flowchart TD
    R0["Start with \(\mathbf{r}^{(0)} = \mathbf{b} - A\mathbf{x}^{(0)}\)"] --> P0["Set \(\mathbf{p}^{(0)} = \mathbf{r}^{(0)}\)"]
    P0 --> C1["Compute \(\alpha_k = \frac{(\mathbf{r}^{(k)})^T\mathbf{r}^{(k)}}{(\mathbf{p}^{(k)})^T A \mathbf{p}^{(k)}}\)"]
    C1 --> C2["Update \(\mathbf{x}^{(k+1)} = \mathbf{x}^{(k)} + \alpha_k \mathbf{p}^{(k)}\)"]
    C2 --> C3["Update \(\mathbf{r}^{(k+1)} = \mathbf{r}^{(k)} - \alpha_k A \mathbf{p}^{(k)}\)"]
    C3 --> C4["Check convergence"]
    C4 -->|"Not converged"| C1
    C4 -->|"Converged"| C5["Solution \(\mathbf{x}\)"]

Error Analysis

Forward Error

$${\epsilon }_f=\frac{\Vert {\mathbf{x}}_{\text{true}}-{\mathbf{x}}_{\text{num}}\Vert }{\Vert {\mathbf{x}}_{\text{true}}\Vert }.$$

Backward Error

$${\epsilon }_b=\frac{\Vert \mathbf{b}-A{\mathbf{x}}_{\text{num}}\Vert }{\Vert \mathbf{b}\Vert }.$$

Relative Error

$${\epsilon }_{\text{rel}}=\frac{\Vert {\mathbf{x}}_{\text{true}}-{\mathbf{x}}_{\text{num}}\Vert }{\Vert {\mathbf{x}}_{\text{true}}\Vert }.$$

Condition Number

A large $\kappa (A)$ indicates that small perturbations in $\mathbf{b}$ or $A$ can cause large changes in $\mathbf{x}$.

Figure 5 – Condition number effect

Comparison of Direct vs Iterative Methods

In the real world

1. Google PageRank – Computes the rank vector $\mathbf{r}$ by solving $(I-\alpha P^T)\mathbf{r}=\mathbf{e}$, where $P$ is the web‑link matrix. The system is large and sparse; Conjugate Gradient or GMRES is used.
2. eSewa transaction settlement – Balancing debits and credits across multiple user accounts can be expressed as $A\mathbf{x}=\mathbf{b}$ where $A$ is a sparse incidence matrix. Gauss‑Seidel with partial pivoting ensures fast settlement.
3. Ncell network optimization – Allocation of bandwidth to base stations is modeled as a linear system derived from flow conservation constraints. LU decomposition with partial pivoting solves the system quickly for real‑time adjustments.

Worked Example – Heat Conduction in a 2×2 Plate

A square metal plate of side $1\text{m}$ is discretized into a $3\times 3$ grid (including boundaries). Boundary temperatures: left and right edges $100^{\circ }\text{C}$, top and bottom edges $0^{\circ }\text{C}$. The interior node temperatures $T_{11},T_{12},T_{21},T_{22}$ satisfy

$$\begin{array}{rl}4T_{11}-T_{12}-T_{21}& =200\\ -T_{11}+4T_{12}-T_{22}& =200\\ -T_{11}+4T_{21}-T_{22}& =200\\ -T_{12}-T_{21}+4T_{22}& =200\end{array}$$

Matrix form $A\mathbf{T}=\mathbf{b}$ with

$$A=\left[\begin{array}{cccc}4& -1& -1& 0\\ -1& 4& 0& -1\\ -1& 0& 4& -1\\ 0& -1& -1& 4\end{array}\right],\mathbf{b}=\left[\begin{array}{c}200\\ 200\\ 200\\ 200\end{array}\right].$$

Using Gaussian elimination (no pivot needed), we obtain

$$\mathbf{T}=\left[\begin{array}{c}100\\ 100\\ 100\\ 100\end{array}\right]{}^{\circ }\text{C}.$$

Figure 6 – Grid and matrix for heat conduction

The solution shows that the interior nodes reach the same temperature as the boundaries due to symmetry.

Exam tip

• Know the algorithmic steps of Gaussian elimination, LU, and iterative methods; be able to write them in pseudocode.
• Practice pivoting: identify when partial pivoting is required and how it changes the algorithm.
• Solve a full system by hand (or with a calculator) to reinforce forward/backward substitution.
• Understand convergence criteria for iterative methods; be prepared to compute the spectral radius or use the relaxation factor $\omega$.
• Error analysis: be able to compute forward, backward, and relative errors, and explain the impact of the condition number.
• Draw a small matrix and perform one elimination step to demonstrate your understanding of row operations.