Sequences and Series
Series expansion
Conditional execution can be combined with a jump to produce a terminating loop. Convince yourself that the following code sums the first 100 terms of the Taylor expansion of exp(0.752) =2.12123825347001. Remember that \(\mathsf{e^x = 1 + x + x^2/2! + x^3/3! + \cdots}\)
const double x = 0.752;
unsigned long int n = 1;
double sum = 1.0;
double term = x;
my_loop:
sum += term;
term *= (x/++n);
if (n <= 100) goto my_loop;
cout << sum << endl;
Brief aside: Why would it be a bad idea to implement the expansion in the following way?
const double x = 0.752;
unsigned long int n = 1;
double sum = 1.0;
double x_pow_n = x;
unsigned long int n_fact = n;
repeat_loop:
sum += x_pow_n/n_fact;
x_pow_n *= x;
n_fact *= ++n;
if (n < 100) goto repeat_loop;
cout << sum << endl;
First, note that 12! = 479001600 and 13! = 6227020800. The largest unsigned long int
that can be represented, however, is 232 = 4294967296. Thus, the variable n_fact will
overflow as n goes from 12 to 13.
Here is a function that implements the Taylor series of \(\mathsf{exp(x)}\) to arbitrary order. For efficiency, the lowest three cases are coded explicitly.
double taylor_exp(double x, unsigned int max_order)
{
if (max_order == 0) return 1.0;
if (max_order == 1) return 1.0+x;
if (max_order == 2) return 1.0+x*(1.0+0.5*x);
double sum = 1.0, term = x;
for (unsigned long int n = 2; n < max_order+2; ++n)
{
sum += term;
term *= x/n;
}
return sum;
}
Pay special attention to the test expression (n < max_order+2). Why does it have this form?
Recursively defined sequences
Suppose we want to compute the Fibonacci sequence. The definition is a recursive one: \(\mathsf{x_0 = 0}\), \(\mathsf{x_1 = 1}\), and \(\mathsf{x_n = x_{n-1} + x_{n-2}}\).
If we need to store all of the terms in the sequence, we might try something like the following.
int x[100];
x[0] = 0;
x[1] = 1;
for (int n = 2; n < 100; ++n)
x[n] = x[n-1] + x[n-2];
The drawback to this approach is that we have to decide in advance how many terms we want to keep and declare an appropriately-size buffer array. Alternatively, if we simply need to produce (rather than store) the sequence, then no array is necessary.
int x_p = 0; // penultimate term
int x_u = 1; // ultimate term
cout << x_p << endl;
while (x_u > x_p)
{
cout << x_u << endl;
const int tmp = x_u + x_p;
x_p = x_i;
x_u = tmp;
}
This sequence will terminate when the int variables overflow: i.e., when the
terms in the Fibonacci sequence become too large to represent as a 32-bit integer type.
Here’s an approach in which the elements of the sequence are generated on demand and stored in a lookup table for later use:
#include <iostream>
using std::cout;
using std::endl;
#include <vector>
using std::vector;
size_t fib(size_t n)
{
// static causes val to persist across function calls
static vector<size_t> val;
if (val.size() == 0)
{
val.push_back(0);
val.push_back(1);
}
while (n+1 > val.size())
{
const size_t xk = val.back();
val.pop_back();
const size_t xkm1 = val.back();
const size_t xkp1 = xkm1 + xk;
val.push_back(xk);
val.push_back(xkp1);
}
return val[n];
}
int main()
{
cout << fib(0) << endl; // populate x0=0, x1=1
cout << fib(1) << endl; // lookup x1
cout << fib(15) << endl; // compute x2 .. x15
cout << fib(10) << endl; // loopup x10
cout << fib(20) << endl; // compute x16 .. x20
cout << fib(40) << endl; // compute x21 .. x40
cout << fib(30) << endl; // lookup x30
cout << fib(50) << endl; // compute x41 .. x50
return 0;
}
Periodicity
Q: What is the period \(\mathsf{P}\) of the repeating sequence
defined by \(\mathsf{x_0 = 5}\), \(\mathsf{x_{n+1} = (2x_n + 7) \!\mod 100}\)?
First, we need to think about how to identify the period. We observe that \(\mathsf{(x_n)}\) is constructed in such a way that each successive element depends only on the value of the previous one, according to a rule \(\mathsf{F(x) = (2x+7) \mod 100}\). This means that \(\mathsf{x_1 = F(x_0)}\), \(\mathsf{x_2 = F(x_1)}\), etc. It also means that if \(\mathsf{x_P = x_0}\), then \(\mathsf{x_{P+1} = F(x_P) = F(x_0) = x_1}\), ensuring that the sequence cycles over the values \(\mathsf{x_0, x_1, \ldots , x_{P-1}}\) endlessly. Hence, the lowest integer \(\mathsf{P}\) such that \(\mathsf{x_P = x_0}\) is the period.
Our algorithm can be summarized as follows:
start from \(\mathsf{x_0}\)
generate subsequent elements one by one, keeping a counter
stop if the current element equals \(\mathsf{x_0}\)
report the counter value
Here’s one possible C++ implementation:
#include <iostream>
using std::cout;
using std::endl;
int next(int x) { return (3*x+7)%100; }
int main()
{
const int x0 = 5;
int period = 1;
int x = x0;
while ( (x = next(x)) != x0 ) ++period;
cout << "The period is " << period << endl;
return 0;
}
It’s useful to mentally trace out the program execution for the extreme cases.
What if the very first call to next assigns x a value that is equal to x0?
Verify that the reported period is 1.
Q: How many unique elements appear in the sequence defined by \(\mathsf{x_0 = 1}\), \(\mathsf{x_1 = 2}\), \(\mathsf{x_{n+1}= (2x_n+ x_{n-1}) \!\mod 50}\)?
Identifying the period here is a little more subtle. Each element in the sequence now depends on the two previous ones. The sequence doesn’t start to repeat until both \(\mathsf{x_P = x_0}\) and \(\mathsf{x_{P+1} = x_1}\).
As a first step, let’s identify the elements of the sequence up to the point where it starts repeating:
#include <iostream>
using std::cout;
using std::endl;
void advance(int &x0, int &x1)
{
const int xx = (2*x1+x0)%50;
x0 = x1;
x1 = xx;
}
int main()
{
const int x0 = 1;
const int x1 = 2;
int x = x0;
int xx = x1;
cout << "0\t" << x << endl;
cout << "1\t" << xx << endl;
int period = 0;
do
{
advance(x,xx);
++period;
cout << period+1 << "\t" << xx << endl;
} while ( !( x == x0 and xx == x1) );
// while (x != x0 or xx != x1);
cout << "The period is " << period << endl;
return 0;
}
The output from this program is
0 1
1 2
2 5
3 12
4 29
5 20
. .
. .
57 48
58 1
59 0
60 1
61 2
The period is 60
We now need to devise a strategy for determining the unique elements. One
possibility is to keep a vector<int> for storage and add to it each
new element of the sequence that it does not already contain. We can add
elements using the push_back method, in which case the vector is dynamically
resized. We just need to write a function that checks whether a given element
has been encountered before.
#include <iostream>
using std::cout;
using std::endl;
#include <cstddef>
using std::size_t;
#include <vector>
using std::vector;
void advance(int &x0, int &x1)
{
const int xx = (2*x1+x0)%50;
x0 = x1;
x1 = xx;
}
vector<int> unique_elements;
bool already_found(int x)
{
for (size_t i = 0; i < unique_elements.size(); ++i)
if (x == unique_elements[i]) return true;
return false;
}
int main()
{
const int x0 = 1;
const int x1 = 2;
unique_elements.push_back(x0);
unique_elements.push_back(x1);
int x = x0;
int xx = x1;
int period = 0;
do
{
advance(x,xx);
++period;
if (!already_found(xx)) unique_elements.push_back(xx);
} while ( !( x == x0 and xx == x1) );
cout << "The period is " << period << endl;
cout << "The number of unique elements is " << unique_elements.size() << endl;
return 0;
}
Here is the output from this program:
The period is 60
The number of unique elements is 30
An alternate strategy is to keep an array of boolean values that flag whether each
element appears in the sequence. Since the sequence is defined \(\mathsf{\textsf{mod}~50}\), an
array with indices 0, 1, …, 49 is sufficient:
#include <iostream>
using std::cout;
using std::endl;
#include <cstddef>
using std::size_t;
#include <vector>
using std::vector;
void advance(int &x0, int &x1)
{
const int xx = (2*x1+x0)%50;
x0 = x1;
x1 = xx;
}
bool found[50];
int main()
{
const int x0 = 1;
const int x1 = 2;
for (int i = 0; i < 50; ++i) found[i] = false;
found[x0] = true;
found[x1] = true;
int x = x0;
int xx = x1;
int period = 0;
do
{
advance(x,xx);
++period;
found[xx] = true;
} while ( !( x == x0 and xx == x1) );
cout << "The period is " << period << endl;
int count = 0;
for (int i = 0; i < 50; ++i)
if (found[i]) ++count;
cout << "The number of unique elements is " << count << endl;
return 0;
}
Here is a version that uses the set<int> data structure:
#include <iostream>
using std::cout;
using std::endl;
#include <cstddef>
using std::size_t;
#include <set>
using std::set;
void advance(int &x0, int &x1)
{
const int xx = (2*x1+x0)%50;
x0 = x1;
x1 = xx;
}
set<int> unique_elements;
int main()
{
const int x0 = 1;
const int x1 = 2;
unique_elements.insert(x0);
unique_elements.insert(x1);
int x = x0;
int xx = x1;
int period = 0;
do
{
advance(x,xx);
++period;
unique_elements.insert(xx);
} while ( !( x == x0 and xx == x1) );
cout << "The period is " << period << endl;
cout << "The number of unique elements is " << unique_elements.size() << endl;
return 0;
}
Finally, we could also try to solve this by piggy-backing on pre-existing UNIX tools. The following program just dumps the elements of the sequence to the terminal.
#include <iostream>
using std::cout;
using std::endl;
void advance(int &x0, int &x1)
{
const int xx = (2*x1+x0)%50;
x0 = x1;
x1 = xx;
}
int main()
{
const int x0 = 1;
const int x1 = 2;
int x = x0;
int xx = x1;
cout << x << endl;
cout << xx << endl;
do
{
advance(x,xx);
cout << xx << endl;
} while ( !( x == x0 and xx == x1) );
return 0;
}
Supposing that this program file is named q2.cpp, we can compile it to an
executable q2 and then pipe its output to uniq and then to wc
$ g++ -o q2 q2.cpp
$ ./q2 | uniq | wc -l
30
Sequences and series
Most interesting computer calculations can be understood as the asymptotic limit of some algebraic process. The goal is to determine when such a process has converged to within the desired accuracy. The relevant mathematical concepts are sequences and series, and the amount of computational work that has been expended is analogous to how many elements along we are.
A sequence \(\mathsf{(a_n) = (a_0, a1, a2, \ldots)}\) tends to a finite limit \(\mathsf{L}\) if its residuals \(\mathsf{r_n = L-a_n}\) can be made arbitrarily small simply by choosing \(\mathsf{n}\) large enough. More formally, we say that
iff for any \(\mathsf{\epsilon > 0}\) we can find \(\mathsf{N(\epsilon)}\) such that \(\mathsf{|r_n| = |L-a_n| < \epsilon}\) for all \(\mathsf{n > N}\).
To say that an infinite series converges is to assert that its sequence of partial sums—\(\mathsf{(s_n)}\) with
—has a finite limit: i.e.,
One can formalize various tests for convergence, but these tend not to be useful. For practical purposes, the interesting question is not whether a sequence converges but how fast. In real-world settings, we typically want to estimate L from a finite data set (knowledge of only a certain number of terms in the sequence). We have a better chance of doing this well if we have some idea of how the sequence approaches L at large n. In other words, what is the functional form of the residuals?
The residuals are typically \(\mathsf{r_n=O(n^{-p})}\) (polynomial, with \(\mathsf{p>0}\)) or \(\mathsf{r_n=O(z^n)}\) (exponential with \(\mathsf{0<z<1}\)). In the common mathematical nomenclature, we classify the rate of convergence according to the asymptotic behaviour of the ratio of adjacent residuals. If that ratio tends to a number in the interval (0,1)
then the convergence is said to be linear. The endpoints of the interval correspond to so-called superlinear (0) and sublinear (1) behaviour. The convergence is called quadratic or cubic if
with \(\mathsf{q = 2}\) or \(\mathsf{q = 3}\).
In this language, polynomial dependence of the residuals is sublinear and is considered slow convergence. Exponential dependence of the residuals is linear and is considered fast.
Let’s use the following notation \(\mathsf{R(n)}\) to denote the functional form of the \(\mathsf{n}\)th residual’s magnitude. The most naive way to approximate the sequence limit is to choose an accuracy threshold and then compute a representative element with \(\mathsf{n}\) sufficiently large:
Depending on the accuracy requirement and how \(\mathsf{R(n)}\) scales, \(\mathsf{N}\) may be very large indeed—perhaps even too large to calculate in reasonable CPU time. The most common scaling behaviours lead to the following requirements:
Keep in mind that \(\mathsf{N(\epsilon) \sim \epsilon^{-1/p}}\) is much more work (slow convergence) than \(\mathsf{N(\epsilon) \sim log(\epsilon)}\) (fast convergence). An example might help put this in context. If we’ve computed a sequence limit to some threshhold \(\mathsf{\epsilon}\) by computing \(\mathsf{a_N}\) and we want to add two more decimal digits of accuracy (\(\mathsf{\epsilon \to \epsilon/100}\)), then the polynomial case requires that we go to \(\mathsf{a_{100N}~(p=1)}\) or \(\mathsf{a_{10N}~(p=2)}\), whereas the exponential case requires \(\mathsf{a_{N+m}}\) for some \(\mathsf{m \ll N}\).
Sequence transformations
A sequence that is slowly converging can often be transformed into another sequence that converges to the same limit but gets there faster. In other words, given a sequence \(\mathsf{(a_n)}\) we want to construct another, call it \(\mathsf{(b_n)}\), such that
and
This can be done in a mathematically trivial way by considering an infinite subsequence, e.g.,
This is nothing more than constructing a new sequence by picking out terms 2, 4, 8, 16 and so on from the old sequence. This doesn’t gain us anything in terms of computational expenditure. But, formally at least, it does have the effect of transforming a sequence with polynomial residuals—e.g.
—into one with exponential residuals:
In the former case, our best strategy so far has been to plot the known values of \(\mathsf{(a_n)}\) against \({n^{-p}}\) and to use curve-fitting techniques to extrapolate those points to the y-axis, where the intercept is an estimate of \(\mathsf{L}\). We also considered a nonlinear sequence transformation called the Aitken process, defined by
What we found was that a convergent sequence with residuals \(\mathsf{r_n = Cn^{-1/p}}\) was mapped to a new sequence with residuals \(\mathsf{|L-A_n| = [C/(p+1)]n^{-1/p}}\). That is to say, the functional form remained polynomial and the convergence was only improved by a \(\mathsf{p+1}\) rescaling of the prefactor. On the other hand, the efficacy of the Aitken process for sequences with exponentially converging residuals is rather dramatic.
In the special case of the geometric series, whose partial sums form a sequence
of the form \(\mathsf{s_n = L + Cz^n}\) (with \(\mathsf{L=1/(1-z)}\) and \(\mathsf{C=z/(1-z)}\)) with no higher-order corrections, the Aitken process is exact. It maps the partial sums \(\mathsf{(s_0, s_1, s_2, \ldots)}\) onto the perfectly converged sequence \(\mathsf{(A_n) = (L, L, L, \ldots)}\). This can be demonstrated by direct calculation:
A sequence \(\mathsf{(a_n)}\) that is close to geometric will behave similarly but its residuals will likely have additional terms that scale as \(\mathsf{z^{2n}}\), \(\mathsf{z^{3n}}\), etc. In that case, the Aitken process removes the part that is leading order in \(\mathsf{z^n}\).
Richardson extrapolation
Attempts to compute derivatives numerically inevitably start from the Taylor series expansion
We know that the finite difference
formally converges to the desired derivative in the limit \(\mathsf{h\to 0}\), but in practice the limitations of the floating-point representation render that limit unworkable. Sufficiently small \(\mathsf{h}\) leads to a subtraction of nearly equal quantities in the numerator and hence catastrophic loss of significance. In general, we should be very careful of any calculation that behaves like 0/0.
So how to improve on \(\mathsf{D^{(1)}(h)}\)? The simplest trick is to use a finite difference that is symmetric about \(\mathsf{x}\), which reduces the error to \(\mathsf{O(h^2)}\). A more systematic approach is to build linear combinations designed to eliminate the errors at the next highest order:
This leads to a heirarchy of approximations to \(\mathsf{f(x)}\). The key point, however, is that the step sizes \(\mathsf{h}\), \(\mathsf{h/\xi}\), \(\mathsf{h/\xi^2}\), \(\mathsf{h/\xi^3}\), … (\(\mathsf{\xi > 1}\)) at which the finite differences are evaluated are not necessarily small. A good choice of \(\mathsf{h}\) in this scheme may be very far from infinitesimal. The progress at each stage arises from systematically cancelling off the errors at increasing order. The results are often arranged in a triangular tableau.
If the Richardson extrapolation succeeds, the results should converge rapidly from left to right. The rightmost columns will have accuracy that far exceeds the input points in the far left column, and \(\mathsf{f'(x)}\) can be extremely well approximated even if none of the step sizes are particularly small.
Of course, in many situations you may have no control over where \(\mathsf{f(x)}\) can be evaluated. You may only have tabulated data and no access to or knowledge of \(\mathsf{f(x)}\) at all. It is still possible to implement the Richardson scheme using discrete points. Given \(\mathsf{N}\) finite differences
defined for the points \(\mathsf{\{ x_0, x_1, \ldots , x_N \}}\), we can construct \(\mathsf{N-1}\) improved differences
More generally, at order \(\mathsf{n}\) there are differences \(\mathsf{D^{(n)}_k}\) defined for each of \(\mathsf{k = 1, 2, \ldots , N+1-n}\).