MATH661

1.Numbers

1.1.Number sets

Most scientific disciplines introduce an idea of the amount of some entity or property of interest. Furthermore, the amount is usually combined with the concept of a number, an abstraction of the observation that the two sets $A =$ {Mary, Jane, Tom} and $B =$ {apple, plum, cherry} seem quite different, but we can match one distinct person to one distinct fruit as in {Mary $\leftrightarrow$ plum, Jane $\leftrightarrow$ apple, Tom $\leftrightarrow$ cherry}. In contrast, we cannot do the same matching of distinct persons to a distinct color from the set {red, green}, and one of the colors must be shared between two persons. Formal definition of the concept of a number from the above observations is surprisingly difficult since it would be self-referential due to the apperance of the numbers “one” and “two”. Leaving this aside, the key concept is that of quantity of some property of interest that is expressed through a number.

Several types of numbers have been introduced in mathematics to express different types of quantities, and the following will be used throughout this text:

$ℕ$: The set of natural numbers, $ℕ = {0, 1, 2, 3, \dots}$ , infinite and countable, $ℕ_{+} = {1, 2, 3, \dots}$ ;
$ℤ$: The set of integers, $ℤ = {0, \pm 1, \pm 2, \pm 3, \dots}$ , infinite and countable;
$ℚ$: The set of rational numbers $ℚ = {p / q, p \in ℤ, q \in ℕ_{+}}$ , infinite and countable;
$ℝ$: The set of real numbers, infinite, not countable, can be ordered;
$ℂ$: The set of complex numbers, $ℂ = {x + i y, x, y \in ℝ}$ , infinite, not countable, cannot be ordered.

These sets of numbers form a hierarchy, with $ℕ \subset ℤ \subset ℚ \subset ℝ \subset ℂ$ . The size of a set of numbers is an important aspect of its utility in describing natural phenomena. The set $S = {Mary, Jane, Tom}$ has three elements, and its size is defined by the cardinal number, $| S | = 3$ . The sets $ℕ, ℤ, ℚ, ℝ, ℂ$ have an infinite number of elements, but the relation

z = {\begin{cases} - n / 2 & for n even \\ (n + 1) / 2 & for n odd \end{cases} .

defines a one-to-one correspondence between $n \in ℕ$ and $z \in ℤ$ , so these sets are of the same size denoted by the transfinite number $ℵ_{0}$ (aleph-zero). The rationals can also be placed into a one-to-one correspondence with $ℕ$ , hence

| ℕ | = | ℤ | = | ℚ | = ℵ_{0} .

In contrast there is no one-to-one mapping of the reals to the naturals, and the cardinality of the reals is $| ℝ | = 𝔠$ (Fraktur-script c). Georg Cantor established set theory and introduced a proof technique known as the diagonal argument to show that $𝔠 = 2^{ℵ_{0}}$ . Intuitively, there are exponentially more reals than naturals.

1.2.Quantification

One of the foundations of the scientific method is quantification, ascribing numbers to phenomena of interest. To exemplify the utility of different types of number to describe natural phenomena, consider common salt (sodium chloride, Fig. 1) which has the chemical formula NaCl with the sodium ions ( ${Na}^{+}$ ) and chloride ions ( ${Cl}^{-}$ ) spatially organized in a cubic lattice, with an edge length $a = 5.6402$ Å (1 Å = $10^{- 10}$ m) between atoms of the same type. Setting the origin of a Cartesian coordinate system $O x y z$ at a sodium atom, the position of some atom within the lattice is

(x, y, z) = (i \frac{a}{2}, j \frac{a}{2}, k \frac{a}{2}) .

Sodium atoms are found positions where $i + j + k$ is even, while chloride atoms are found at positions where $i + j + k$ is odd. The Cartesian coordinates $(x, y, z)$ describe some arbitrary position in space, which is conceptualized as a continuum and placed into one-to-one correspondence with $ℝ^{3}$ . A particular lattice position can be specified simply through a label consisting of three integers $(i, j, k) \in ℤ^{3}$ . The position can be recovered through a scaling operation

(x, y, z) = \frac{a}{2} (i, j, k),

and the number $a / 2 \in ℝ$ that modifies the length scale from $1$ to $a / 2$ , it is called a scalar.

Figure 1. Left: Polycrystalline sodium chloride. Right: Cubic lattice structure of a single sodium chloride crystal

1.3.Computer number sets

A computer has a finite amount of memory, hence cannot represent all numbers, but rather subsets of the above number sets. Current digital computers internally use numbers represented through binary digits, or bits. Many computer number types are defined for specific purposes, and are often encountered in applications such as image representation or digital data acquisition. Here are the main types.

Subsets of $ℕ$

The number types uint8, uint16, uint32, uint64 represent subsets of the natural numbers (unsigned integers) using 8, 16, 32, 64 bits respectively. An unsigned integer with $b$ bits can store a natural number in the range from $0$ to $2^{b} - 1$ . Two arbitrary natural numbers, written as $\forall i, j \in ℕ$ can be added and will give another natural number, $k = i + j \in ℕ$ . In contrast, addition of computer unsigned integers is only defined within the specific range $0$ to $2^{b} - 1$ . If $k > 2^{b} - 1$ , the result might be displayed as the maximum possible value or as $k mod 2^{b}$ .

∴	i=UInt8(15); j=UInt8(10); k=i+j

$25$

∴	i=UInt8(150); j=UInt8(200); k=i+j; [k i+j mod(350,256)]

$[\begin{array}{ccc} 94 & 94 & 94 \end{array}]$ (1)

∴	i=UInt8(150); j=UInt8(200); k=i-j; [k i-j mod(-50,256)]

$[\begin{array}{ccc} 206 & 206 & 206 \end{array}]$ (2)

∴	typeof(i-j)

UInt8

∴

Subsets of $ℤ$

The number types int8, int16, int32, int64 represent subsets of the integers. One bit is used to store the sign of the number, so the subset of $ℤ$ that can be represented is from $1 - 2^{b - 1}$ to $2^{b - 1} - 1$ .

∴	i=Int8(15); j=Int8(21); k=i+j

$36$

∴	i=Int8(100); j=Int8(101); k=i+j; [k i+j mod(201,128)-128]

$[\begin{array}{ccc} - 55 & - 55 & - 55 \end{array}]$ (3)

∴

typeof(k)

Int8

∴	[typemin(Int8) typemax(Int8)]

$[\begin{array}{cc} - 128 & 127 \end{array}]$ (4)

∴

Subsets of $ℚ, ℝ, ℂ$

Computers approximate the real numbers through the set $𝔽$ of floating point numbers. Floating point numbers that use $b = 32$ bits are known as single precision, while those that use $b = 64$ are double precision. A floating point number $x \in 𝔽$ is stored internally as $x = \pm . B_{1} B_{2} \dots B_{m} \times 2^{\pm b_{1} b_{2} \dots b_{e}}$ where $B_{i}$ , $i = 1, \dots, m$ are bits within the mantissa of length $m$ , and $b_{j}$ , $j = 1, \dots, e$ are bits within the exponent, along with signs $\pm$ for each. The default number type is usually double precision, more concisely referred to Float64. Common irrational constants such as $e$ , $π$ are predefined as irrationals and casting to Float64 or Float32 gives floating point approximation. Unicode notation is recognized. Specification of a decimal point indicates a floating point number; its absence indicates an integer.

∴

pi

$π$

∴	typeof(pi)

Irrational{:π}

∴	[Float32(pi) Float64(pi) Float64(π)]

$[\begin{array}{ccc} 3.1415927410125732 & 3.141592653589793 & 3.141592653589793 \end{array}]$ (5)

∴	a=2.3; b=2; c=3.; [typeof(a) typeof(b) typeof(c)]

DataType[Float64 Int64 Float64]

∴

The approximation of the reals $ℝ$ by the floats $𝔽$ is characterized by: floatmax(), the largest float, floatmin the smallest positive float, and eps() known as machine epsilon. Machine epsilon highlights the differences between floating point and real numbers since it is defined as the largest number $ϵ \in 𝔽$ that satisfies $1 + ϵ = 1$ . If $ε \in ℝ$ of course $1 + ε = 1$ implies $ε = 0$ , but floating points exhibit “granularity”, in the sense that over a unit interval there are small steps that are indistinguishable from zero due to the finite number of bits available for a float. The granularity of double precision expressed by machine epsilon is sufficient to represent natural phenomena, and floating point errors can usually be kept under control,

∴	[floatmin(Float32) floatmax(Float32) eps(Float32)]

$[\begin{array}{ccc} 1.1754944 e - 38 & 3.4028235 e 38 & 1.1920929 e - 7 \end{array}]$ (6)

∴	[floatmin(Float64) floatmax(Float64) eps(Float64)]

$[\begin{array}{ccc} 2.2250738585072014 e - 308 & 1.7976931348623157 e 308 & 2.220446049250313 e - 16 \end{array}]$ (7)

∴

Keep in mind that perfect accuracy is a mathematical abstraction, not encountered in nature. In fields such as sociology or psychology 3 digits of accuracy are excellent, in mechanical engineering this might increase to 6 digits, or in electronic engineering to 8 digits. The most precisely known physical constant is the Rydberg constant known to 12 digits, hence a mathematical statement such as

x = 2.6309283450461248350319486319845

is unlikely to have any real significance, while

x = 2.631 \pm 0.0005

is much more informative.

Within the reals certain operations are undefined such as $1 / 0$ . Special float constants are defined to handle such situations: Inf is a float meant to represent infinity, and NaN (“not a number”) is meant to represent an undefinable result of an arithmetic operation.

∴	[1/0 -1.0/0.0 1/Inf -1/Inf Inf/Inf]

$[\begin{array}{ccccc} \infty & - \infty & 0.0 & - 0.0 & N a N \end{array}]$ (8)

∴

Complex numbers $z \in ℂ$ are specified by two reals, in Cartesian form as $z = x + i y$ , $x, y \in ℝ$ or in polar form as $z = ρ e^{i θ}$ , $ρ, θ \in ℝ$ , $ρ ⩾ 0$ . The computer type complex is similarly defined from two floats and the additional constant I is defined to represent $\sqrt{- 1} = i = e^{i π / 2}$ . Functions are available to obtain the real and imaginary parts within the Cartesian form, or the absolute value and argument of the polar form.

∴	z1=1+1im; z2=1-im; [z1+z2 z1/z2]

$[\begin{array}{cc} 2.0 + 0.0 i & - 0.0 + 1.0 i \end{array}]$ (9)

∴	[real(z1) real(z2) real(z1+z2) real(z1/z2)]

$[\begin{array}{cccc} 1.0 & 1.0 & 2.0 & - 0.0 \end{array}]$ (10)

∴	[imag(z1) imag(z2) imag(z1+z2) imag(z1/z2)]

$[\begin{array}{cccc} 1.0 & - 1.0 & 0.0 & 1.0 \end{array}]$ (11)

∴	[abs(z1) abs(z2) abs(z1+z2) abs(z1/z2)]

$[\begin{array}{cccc} 1.4142135623730951 & 1.4142135623730951 & 2.0 & 1.0 \end{array}]$ (12)

∴	[angle(z1) angle(z2) angle(z1+z2) angle(z1/z2)]

$[\begin{array}{cccc} 0.7853981633974483 & - 0.7853981633974483 & 0.0 & 1.5707963267948966 \end{array}]$ (13)

∴

2.Approximation

2.1.Axiom of floating point arithmetic

The reals

ℝ

form an algebraic structure known as a field

(ℝ, +, \cdot)

. The set of floats together with floating point addition and multiplication are denoted as

(𝔽, \oplus, ⊙)

. Operations with floats do not have the same properties as the reals, but are assumed to have a relative error bounded by machine epsilon

where

fl (x) \in 𝔽

is the floating point representation of

x \in ℝ

. The above is often used in error analysis as

Computer number sets are a first example of approximation: replacing some complicated object with a simpler one. It is one of the key mathematical ideas studied throughout this text.

2.2.Cummulative floating point operations

Care should be exercised about the cummulative effect of many floating point operations. An informative example is offered by Zeno's paradox of motion, that purports that fleet-footed Achilles could never overcome an initial head start of

D = 2

given to the lethargic Tortoise since, as stated by Aristotle:

The above is often formulated by considering that the first half of the initial head start must be overcome, then another half and so on. The distance traversed after

N

such steps is

Calculus resolves the paradox by rigorous definition of the limit

D = {lim}_{N \to \infty} D_{N} = 2

and definition of velocity as

v (t) = {lim}_{δ t \to 0} (D (t + δ t) - D (t)) / δ t

δ t = 1 / N

D (t) = 2 [1 - {(1 / 2)}^{t / δ t}]

Undertake a numerical invesigation and consider two scenarios, with increasing or decreasing step sizes

∴	N=10; D=2.0 .^ (0:-1:-N); C=2.0 .^ (-N:1:0); sum(D)==sum(C)

$t r u e$

∴	N=20; D=2.0 .^ (0:-1:-N); C=2.0 .^ (-N:1:0); sum(D)==sum(C)

$t r u e$

Irrespective of the value for

N

D_{N} = C_{N}

in floating point arithmetic. Recall however that computers use binary representations internally, so division by powers of two might have unique features (indeed, it corresponds to a bit shift operation). Try subdividing the head start by a different number, perhaps

π

to get an “irrational” numerical investigation of Zeno's paradox of motion. Define now the distance

S_{N}

traversed by step sizes that are scaled by

1 / π

starting from one to

T_{N}

, traversed by step sizes scaled by

π

starting from

π^{- N}

Again, in the reals the above two expressions are equal,

S_{N} = T_{N}

, but this is no longer verified computationally for all

N

, not even within a tolerance of machine epsilon.

∴	fpi=Float64(pi);

∴	N=10; S=fpi .^ (0:-1:-N); T=fpi .^ (-N:1:0); sum(S)==sum(T)

$t r u e$

∴	N=20; S=fpi .^ (0:-1:-N); T=fpi .^ (-N:1:0); sum(S)==sum(T)

$f a l s e$

∴	sum(S)-sum(T)<eps(Float64)

$f a l s e$

This example gives a first glimpse of the steps that need to be carried out in addition to mathematical analysis to fully characterize an algorithm. Since

S_{N} \neq T_{N}

, a natural question is whether one is more accurate than the other. For some arbitrary ratio

a

, the exact value is known

and can be used to evaluate the errors

| S_{N} - E_{N} |

| T_{N} - E_{N} |

∴	function E(N,a) (1-(1/a)^(N+1))/(1-(1/a)) end;

∴	function εs(N,a) S=a .^ (0:-1:-N) abs(sum(S)-E(N,a)) end;

∴	function εt(N,a) T=a .^ (-N:1:0) abs(sum(T)-E(N,a)) end;

∴	n=30; errs=zeros(Float64,n); errt=zeros(Float64,n);

∴	for i=1:n errs[i]=εs(N,fpi); errt[i]=εt(N,fpi); end

∴	clf(); plot(1:n,errs,1:n,errt,marker="o"); title("Summation error"); grid("on"); xlabel("n"); ylabel("εs,εt"); legend(["εs"; "εt"]);

Note that errors are about the size of machine epsilon for

S_{N}

, but are zero for

T_{N}

, it seems that the summation ordering

T_{N} = a^{- N} + a^{- N + 1} + \dots + 1

gives the exact value. A bit of reflection reinforces this interpretation: first adding small quantities allows for carry over digits to be accounted for.

This example is instructive beyond the immediate adage of “add small quantities first”. It highlights the blend of empirical and analytical approaches that is prevalent in scientific computing.

3.Successive approximations

3.1.Sequences in

ℝ

Single values given by some algorithm are of little value in the practice of scientific computing. The main goal is the construction of a sequence of approximations

{x_{n}}_{n \in ℕ}

that enables assessment of the quality of an approximation. Recall from calculus that

{x_{n}}_{n \in ℕ}

converges to

x

| x_{n} - x |

can be made as small as desired for all

n

beyond some threshold. In precise mathematical language this is stated through:

Definition. ${x_{n}}_{n \in ℕ}$ converges to $x$ if $\forall ε > 0$ , $\exists N (ε)$ such that $| x_{n} - x | < ε$ for $n > N (ε)$ .

Though it might seem natural to require a sequence of approximations to converge to an exact solution

x

Both situations arise when approximating numbers and serve as useful reference points when considering approximation other types of mathematical objects such as functions. For example, the number

π

is readily defined in geometric terms as the ratio of circle circumference to diameter, but can only be approximately expressed by a rational number, e.g.,

π ≅ 22 / 7

. The exact value of

π

is only obtained as the limit of an infinite number of operations with rationals. There are many such infinite representations, one of which is the Leibniz series

Rather, the Leibniz series should be understood as an algorithm, i.e., a sequence of elementary operations that leads to succesively more accurate approximations of

π / 4

Complex analysis provides a convergence proof starting from properties of the

\arctan

function

For

| z | < 1

the sequence

S_{n} = \sum_{k = 0}^{n} {(- z^{2})}^{k}

of partial sums of a geometric series converges uniformly

This elegant result does not address however the points raised above: if

π

were not known, how could the convergence of the sequence

{L_{n}}_{n \in ℕ}

be assessed? A simple numerical experiment indicates that the familiar value of

π

is only recovered for large

n

, with 10000 terms insufficient to ensure five significant digits.

∴	function L(n) L=1.0; s=-1.0 for k=1:n L += s/(2k+1); s = -s end return 4L end

L

∴	[L(100) L(1000) L(10000) Float64(π)]

$[\begin{array}{cccc} 3.1514934010709914 & 3.1425916543395442 & 3.1416926435905346 & 3.141592653589793 \end{array}]$ (14)

3.2.Cauchy sequences

Instead of evaluating distance to an unknown limit, as in

| L_{n} - π | < ε

, one could evaluate if terms get closer to one another as in

| L_{n} - L_{m} | < ε

, a condition that can readily be checked in an algorithm.

Definition. ${x_{n}}_{n \in ℕ}$ is a Cauchy sequence if $\forall ε > 0$ , $\exists N (ε)$ such that $| x_{n} - x_{m} | < ε$ for all m, $n > N (ε)$ .

Note that the distance between any two terms after the threshold

N (ε)

must be smaller than an arbitrary tolerance

ε

. For example the sequence

a_{n} = \sqrt{n}

is not a Cauchy sequence even though the distance between successive terms can be made arbitrarily small

Verification of decreasing successive distance is therefore a necessary but not sufficient condition to assess whether a sequence is a Cauchy sequence. Furthermore, the distance between successive iterates is not necessarily an indication of the distance to the limit. Reprising the Leibniz example, successive terms can be further apart than the distance to the limit, though terms separated by 2 are closer than the distance to the limit (a consequence of the alternating Leibniz series signs)

∴	n=1000; [log10(abs(L(n)-L(n-1))) log10(abs(L(n)-π))]

$[\begin{array}{cc} - 2.6991870973082537 & - 3.000434185835426 \end{array}]$ (15)

∴	[log10(abs(L(n)-L(n-2))) log10(abs(L(n)-π))]

$[\begin{array}{cc} - 5.698969895788488 & - 3.000434185835426 \end{array}]$ (16)

Another question is whether a Cauchy sequence is itself convergent. For sequences of reals this is true, but the Leibniz sequence furnishes a counterexample since it contains rationals and converges to an irrational. Such aspects that arise in number approximation sequences become even more important when considering approximation sequences composed of vectors or functions.

3.3.Sequences in

𝔽

Consideration of floating point arithmetic indicates adaptation of the mathematical concept of convergence is required in scientific computation. Recall that machine epsilon

ϵ

is that largest number such that

1 + ϵ = 1

is true, and characterizes the granularity of the floating point system. A reasonable adaptation of the notion of convergence might be:

Definition. ${x_{n}}_{n \in ℕ}$ , $x_{n} \in 𝔽$ converges to $x \in 𝔽$ if $\forall ε > ϵ$ , $\exists N (ε)$ such that $| x_{n} - x | < ε$ for $n > N (ε)$ .

What emerges is the need to consider a degree of uncertainty in an approximating sequence. If the uncertainty can be bounded to the intrinsic granularity of the number system, a good approximation is obtained.

Number Approximation

1.Numbers

1.1.Number sets

1.2.Quantification

1.3.Computer number sets

2.Approximation

2.1.Axiom of floating point arithmetic

2.2.Cummulative floating point operations

3.Successive approximations

3.1.Sequences in $ℝ$

3.2.Cauchy sequences

3.3.Sequences in $𝔽$

Summary.