Unit 4: Numerical Linear Algebra & Multi-Variable Non-Linear Systems

1. Gaussian Elimination and Computational Complexity

Solving a linear system $A\vec{x} = \vec{b}$ with $A \in \mathbb{R}^{n \times n}$ and $\vec{b} \in \mathbb{R}^n$ is the central computational primitive of scientific computing.

Naive Gaussian Elimination reduces the augmented matrix $[A \mid \vec{b}]$ to upper triangular form $[U \mid \vec{c}]$ through row operations:

Computational Operation Count (FLOPs):

• Forward Elimination: $\sum_{k=1}^{n-1} [2(n - k)(n - k + 1)] = \frac{2}{3} n^3 + \mathcal{O}(n^2)$ floating point operations.
• Backward Substitution: $\sum_{k=1}^n [2(n - k) + 1] = n^2 + \mathcal{O}(n)$ floating point operations.

The cubic complexity $\frac{2}{3} n^3$ dominates overall execution time.

2. Numerical Instability and Pivoting Strategies

If a pivot element $a_{kk}^{(k)}$ is zero, the algorithm crashes (division by zero). More insidiously, if $|a_{kk}^{(k)}| \ll 1$, the multipliers $m_{ik} = a_{ik} / a_{kk}$ become gigantic ($|m_{ik}| \gg 1$), causing catastrophic exponential growth of round-off errors and destroying numerical stability.

Partial Pivoting Algorithm:

At stage $k$, search the active column $k$ from row $k$ down to row $n$ for the element with the maximum absolute magnitude:

Interchange rows $k$ and $p$: $R_k \leftrightarrow R_p$. This guarantees that all multipliers satisfy $|m_{ik}| \le 1$, bounding error amplification and ensuring backwards stability.

Scaled Partial Pivoting:

When row elements differ by orders of magnitude, row scaling factor $s_i = \max_{1 \le j \le n} |a_{ij}|$ is computed. The pivot row $p$ is selected to maximize the relative ratio:

3. Matrix Factorizations: LU and Cholesky

The LU Decomposition ($PA = LU$):

Gaussian elimination with partial pivoting factors a permutation of $A$ into:

where:

• $P$ is a permutation matrix recording row interchanges.
• $L$ is a unit lower-triangular matrix with $1$'s on the main diagonal and the multipliers $m_{ik}$ below the diagonal.
• $U$ is the upper-triangular matrix produced at the termination of forward elimination.

Solving $A\vec{x} = \vec{b}$ reduces to two sequential $\mathcal{O}(n^2)$ triangular solves:

1. Permute right-hand side: $\vec{b}^* = P\vec{b}$.
2. Forward solve: $L\vec{y} = \vec{b}^*$.
3. Backward solve: $U\vec{x} = \vec{y}$.

Cholesky Factorization for Symmetric Positive-Definite (SPD) Matrices:

If $A = A^T$ and $\vec{x}^T A \vec{x} > 0$ for all $\vec{x} \ne \vec{0}$, then $A$ admits a unique factorization:

where $L$ is a lower triangular matrix with strictly positive diagonal entries:

Properties: Cholesky factorization requires $\frac{1}{3} n^3$ FLOPs (half the work of standard LU) and is unconditionally numerically stable without any pivoting!

1. Matrix Splitting and General Iterative Framework

When $n$ is very large (e.g. $n = 10^5$ to $10^7$ in finite element and CFD grid models) and $A$ is sparse, direct methods ($\mathcal{O}(n^3)$) are prohibitive due to fill-in. Iterative methods generate a sequence of approximations $\vec{x}^{(k)} \to \vec{x}^*$ with sparse matrix-vector products ($\mathcal{O}(n)$ per iteration).

Split $A$ into:

where $D$ is the diagonal, $-L$ is the strictly lower triangular part, and $-U$ is the strictly upper triangular part of $A$.

The general stationary linear iterative scheme has the form:

where $T$ is the iteration matrix.

2. The Jacobi Method

The Jacobi method updates each coordinate simultaneously using values from the previous iteration $k$:

In matrix form:

3. The Gauss-Seidel Method

The Gauss-Seidel method immediately utilizes updated components $x_1^{(k+1)}, \dots, x_{i-1}^{(k+1)}$ as soon as they are computed:

