Chapter 1 What is numerical analysis?

Much of today’s science and engineering depends on large-scale calculations performed with computers. These calculations find solutions or approximate solutions to mathematical models and enable scientists and engineers to predict behaviours of interest. A prime example is the weather: PDE models are used to predict the weather based on recent observations. The accuracy of the predictions depend on many things, including the accuracy of the numerical solution to the PDEs. In this unit, we explore the theoretical basis for these numerical methods, especially their reliability and efficiency. The name given to this subject is Numerical Analysis.

Numerical analysis is half theory and half practice. We want to prove that algorithms work with rigorous mathematical analysis and to implement them. An essential part of the course will be for you to implement and use the methods yourself. You will use Python to do this.

1.1 Examples

1.1.1 Quadratic-equation formula

The following formula for the two roots of a quadratic equation is well known: \[ x_\pm=\frac{-b\pm\sqrt{b^2-4ac}}{2a} \] gives the solutions \(x\) of \[ a x^2 + b x + c=0,\qquad \text{for given $a,b,c$.} \] For a pure mathematician, this could be the end of the story: Python provides built-in routines for evaluating such expressions and hence we can input \(a,b,c\) and find the roots.

It is not always so simple: computers work with finite-precision arithmetic and only store finitely many different numbers. Any number that is too long (\(\pi\), \(\sqrt{2}\),… written in base 10) or too big (above \(10^{308}\) or too small \(10^{-323}\) on my machine) causes problems. We focus on numbers that are too long: in Python and many other computing environments, real numbers are stored to 16-significant figures. That is, any number with more than 16 digits (excluding the exponent) in base 10 is rounded (by chopping or choosing the nearest) to 16 digits. Then in Python, typing import numpy as np then np.pi, we see \[ \pi \quad \text{is replaced by $3.141592653589793$.} \] Of course, computers work in base 2 and the principle there is similar.

At first, this appears like a minor irritation as the error caused is so small relative to the size of \(\pi\). However, when performing long computations in finite-precision arithmetic, the effects can accumulate to cause a catastrophic loss of accuracy. For example, consider computing \(x_{\pm}\) in the case \(b=10^6\), \(a=10^{-3}\) and \(c=10^{-3}\): \[ x_\pm= \frac{-10^6 \pm \sqrt{10^{12}-4 \times 10^{-6}}}{2\times 10^{-3}}. \] Working to 16-significant figures, the square root evaluates to \(\pm10^6\) because \(10^{12}-4\times 10^{-6}=10^{12}\) to 16 s.f. (s.f. denotes significant figures). Hence, the computed values of the roots are \(x^c_\pm = -10^9, 0\). In fact, the exact answers are \(x_\pm = -10^{9}, -10^{-9}\) (to 16 s.f.). There are no correct digits in the \(x_+^c\). Indeed, \[ \text{the absolute error in $x_+$ is } \left|x_+-x_+^c\right|\approx 10^{-9}, \] but
\[ \text{the relative error in $x_+$ is } \frac{\left|x_+-x_+^c\right|}{\left|x_+\right|}\approx 1. \] Using the relative error, we scale the error relative to what we are trying to compute and find the error in computing \(x_+\) is unacceptably large.

How can the quadratic-equation formula be evaluated accurately in this case? Numerical analysis provides ways of improving algorithms so they are less sensitive to the effects of rounding error, without the need to change computing environment. Indeed, now that single-processor computing speed is no longer increasing dramatically, it is becoming more important to exploit good algorithms. Sixteen figures is enough to represent the answer and the algorithm can be adjusted to avoid the problematic cancellation \(-10^6+\sqrt{10^{12}+\text{neglected}}\) and find the correct answer in Python. In this case, we note that one of the roots \(x_\pm\) is evaluated accurately and the second root can be computed accurately by exploiting the identity \(x_+ x_- =c/a\) for the product of the roots.

1.1.2 Linear equations

To convince you that the previous example is not overly contrived, consider the linear system of equations: \[ \begin{bmatrix} \epsilon &1\\ 0& 1 \end{bmatrix}\mathbf{x}=\mathbf{b},\qquad \mathbf{b} := \begin{bmatrix} 1\\1 \end{bmatrix} \] for a known small number \(0<\epsilon\ll 1\). We are interested in determining \(\mathbf{x}\in \mathbb{R}^2\) and it is easy to show that \[ \mathbf{x}= \begin{bmatrix} 0\\ 1 \end{bmatrix}. \] Imagine that there has been rounding error and the vector \(\mathbf{b}\) is actually stored as \([1+\delta,1]^T\) (the \(T\) denotes transpose): \[\begin{equation} \tag{1.1} \begin{bmatrix} \epsilon &1\\ 0& 1 \end{bmatrix}\mathbf{x}=\mathbf{b},\qquad \mathbf{b}:= \begin{bmatrix} 1+\delta\\1 \end{bmatrix}. \end{equation}\] Again solving the linear system, we find \[ \mathbf{x}= \begin{bmatrix} \delta/\epsilon\\ 1 \end{bmatrix}. \] In the case that \(0<\epsilon\ll \delta\) (\(\epsilon\) is much smaller than \(\delta\)), there is a large change to the solution \(\mathbf{x}\). This system and its solution \(\mathbf{x}\) is highly sensitive to small changes in input data (as modelled by \(\delta\)). This is a simple example of an system of equations and these arise widely in mathematical modelling and are particularly hard to solve accurately using numerical methods. The perturbation represented by \(\delta\) always occurs in numerical simulations due to rounding error.

In contrast to the quadratic-equation example, the ill-conditioning here is fundamental to the underlying equations and is not a consequence of the method of solution. Numerical analysis can help identify numerically stable algorithms that are less susceptible to the effects of rounding error, but, if there is an instability in the underlying model, even a good algorithm will produce wrong answers. This is why well posedness (existence, uniqueness, and continuity of solutions with respect to parameters) is studied in modules on differential equations.