Iterative Methods - Executable Engineering: Interactive Numerical Methods in Python

import numpy as np
import matplotlib.pyplot as plt

try:
    import executable_engineering as exe
except ImportError:
    %pip install -q executable_engineering
    import executable_engineering as exe

Direct vs Iterative Solvers¶

Gaussian Elimination and $\mathbf{L}\mathbf{U}$ decomposition are direct solvers in that they solve systems in a fixed number of operations. If we could calculate in infinite precision, these answers would be exact. Direct solvers generally scale $\approx \mathcal{O}(n^3)$ and are therefore impractical for massive systems.

We will now talk about iterative / indirect methods which start with an initial guess and iterate (hopefully) towards a solution. Each iteration is typically faster and more memory efficient than direct solvers, but convergence is only guaranteed under certain limited circumstances. Even then, the rate of convergence can still make these techniques impractical.

Properties of Iterative Solvers¶

Iterative techniques have some very useful properties:

They are self-correcting in that roundoff (or arithmetic) errors are corrected through further cycles.
They are vastly more conducive to sparse matrix storage.

Only the solution vector $\mathbf{x}$ is altered during these computations. The coefficient matrix $\mathbf{A}$ is not altered during any of these steps!

Matrix-Free Solvers: In fact, there is a class of solvers called matrix free solvers that avoid explicitly writing out $\mathbf{A}$ entirely! This is useful if the matrix-vector product $\mathbf{A}\mathbf{x}$ can be calculated in a more effective manner on-the-fly, saving massive amounts of memory.

Quantifying Convergence¶

Each iteration will (hopefully) improve the answer asymptotically, so we estimate the error and stop the algorithm when the answer stops getting better. As discussed in error analysis, this involves defining tolerances for computation.

The system has converged when $\mathbf{A}\mathbf{x}=\mathbf{b}$ . We define the residual vector, $\mathbf{r} = \mathbf{A}\mathbf{x}-\mathbf{b}$ , such that when $\mathbf{r}\rightarrow 0$ , the system is converged.

Let’s use the norm (some measure of magnitude) of $\mathbf{r}$ and say the system has converged when $||\mathbf{r}|| < \text{tolerance}$ .

Questions for thought:
Should we use absolute tolerance or relative tolerance?
What other tolerance could we imagine using? (Hint: We are actually interested in $\mathbf{x}$ ...)