MATH566

MATH566 Lesson 3: Programming for numerical analysis

Overview

Types of numerical approximation and their coding
Array operations

Aspects of numerical approximation

Numerical methods constructs approximation sequences ${y_{n}}_{n \in ℕ}$ to converge to the solution of a mathematical problem $y = P (x)$ , ${lim}_{n \to \infty} y_{n} = y$ .
Numerical analysis studies the convergence behavior: conditions for convergence, order and rate of convergence.
Computation theory constructs algorithms within some model of computation and studies the behavior of these algorithms, in particular the computational effort required for each step of the algorithm

Numerical method example: Newton's method for square root

y_{1} = 1, y_{n + 1} = \frac{1}{2} (y_{n} + \frac{x}{y_{n}}),

{lim}_{n \to \infty} y_{n} = \sqrt{x} = P (x) .

f (y) = R (x, y) = x - y^{2}, r_{n} = R (x, y_{n})

L_{n} (y) = f^{'} (y_{n}) (y - y_{n}) + f (y_{n})

L_{n} (y_{n + 1}) = 0 \Rightarrow y_{n + 1} = \frac{1}{2} (y_{n} + \frac{a}{y_{n}})

Numerical analysis example: Newton's method for square root

y_{1} = 1, y_{n + 1} = \frac{1}{2} (y_{n} + \frac{a}{y_{n}})

Suppose $y_{n} = \sqrt{x} (1 + δ_{n})$ , signifying that $y_{n}$ is “ $δ_{n}$ -close” to $\sqrt{x}$ . Then, from $y_{n + 1} - y_{n} = \frac{1}{2} (x / y_{n} - y_{n})$ obtain

y_{n + 1} - y_{n} = \frac{1}{2} [\frac{x}{\sqrt{x} (1 + δ_{n})} - \sqrt{x} (1 + δ_{n})] = \frac{\sqrt{x}}{2} [\frac{1}{1 + δ_{n}} - 1 - δ_{n}] \Rightarrow

y_{n + 1} - y_{n} = \frac{\sqrt{x}}{2} [(1 - δ_{n} + δ_{n}^{2} - δ_{n}^{3} + \dots) - 1 - δ_{n}] = δ_{n} \sqrt{x} \frac{2 - δ_{n}}{2 (1 - δ_{n})} ⩽ δ_{n} \sqrt{x} .

For $\forall ε > 0$ choose $N_{ε}$ such that $| δ_{N_{ε}} | < ε / \sqrt{x}$ in which case $| y_{n + 1} - y_{n} | < ε$ , hence ${y_{n}}_{n \in ℕ}$ is a Cauchy sequence in $ℝ$ , and converges.

Computation theory example: Newton's method

Newton's method for square root extraction: $f (y) = a - y^{2} = 0$
$y_{1} = 1, y_{n + 1} = \frac{1}{2} (y_{n} + \frac{a}{y_{n}})$
- At each iteration carry out 1 addition and two multiplications. Computational effort is $𝒪 (2 n)$ FLOPs (floating point operations)

Newton's method for finding zeros of $f$ , $f (y) = 0$
$y_{n + 1} = y_{n} - \frac{f (y_{n})}{f^{'} (y_{n})} .$
- At each iteration carry out 1 evaluation of $f$ and 1 evaluation of $f^{'}$ . Computational effort is $𝒪 (n)$ $f$ -evals, and $𝒪 (n)$ $f^{'}$ -evals.

Approaches to approximation sequence construction

Additive corrections
$\frac{π}{4} = 1 - \frac{1}{3} + \frac{1}{5} - \frac{1}{7} + \frac{1}{9} - \dots .,$ $L_{n} = \sum_{k = 0}^{n} \frac{{(- 1)}^{k}}{2 k + 1}, L_{n} = L_{n - 1} + \frac{{(- 1)}^{n}}{2 n + 1} \to \frac{π}{4} .$
Multiplicative corrections
$S_{n} = (\frac{2}{1} \cdot \frac{2}{3}) \cdot (\frac{4}{3} \cdot \frac{4}{5}) \cdot (\frac{6}{5} \cdot \frac{6}{7}) \dots,$ $S_{n} = \prod_{k = 1}^{n} \frac{4 k^{2}}{4 k^{2} - 1} = S_{n - 1} \cdot \frac{4 n^{2}}{4 n^{2} - 1}, S_{n} \to \frac{π}{2} .$

Approaches to approximation sequence construction

Continued fractions
$π + 3 = 6 + \frac{1^{2}}{6 +},$ $F_{n} = b_{0} + \overset{n}{\underset{k = 1}{K}} \frac{a_{k}}{b_{k}} = b_{0} + \frac{a_{1}}{b_{1} +},$ $y_{n} = 0, y_{n - 1} = \frac{a_{n - 1}}{b_{n - 1} + y_{n}}, F_{n} = b_{0} + y_{1} .$

Array operations: Matrix addition and scaling

Definition. The sum of matrices $𝑨, 𝑩 \in ℝ^{m \times n}$

$𝑨 = [a_{i, j}], 𝑩 = [b_{i, j}]$

is the matrix $𝑪 = 𝑨 + 𝑩$ with components

$𝑪 = [c_{i, j}], c_{i, j} = a_{i, j} + b_{i, j}$

Definition. The scalar multiplication of matrix $𝑨 \in ℝ^{m \times n}$ by $α \in ℝ$ is $𝑩 = α 𝑨$

$𝑨 = [a_{i, j}], 𝑩 = [b_{i, j}] = [α a_{i, j}]$

Array operations: matrix multiplication

Definition. Consider matrices $𝑨 = [\begin{array}{lll} 𝒂_{1} & \dots & 𝒂_{n} \end{array}] \in ℝ^{m \times n}$ , and $𝑿 = [\begin{array}{lll} 𝒙_{1} & \dots & 𝒙_{p} \end{array}] \in ℝ^{n \times p} .$ The matrix product $𝑩 = 𝑨 𝑿$ is a matrix $𝑩 = [\begin{array}{lll} 𝒃_{1} & \dots & 𝒃_{p} \end{array}] \in ℝ^{m \times p}$ with column vectors given by the matrix vector products

$𝒃_{k} = 𝑨 𝒙_{k}, for k = 1, 2 \dots, p .$

A matrix-matrix product is simply a set of matrix-vector products, and hence expresses multiple linear combinations in a concise way.
The dimensions of the matrices must be compatible, the number of rows of $𝑿$ must equal the number of columns of $𝑨$ .
A matrix-vector product is a special case of a matrix-matrix product when $p = 1$ .
We often write $𝑩 = 𝑨 𝑿$ in terms of columns as
$[\begin{array}{lll} 𝒃_{1} & \dots & 𝒃_{p} \end{array}] = 𝑨 [\begin{array}{lll} 𝒙_{1} & \dots & 𝒙_{p} \end{array}] = [\begin{array}{lll} 𝑨 𝒙_{1} & \dots & 𝑨 𝒙_{p} \end{array}]$

Array operations: componentwise matrix multiplication

Definition. Consider matrices $𝑨 = [\begin{array}{l} a_{i, j} \end{array}] \in ℝ^{m \times n}$ , and $𝑿 = [x_{i, j}] \in ℝ^{n \times p} .$ The matrix product $𝑩 = 𝑨 𝑿 = [b_{i, j}]$ is a matrix $𝑩 \in ℝ^{m \times p}$ with components

$b_{i, j} = a_{i, 1} x_{1, j} + a_{i, 2} x_{2, j} + \dots + a_{i, n} x_{n, j} = \sum_{k = 1}^{n} a_{i, k} b_{k, j}$

$𝑨 𝑿 = [\begin{array}{llll} a_{1, 1} & a_{1, 2} & \dots & a_{1, n} \\ a_{2, 1} & a_{2, 2} & \dots & a_{2, n} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ a_{m, 1} & a_{m, 2} & a_{m, n} \end{array}] [\begin{array}{llll} x_{1, 1} & x_{1, 2} & \dots & x_{1, n} \\ x_{2, 1} & x_{2, 2} & \dots & x_{2, n} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ x_{m, 1} & x_{m, 2} & x_{m, n} \end{array}]$

$b_{2, 1} = a_{2, 1} x_{1, 1} + a_{2, 2} x_{2, 1} + \dots + a_{2, n} x_{n, 1}$

Array operations: Transpose of matrix columns

$𝑨 \in ℝ^{m \times n}$ contains $n$ column vectors with $m$ components each
$𝑨 = [\begin{array}{llll} 𝒂_{1} & 𝒂_{2} & \dots & 𝒂_{n} \end{array}]$
The transpose switches rows and columns
$𝑨^{T} = [\begin{array}{l} 𝒂_{1}^{T} \\ 𝒂_{2}^{T} \\ ⋮ \\ 𝒂_{n}^{T} \end{array}] \in ℝ^{n \times m}$
has $n$ rows and $m$ columns

Array operations: block matrix operations

Matrix block addition
$𝑨 = [\begin{array}{ll} 𝑩 & 𝑪 \\ 𝑪 & 𝑩 \end{array}], 𝑫 = [\begin{array}{ll} 𝑬 & 𝑭 \\ 𝑭 & 𝑬 \end{array}],$ $𝑨 + 𝑫 = [\begin{array}{ll} 𝑩 & 𝑪 \\ 𝑪 & 𝑩 \end{array}] + [\begin{array}{ll} 𝑬 & 𝑭 \\ 𝑭 & 𝑬 \end{array}] = [\begin{array}{ll} 𝑩 + 𝑬 & 𝑪 + 𝑭 \\ 𝑪 + 𝑭 & 𝑩 + 𝑬 \end{array}]$
Matrix block multiplication
$𝑨 𝑫 = [\begin{array}{ll} 𝑩 & 𝑪 \\ 𝑪 & 𝑩 \end{array}] [\begin{array}{ll} 𝑬 & 𝑭 \\ 𝑭 & 𝑬 \end{array}] = [\begin{array}{ll} 𝑩 𝑬 + 𝑪 𝑭 & 𝑩 𝑭 + 𝑪 𝑬 \\ 𝑪 𝑬 + 𝑩 𝑭 & 𝑪 𝑭 + 𝑩 𝑬 \end{array}]$
Matrix block transposition
$𝑴 = [\begin{array}{ll} 𝑼 & 𝑽 \\ 𝑿 & 𝒀 \end{array}], 𝑴^{T} = [\begin{array}{ll} 𝑼^{T} & 𝑿^{T} \\ 𝑽^{T} & 𝒀^{T} \end{array}]$