A few notes on computation

Section 2.11 A few notes on computation

Many of the results in these notes have been illustrated on the computer and some of the exercises require a computational approach. Whenever using the computer, it is always wise to examine the results critically. Here's a simple numerical example where things clearly go awry. In it, we are iterating the function

f (x) = x^{2} - 9.1 x + 1

from the fixed point

x_{0} = 0.1 .

We should generate just a constant sequence.

xxxxxxxxxx
 
x = 0.1
for i in range(20):
    x = x**2 - 9.1*x+ 1
    print(x)

🔗

Uh-oh!

🔗

It must be understood that this is a simple consequence of the nature of floating point arithmetic. Part of the issue is that the decimal number

0.1

1 / 10

is not exactly representable in binary. In fact,

\frac{1}{10} = 0_{\dot{2}} 0 \overset{―}{0011} = \frac{1}{2} \sum_{k = 1}^{\infty} \frac{3}{16^{k}} .

🔗

Thus, the computer must introduce round-off error in the computation. Furthermore,

0.1

is a repelling fixed point of the function. Thus, that round-off error is magnified with each iteration. Our study of dynamics has illuminated a critical issue in numerical computation!

🔗

In the terminology of numerical analysis, computation near an attractive fixed point is stable while computation near a repelling fixed point is unstable. Generally speaking, stable computation is trustworthy while unstable computation is not. The implication for the pictures that we see here is that illustration of attractive behavior should be just fine. In Figure 2.3.1, for example, images (b) and (c) illustrate attraction to a fixed point that not only involves stable computation but also agrees with our theoretical development. We are happy with those figures. The cobweb plot shown in Figure 2.3.1 (d), however, should frankly be viewed with some suspicion.

🔗

The same can be said for the bifurcation diagram in Figure 2.7.2. In much of that image, we see a gray smear indicating chaos. How can we trust that? Well, first, theory tells us that there really is chaos. That is, there are orbits that are dense in some interval for many

c

values. Furthermore, much of the image shows attractive regions and we can be confident in that portion.

🔗

In fact, in many of the images that we will generate later - Julia sets, the Mandelbrot set, and similar images - the stable region dominates. Thus, we can be confident in overall image because the unstable region is the complement of the stable region. We might not be confident in computations involving some particular point, but we can be confident in the overall picture. (This will, perhaps, be more clear as we move into complex dynamics.)

🔗

Nonetheless, sometimes we want to experiment with genuinely unstable dynamics. One way to improve our confidence in these kinds of computations is to use high precision numbers. Consider, for example, the cobweb plot of the doubling map shown in Figure 2.9.1. A naive approach to generate the first few terms of an orbit associated with the doubling map might be as follows:

xxxxxxxxxx
 
import numpy as np
x = 1/np.pi
for i in range(55):
    x = 2*x%1
    print(x)

🔗

We've reached the fixed point zero and now we're stuck! Even if we iterate 1000 times, we'll generate a cobweb plot that looks like Figure 2.11.1

🔗
Figure 2.11.1. An inaccurate cobweb plot

🔗

The cobweb plot shown in Figure 2.9.1 was generated using the mpmath multi-precision library for Python with code that looked something like so:

xxxxxxxxxx
 
from mpmath import mp
mp.prec = 1000
x = 1/mp.pi
for i in range(1000):
    x = 2*x%1
    print(float(x))

🔗

While the truncated output looks the same, note that this was after 1000 iterates. This behavior makes perfect sense if you understand that the doubling map loses one bit of precision with every iterate.