In matrix form:

4. Successive Over-Relaxation (SOR)

To accelerate convergence, SOR computes a weighted average of the previous iterate and the Gauss-Seidel update with relaxation parameter $\omega \in (0, 2)$:

In matrix form:

Theorem 4.1 (Ostrowski-Reich Theorem):

If $A$ is a symmetric positive-definite matrix and $0 < \omega < 2$, then the SOR method converges for any initial guess $\vec{x}^{(0)}$.

Theorem 4.2 (Optimal Relaxation Parameter $\omega_{\text{opt}}$):

For a consistently ordered tridiagonal matrix $A$, the optimal relaxation parameter $\omega_{\text{opt}}$ that minimizes the spectral radius $\rho(T_\omega)$ is:

where $\rho(T_J)$ is the spectral radius of the Jacobi iteration matrix. The corresponding optimal spectral radius is:

1. The Spectral Radius of a Matrix

Definition 4.1:

The spectral radius $\rho(M)$ of a square matrix $M \in \mathbb{C}^{n \times n}$ is the maximum absolute value of its eigenvalues:

Theorem 4.3 (Fundamental Convergence Criterion for Linear Iterative Schemes):

The iterative scheme $\vec{x}^{(k+1)} = T \vec{x}^{(k)} + \vec{c}$ converges to the unique solution $\vec{x}^* = (I - T)^{-1}\vec{c}$ for any initial vector $\vec{x}^{(0)}$ if and only if:

Rigorous Proof: Let $\vec{e}^{(k)} = \vec{x}^{(k)} - \vec{x}^*$ denote the error vector at iteration $k$. Subtracting $\vec{x}^ = T \vec{x}^ + \vec{c}$ from the recurrence:

By induction:

For the sequence to converge for every arbitrary $\vec{e}^{(0)}$, we must have:

By the Jordan Canonical Form theorem, $T = P J P^{-1}$ where $J = \text{diag}(J_1, \dots, J_m)$. The powers are $T^k = P J^k P^{-1}$. Each Jordan block of eigenvalue $\lambda$ satisfies:

As $k \to \infty$, $J_i^k \to \mathbf{0}$ if and only if $|\lambda_i| < 1$ for all eigenvalues $\lambda_i$. Hence $\lim_{k \to \infty} T^k = \mathbf{0} \iff \max |\lambda_i| = \rho(T) < 1$. $\blacksquare$

2. Strictly Diagonally Dominant Matrices

Definition 4.2:

A matrix $A \in \mathbb{R}^{n \times n}$ is strictly diagonally dominant (SDD) if for every row $i$:

Theorem 4.4:

If $A$ is strictly diagonally dominant, then both the Jacobi and Gauss-Seidel iterative methods converge unconditionally for any initial guess $\vec{x}^{(0)}$. Furthermore, Gauss-Seidel converges at least twice as fast as Jacobi:

1. Vector Formulation of Non-Linear Systems

Consider a coupled system of $n$ nonlinear equations in $n$ unknown variables:

Definition 4.3 (The Jacobian Matrix):

The Jacobian Matrix $J(\vec{x}) \in \mathbb{R}^{n \times n}$ contains the first partial derivatives of $\vec{F}$:

2. The Multi-Dimensional Newton-Raphson Scheme

Expanding $\vec{F}(\vec{x})$ about the current vector iterate $\vec{x}^{(k)}$ using the multivariate Taylor theorem:

Setting $\vec{F}(\vec{x}) = \vec{0}$ and defining the update step $\Delta \vec{x}^{(k)} = \vec{x}^{(k+1)} - \vec{x}^{(k)}$:

Key Algorithmic Principle: One must never compute the matrix inverse $J^{-1}$ explicitly! Instead, at each iteration, solve the linear system $J \Delta \vec{x} = -\vec{F}$ using LU decomposition with partial pivoting.

Theorem 4.5 (Quadratic Convergence in $\mathbb{R}^n$):

If $\vec{F} \in C^2(\mathbb{R}^n)$, $J(\vec{x}^)$ is non-singular at the root $\vec{x}^$, and $\vec{x}^{(0)}$ is chosen within a sufficiently small ball $\|\vec{x}^{(0)} - \vec{x}^*\| \le \delta$, then the multivariate Newton sequence converges quadratically: