Chapter 1 — Mathematical Preliminaries and Error Analysis

This chapter collects the tools the rest of the course is built on. The first half is calculus: the mean value theorems and Taylor’s theorem, which is where nearly every error bound in these notes comes from. The second half is about the machine, where a number is stored in finitely many bits and arithmetic is only approximate. Finally we combine the two, asking how those small errors grow as an algorithm runs and how fast an approximation converges.

1.1 Review of Calculus

Numerical analysis is built almost entirely out of a handful of theorems from calculus. This section collects the ones we will use over and over: the mean value theorems, which let us convert a difference into a derivative, and Taylor’s theorem, which is the source of nearly every error estimate in the course.

Limits and continuity

Limits and Continuity

Definition 1.1

A function \(f\) defined on a set \(X\) of real numbers has the limit \(L\) at \(x_0\), written \(\lim_{x\rightarrow x_0} f(x) = L\), if, given any real number \(\varepsilon>0\), there exists a real number \(\delta>0\) such that \[ |f(x)-L|<\varepsilon,\quad\text{ whenever }\quad x\in X\text{ and } 0<|x-x_0|<\delta. \]

Definition 1.2

Let \(f\) be a function defined on a set \(X\) of real numbers and \(x_0\in X\). Then \(f\) is continuous at \(x_0\) if \[ \lim_{x\rightarrow x_0} f(x) = f(x_0). \] The function \(f\) is continuous on the set \(X\) if it is continuous at each number in \(X\).

The definition says: to force \(f(x)\) within \(\varepsilon\) of \(L\), it is enough to force \(x\) within \(\delta\) of \(x_0\). The picture is a pair of bands. Given the horizontal band of half-width \(\varepsilon\) about \(L\), we must be able to find a vertical band of half-width \(\delta\) about \(x_0\) whose part of the graph lies inside the horizontal one.

Continuity at \(x_0\) asks for something slightly stronger than the limit existing: the limit must exist and equal the value \(f(x_0)\). Both can fail. In Figure 1.2 the one-sided limits disagree, so no limit exists at \(x_0\) and \(f\) cannot be continuous there, whatever value we assign at \(x_0\).

Limits of sequences

Limits of Sequences

Definition 1.3

Let \(\{x_n\}_{n=1}^\infty\) be an infinite sequence of real of complex numbers. The sequence \(\{x_n\}_{n=1}^\infty\) has the limit \(x\) if, for any \(\varepsilon>0\), there exists a positive integer \(N(\varepsilon)\) such that \(|x_n-x|<\varepsilon\), whenever \(n>N(\varepsilon)\). The notation \[ \lim_{n\rightarrow\infty} x_n = x,\text{ or } x_n\rightarrow x\text{ as } n\rightarrow\infty, \] means that the sequence \(\{x_n\}_{n=1}^\infty\) converges to \(x\).

Theorem 1.1

If \(f\) is a function defined on a set \(X\) of real numbers and \(x_0\in X\), then the following statements are equivalent:

\(f\) is continuous at \(x_0\);
If the sequence \(\{x_n\}_{n=1}^\infty\) in \(X\) converges to \(x_0\), then \(\lim_{n\rightarrow\infty} f(x_n) = f(x_0)\).

The sequence version has the same shape, with \(N(\varepsilon)\) playing the role of \(\delta\): given a tolerance \(\varepsilon\), all terms from some index onwards must lie within \(\varepsilon\) of the limit. The notation \(N(\varepsilon)\) is a reminder that the index generally grows as \(\varepsilon\) shrinks.

Boundedness alone is not enough. The sequence \(x_n = 1, 0, 1, 0, \ldots\) stays in \([0,1]\) but never settles down: for \(\varepsilon = \tfrac{1}{4}\) there is no index past which all terms agree to within \(\varepsilon\).

The theorem in the box above is what lets us move between the two notions, and we will use it constantly: to show that an iterative method converges to a root, we show the iterates converge, and continuity then delivers convergence of the function values.

Derivatives

Definition 1.4

Let \(f\) be a function defined in an open interval containing \(x_0\). The function \(f\) is differentiable at \(x_0\) if \[ f'(x_0) = \lim_{x\rightarrow x_0} \frac{f(x)-f(x_0)}{x-x_0} \] exists. The number \(f'(x_0)\) is called the derivative of \(f\) at \(x_0\). A function that has a derivative at each number in a set \(X\) is differentiable on \(X\).

Theorem 1.2

If the function \(f\) is differentiable at \(x_0\), then \(f\) is continuous at \(x_0\).

The quotient \(\bigl(f(x)-f(x_0)\bigr)/(x-x_0)\) is the slope of the secant line through the two points. As \(x\) approaches \(x_0\) the secants rotate towards the tangent, and the derivative is the slope of that limiting line.

Differentiability implies continuity, but not conversely — \(|x|\) at the origin is the standard counterexample.

Derivative theorems

Derivative Theorems

Theorem 1.3 (Rolle’s Theorem)

Suppose \(f\in C[a,b]\) and \(f\) is differentiable on \((a,b)\). If \(f(a)=f(b)\), then a number \(c\) in \((a,b)\) exists with \(f'(c) = 0\).

Theorem 1.4 (Mean Value Theorem)

If \(f\in C[a,b]\) and \(f\) is differentiable on \((a,b)\), then a number \(c\) in \((a,b)\) exists with \[ f'(c) = \frac{f(b)-f(a)}{b-a}. \]

Theorem 1.5 (Extreme Value Theorem)

If \(f\in C[a,b]\), then \(c_1,c_2\in [a,b]\) exist with \(f(c_1)\le f(x) \le f(c_2)\), for all \(x\in [a,b]\). In addition, if \(f\) is differentiable on \((a,b)\), then the numbers \(c_1\) and \(c_2\) occur either at the endpoints of \([a,b]\) or where \(f'\) is zero.

These three are the workhorses. Rolle’s theorem is the special case of the mean value theorem in which the chord is horizontal.

The mean value theorem is the same statement with the picture tilted: the chord from \((a,f(a))\) to \((b,f(b))\) has slope \(\bigl(f(b)-f(a)\bigr)/(b-a)\), and some tangent is parallel to it.

The extreme value theorem is what guarantees that the error bounds we write down are finite: when we bound \(|f^{(n+1)}(\xi)|\) over an interval, it is this theorem that says the maximum is actually attained.

Integrals

Definition 1.5

The Riemann integral of the function \(f\) on the interval \([a,b]\) is the following limit, provided it exists: \[ \int_a^b f(x)\, dx = \lim_{\max \Delta x_i\rightarrow 0} \sum_{i=1}^n f(z_i) \Delta x_i, \] where the numbers \(x_0,x_1,\ldots,x_n\) satisfy \(a=x_0\le x_1 \le\cdots \le x_n = b\), and where \(\Delta x_i = x_i-x_{i-1}\), for each \(i=1,2,\ldots,n\), and \(z_i\) is arbitrarily chosen in the interval \([x_{i-1},x_i]\).

The Riemann integral is the limit of sums of rectangle areas. On each subinterval \([x_{i-1}, x_i]\) we pick any point \(z_i\), use \(f(z_i)\) as the height, and add up \(f(z_i)\,\Delta x_i\). The partition need not be uniform; what matters is that the largest subinterval width tends to zero.

Every quadrature rule in Chapter 4 is a finite version of this picture, with the points \(z_i\) and the weights chosen cleverly rather than arbitrarily.

Integrals

Theorem 1.6 (Weighted Mean Value Theorem for Integrals)

Suppose \(f\in C[a,b]\), the Riemann integral of \(g\) exists on \([a,b]\), and \(g(x)\) does not change sign on \([a,b]\). Then there exists a number \(c\) in \((a,b)\) with \[ \int_a^b f(x)g(x)\, dx = f(c)\int_a^b g(x)\,dx. \]

The special case \(g(x) \equiv 1\) is worth remembering on its own: \[ \int_a^b f(x)\,dx = f(c)\int_a^b 1\,dx = (b-a)\,f(c). \] Geometrically, the area under \(f\) equals the area of a rectangle of width \(b-a\) whose height \(f(c)\) is a value the function actually takes.

Generalizations

Theorem 1.7 (Generalized Rolle’s Theorem)

Suppose \(f\in C[a,b]\) is \(n\) times differentiable on \((a,b)\). If \(f(x)\) is zero at the \(n+1\) distinct numbers \(x_0,\ldots,x_n\) in \([a,b]\), then a number \(c\) in \((a,b)\) exists with \(f^{(n)}(c) = 0\).

Theorem 1.8 (Intermediate Value Theorem)

If \(f\in C[a,b]\) and \(K\) is any number between \(f(a)\) and \(f(b)\), then there exists a number \(c\) in \((a,b)\) for which \(f(c) = K\).

Generalized Rolle is proved by applying ordinary Rolle repeatedly. If \(f\) vanishes at \(n+1\) points, then between each consecutive pair \(f'\) has a zero — that is \(n\) zeros of \(f'\). Applying the same argument to \(f'\) gives \(n-1\) zeros of \(f''\), and after \(n\) steps one zero of \(f^{(n)}\) is left.

This is exactly the argument used in Chapter 3 to derive the error formula for polynomial interpolation.

The intermediate value theorem is the basis of the bisection method: if \(f(a)\) and \(f(b)\) have opposite signs, then \(K=0\) lies between them and a root must exist in \((a,b)\).

Taylor polynomials

Taylor Polynomials

Theorem 1.9 (Taylor’s Theorem)

Suppose \(f\in C^n[a,b]\), that \(f^{(n+1)}\) exists on \([a,b]\), and \(x_0\in [a,b]\). For every \(x\in [a,b]\), there exists a number \(\xi(x)\) between \(x_0\) and \(x\) with \(f(x)=P_n(x)+R_n(x)\), where \[ \begin{aligned} P_n(x) &= f(x_0) + f'(x_0)(x-x_0) + \frac{f''(x_0)}{2!}(x-x_0)^2 + \\ &+ \cdots + \frac{f^{(n)}(x_0)}{n!}(x-x_0)^n = \sum_{k=0}^n \frac{f^{(k)}(x_0)}{k!}(x-x_0)^k \end{aligned} \] and \[ R_n(x) = \frac{f^{(n+1)}(\xi(x))}{(n+1)!}(x-x_0)^{n+1}. \]

Taylor’s theorem is the single most important tool in the course. The polynomial \(P_n\) matches \(f\) and its first \(n\) derivatives at \(x_0\); the remainder \(R_n\) says exactly how much is lost, and it is the remainder — not the polynomial — that produces our error bounds.

Two features of the remainder deserve attention. The factor \((x-x_0)^{n+1}\) means the approximation degrades quickly away from \(x_0\), and the factor \(1/(n+1)!\) means it improves quickly as \(n\) grows. The unknown point \(\xi(x)\) depends on \(x\); in practice we never find it, we only bound \(|f^{(n+1)}(\xi)|\) over the interval of interest.

Example 1.1

Find the maximum error in the approximation \[ \sin x \approx x - \frac{x^3}{3!} + \frac{x^5}{5!}, \qquad -0.3 \le x \le 0.3. \]

Differentiating \(f(x)=\sin x\) repeatedly at \(x_0=0\):

\(k\)	\(f^{(k)}(x)\)	\(f^{(k)}(0)\)
\(0\)	\(\sin x\)	\(0\)
\(1\)	\(\cos x\)	\(1\)
\(2\)	\(-\sin x\)	\(0\)
\(3\)	\(-\cos x\)	\(-1\)
\(4\)	\(\sin x\)	\(0\)
\(5\)	\(\cos x\)	\(1\)
\(6\)	\(-\sin x\)	\(0\)
\(7\)	\(-\cos x\)	—

so the approximation is the Taylor polynomial \(P_5(x)\).

Because \(f^{(6)}(0)=0\), the degree-6 term vanishes and \(P_6 = P_5\). Using \(n=6\) rather than \(n=5\) therefore costs nothing and gives a sharper bound, since the remainder then carries \(x^7\) instead of \(x^6\): \[ P_6(x) = x - \frac{x^3}{3!} + \frac{x^5}{5!}, \qquad R_6(x) = \frac{f^{(7)}(\xi)}{7!}\,x^7 = \frac{-\cos \xi}{7!}\,x^7 . \]

Bounding \(|\cos\xi| \le 1\) and \(|x| \le 0.3\): \[ |R_6(x)| = \frac{|\cos\xi|}{7!}\,|x|^7 \le \frac{1}{7!}\,|x|^7 \le \frac{0.3^7}{7!} = 4.3\times 10^{-8}. \]

Taking \(n=6\) instead of \(n=5\) improved the bound by a factor of about \(0.3\) for no extra work — a small illustration of a general habit worth forming: look for vanishing terms before writing down the remainder.

1.2 Round-off Errors and Computer Arithmetic

Floating point numbers

Writing \(123.456\) digit by digit, with place values \(10^2, 10^1, 10^0, 10^{-1}, 10^{-2}, 10^{-3}\), is fixed point: the decimal point sits in one place. The same number in floating point, or scientific notation, is \[ +1.23456 \times 10^{2}. \]

Computers do this in base \(2\) rather than base \(10\). An integer in binary is read off with place values that are powers of two: \[ 10100101_2 = 2^7 + 2^5 + 2^2 + 2^0 = 128 + 32 + 4 + 1 = 165_{10}. \]

Fractional binary digits carry negative powers, and the same scientific notation applies: \[ -1.0100_2 \times 2^{111_2} = -\bigl(2^0 + 2^{-2}\bigr) \times 2^{\,2^2+2^1+2^0} = -(1.25)\cdot 2^{7} = -160 . \]

IEEE floating point numbers

IEEE Floating Point Numbers

Long real (double precision) format

Widely adopted standard
Default data type in MATLAB, “double” in C
Base 2, 1 sign bit, 11 exponent bits, 52 significand bits: \[ \begin{array}{c|c|c} \text{x} & \text{xxxxxxxxxxx} & \text{xxx}\cdots\text{ (52 total bits) }\cdots\text{xxx} \\ \hline s & c & f \end{array} \]
Represented number: \[ (-1)^s 2^{c-1023}(1+f) \]

A double is 64 bits: one sign bit \(s\), eleven exponent bits \(c\), and fifty-two significand bits \(f\), representing \[ (-1)^s\, 2^{\,c-1023}\,(1+f). \] The exponent is stored with a bias of \(1023\) rather than a sign, and the leading \(1\) of the mantissa is implicit — it costs no storage because a normalized binary mantissa always starts with \(1\).

Example 1.2

The number \(3.0\).

\(s\)	\(c\) (11 bits)	\(f\) (52 bits)
`0`	`10000000000`	`1000` \(\cdots\) `0`

Here \(c = 2^{10} = 1024\) and \(f = 2^{-1}\), so the value is \[ (-1)^0\, 2^{\,1024-1023}\,\bigl(1 + 2^{-1}\bigr) = 1 \cdot 2 \cdot \tfrac{3}{2} = 3.0 . \]

Example 1.3

Machine epsilon, from the neighbours of \(1.0\).

The representation of \(1.0\) has \(c = 01111111111_2 = 1023\) and \(f=0\): \[ (-1)^0 \cdot 2^{\,1023-1023}\cdot(1+0) = 1.0 . \]

The next larger number sets the last significand bit, \(f = 2^{-52}\): \[ (-1)^0\,2^{\,1023-1023}\bigl(1 + 2^{-52}\bigr) = 1 + 2^{-52}, \] so the gap above \(1\) is \[ 2^{-52} \approx 2\times 10^{-16} = \varepsilon . \]

The next smaller number drops the exponent by one and fills the significand with ones, \(c = 1022\): \[ (-1)^0\,2^{\,1022-1023}\bigl(1 + 2^{-1} + 2^{-2} + \cdots + 2^{-52}\bigr) = 2^{-1}\bigl(2 - 2^{-52}\bigr) = 1 - 2^{-53}, \] so the gap below \(1\) is \(2^{-53} = \varepsilon/2\) — half the gap above. The spacing changes at every power of two.

Example 1.4

The gap near \(2^{20}\).

For \(2^{20}\) the exponent field is \(c = 1043\), since \(2^{1043-1023} = 2^{20}\). The next larger number is \[ 2^{20}\bigl(1 + 2^{-52}\bigr) = 2^{20} + 2^{-32}, \] so the gap there is \[ 2^{-32} \approx 2\times 10^{-10}, \] about \(10^{10}\) times larger than the gap at \(1\).

The gaps between representable numbers scale with magnitude. What stays roughly constant is the relative accuracy: about 16 decimal digits, everywhere in the range. This is the single most important fact about floating-point arithmetic.

Decimal floating-point numbers

Value	\(c\)	\(f\)
\(\pm 0\)	all zeros	all zeros
\(\pm\infty\)	all ones	all zeros
`NaN`	all ones	nonzero

Decimal Floating-Point Numbers

Base-10 Floating-Point

For simplicity we study \(k\)-digit decimal machine numbers: \[ \pm0.d_1d_2\ldots d_k \times 10^n, \quad 1\le d_1\le 9, \quad 0\le d_i\le 9 \]
Any positive number within the range can be written: \[ y = 0.d_1d_2\cdots d_kd_{k+1}d_{k+2}\ldots \times 10^n \]
Two ways to represent \(y\) with \(k\) digits:
- Chopping: Chop off after \(k\) digits: \[ fl(y) = 0.d_1d_2\ldots d_k \times 10^n \]
- Rounding: Add \(5 \times 10^{n-(k+1)}\) and chop: \[ fl(y) = 0.\delta_1\delta_2\ldots\delta_k \times 10^n \]

For hand calculation it is easier to work in base 10 with \(k\) digits. The two ways of cutting a number down to \(k\) digits are chopping (discard the rest) and rounding (add \(5\times 10^{n-(k+1)}\), then chop).

Example 1.5

Five-digit representations of \(\pi = 3.1415926535\ldots\)

Rounding to five digits: \[ fl(\pi) = 0.31416 \times 10^{1}. \] Chopping to five digits: \[ fl(\pi) = 0.31415 \times 10^{1}. \]

Errors and significant digits

Errors and Significant Digits

Definition 1.6

If \(p^*\) is an approximation to \(p\), the absolute error is \(|p-p^*|\), and the relative error is \(|p-p^*|/|p|\), provided that \(p\ne 0\).

Definition 1.7

The number \(p^*\) is said to approximate \(p\) to \(t\) significant digits (or figures) if \(t\) is the largest nonnegative integer for which \[ \frac{|p-p^*|}{|p|} \le 5\times 10^{-t}. \]

Example 1.6

Absolute versus relative error.

Take \(p = 1\) and \(p^* = 1 + 10^{-3} = 1.001\). Then \[ \text{absolute error} = |p - p^*| = 10^{-3}, \qquad \text{relative error} = \frac{|p-p^*|}{|p|} = 10^{-3}, \] and since \(10^{-3} \le 5\times 10^{-3}\), the approximation is good to \(3\) significant digits.

Floating point operations and cancellation

Floating Point Operations

Finite-Digit Arithmetic

Machine addition, subtraction, multiplication, and division: \[ \begin{aligned} x \oplus y = fl(fl(x)+fl(y)), \quad & x \otimes y = fl(fl(x)\times fl(y)) \\ x \ominus y = fl(fl(x)-fl(y)), \quad & x \oslash y = fl(fl(x) / fl(y)) \end{aligned} \]
“Round input, perform exact arithmetic, round the result”

Cancellation

Common problem: Subtraction of nearly equal numbers: \[ \begin{aligned} fl(x) &= 0.d_1d_2\ldots d_p\alpha_{p+1}\alpha_{p+2}\ldots\alpha_k \times 10^n \\ fl(y) &= 0.d_1d_2\ldots d_p\beta_{p+1}\beta_{p+2}\ldots\beta_k \times 10^n \end{aligned} \] gives fewer digits of significance: \[ fl(fl(x)-fl(y)) = 0.\sigma_{p+1}\sigma_{p+2}\ldots\sigma_k \times 10^{n-p} \]

The rule is: round the inputs, do exact arithmetic, round the result. Each individual operation is therefore about as accurate as it can be. Trouble comes from how errors combine — above all in subtraction of nearly equal numbers.

Example 1.7

Catastrophic cancellation.

Work with six-digit decimal machine numbers, and take \[ a = 0.1234564, \qquad b = 0.1234546 . \] Rounding each to six digits, \[ fl(a) = 0.123456, \qquad fl(b) = 0.123455 . \]

The machine subtraction is \[ a \ominus b = fl\bigl(fl(a) - fl(b)\bigr) = fl(0.123456 - 0.123455) = fl(0.000001) = 10^{-6}, \] while the true difference is \[ a - b = 1.8 \times 10^{-6}. \]

The absolute error is only \[ |(a-b) - (a \ominus b)| = 8 \times 10^{-7}, \] which looks tiny. But the relative error is \[ \frac{|1.8\times 10^{-6} - 10^{-6}|}{1.8\times 10^{-6}} = 0.44, \] that is \(44\%\): the result has one significant digit, though both inputs had six.

Nothing went wrong in the subtraction itself; it was exact. The leading digits of \(a\) and \(b\) agreed and cancelled, promoting the inputs’ round-off error from the last digit to the first.

Example 1.8

The order of operations matters.

Exactly, \[ (96384 + 26.678) - 96410 = 0.678 . \]

With 5-digit rounding, \(96384 + 26.678 = 96410.678\) rounds to \(96411\), so \[ 96411 - 96410 = 1, \] an error of about \(47\%\). With 5-digit chopping it rounds to \(96410\), and \[ 96410 - 96410 = 0, \] losing the answer entirely.

Rewriting the expression to subtract the nearly equal numbers first, \[ (96384 - 96410) + 26.678 = 0.678, \] and now with 5-digit rounding or chopping alike, \[ -26 + 26.678 = 0.678, \] which is exact. The two forms are algebraically identical and numerically worlds apart.

1.3 Algorithms and Convergence

Error growth

Error Growth and Stability

Definition 1.8

Suppose \(E_0>0\) is an initial error, and \(E_n\) is the error after \(n\) operations.

\(E_n\approx CnE_0\): linear growth of error
\(E_n\approx C^nE_0\): exponential growth of error

Stability

Stable algorithm: Small changes in the initial data produce small changes in the final result
Unstable or conditionally stable algorithm: Large errors in final result for all or some initial data with small errors

Example 1.9

Linear growth.

Consider computing \[ S_n = \sum_{i=1}^n \frac{1}{i}, \] and suppose each division introduces an error \(\varepsilon\). The computed sum is \[ \hat{S}_n = \sum_{i=1}^n \left(\frac{1}{i} + \varepsilon\right) = S_n + n\varepsilon, \] so the absolute error is \[ E_n = |\hat{S}_n - S_n| = n\varepsilon . \] The error grows linearly in the number of operations. This is the benign case, and about the best one can hope for.

Example 1.10

Exponential growth.

Consider the linear recurrence relation (or difference equation) \[ p_n = \tfrac{1}{2}\bigl(3p_{n-2} - 5p_{n-1}\bigr), \] with \(p_0, p_1\) given as initial conditions. Its general solution is \[ p_n = c_1\left(\tfrac{1}{2}\right)^n + c_2(-3)^n, \] with \(c_1, c_2\) determined by the initial conditions (see Math 10B or Math 55 for the standard theory).

Choosing \(c_1 = 1\), \(c_2 = 0\) gives the true solution \[ p_n = \frac{1}{2^n} \longrightarrow 0 . \]

Now perturb the initial data slightly, so that \(c_1 = 1\) but \(c_2 = \varepsilon\) for some tiny \(\varepsilon\): \[ \hat{p}_n = \frac{1}{2^n} + \varepsilon\,(-3)^n . \] The error is \[ E_n = |\hat{p}_n - p_n| = \varepsilon\,3^n, \] which grows exponentially. However small \(\varepsilon\) is, the \((-3)^n\) term eventually swamps the decaying solution: with \(\varepsilon = 10^{-16}\) the two terms are comparable by about \(n = 17\), and after that the computed answer is pure noise.

The recurrence is an unstable algorithm for computing \(2^{-n}\), even though the underlying problem is perfectly well behaved.

The distinction matters because the algorithm, not the problem, is at fault here. A stable algorithm keeps small input changes producing small output changes; an unstable one amplifies them, and no amount of precision rescues it.

Rate of convergence for sequences

Rate of Convergence (Sequences)

Definition 1.9

Suppose \(\{\beta_n\}_{n=1}^\infty\) is a sequence converging to zero, and \(\{\alpha_n\}_{n=1}^\infty\) converges to a number \(\alpha\). If a positive constant \(K\) exists with \[ |\alpha_n-\alpha|\le K|\beta_n|,\quad\text{for large }n, \] then we say that \(\{\alpha_n\}_{n=1}^\infty\) converges to \(\alpha\) with rate of convergence \(O(\beta_n)\), indicated by \(\alpha_n = \alpha + O(\beta_n)\).

Polynomial rate of convergence

Normally we will use \[ \beta_n = \frac{1}{n^p}, \] and look for the largest value \(p>0\) such that \(\alpha_n = \alpha + O(1/n^p)\).

Knowing that \(\alpha_n \to \alpha\) is rarely enough; we want to know how fast. The \(O\) notation compares the error against a reference sequence \(\beta_n\), usually \(1/n^p\), and we look for the largest \(p\) that works.

Example 1.11

Reading off the rate directly.

For \[ \alpha_n^1 = 3 + \frac{1}{n}, \qquad \alpha_n^2 = 3 + \frac{1}{n^2}, \] both converge to \(3\), with \[ \alpha_n^1 = 3 + O\!\left(\tfrac{1}{n}\right), \qquad \alpha_n^2 = 3 + O\!\left(\tfrac{1}{n^2}\right). \] The second is far better: to gain one decimal digit the first needs ten times as many terms, the second only about three times as many.

Example 1.12

Using Taylor’s theorem to find a rate.

Let \(\alpha_n = \cos\frac{1}{n^2}\), which tends to \(\alpha = 1\). Expanding \(\cos\) about \(0\) with Taylor’s theorem, \[ \left|\alpha_n - 1\right| = \left|\cos\tfrac{1}{n^2} - 1\right| = \left|\frac{f''(\xi)}{2!}\left(\frac{1}{n^2}\right)^{2}\right| \le \frac{1}{2}\cdot\frac{1}{n^4}, \] using \(|f''(\xi)| = |\cos\xi| \le 1\). The rate of convergence is therefore \(O\!\left(1/n^4\right)\).

Example 1.13

A rational sequence.

Let \[ \alpha_n = \frac{n^3 + n - 1}{2n^3} \longrightarrow \frac{1}{2}. \] Then \[ \left|\alpha_n - \frac{1}{2}\right| = \left|\frac{n^3 + n - 1 - n^3}{2n^3}\right| = \frac{n-1}{2n^3} \le \frac{1}{2n^2}, \] so \[ \alpha_n = \frac{1}{2} + O\!\left(\frac{1}{n^2}\right). \]

Rate of convergence for functions

Rate of Convergence (Functions)

Definition 1.10

Suppose that \(\lim_{h\rightarrow 0} G(h) = 0\) and \(\lim_{h\rightarrow 0} F(h) = L\). If a positive constant \(K\) exists with \[ |F(h)-L| \le K | G(h) |,\quad\text{for sufficiently small }h, \] then we write \(F(h) = L+O(G(h))\).

Polynomial rate of convergence

Normally we will use \[ G(h) = h^p, \] and look for the largest value \(p>0\) such that \(F(h) = L+O(h^p)\).

The same idea applies as \(h \to 0\), with \(G(h) = h^p\) as the reference. This is the form we will use for every discretization method in the course: a difference formula, a quadrature rule or an ODE solver is \(O(h^p)\), and \(p\) is its order of accuracy.

Example 1.14

A first-order difference quotient.

Let \(F(h) = (e^h - 1)/h\), which tends to \(1\) as \(h \to 0\). Expanding the exponential, \[ \frac{e^h - 1}{h} = \frac{\bigl(1 + h + \tfrac{h^2}{2} + \cdots\bigr) - 1}{h} = 1 + \frac{h}{2} + \cdots = 1 + O(h). \] The convergence is first order: halving \(h\) halves the error.