Voyages in Dragon Country

by Sylvia Camacho

mathematics and statistics programming techniques

Transcribed from the printed issue; not yet reviewed. Read it in the PDF of the issue, page 60

Marriage to the editor of Vector confers some privileges, one of which is to see the copy before it goes to print. Thus I am able to refer to Gérard Langlet’s letter, included in this issue. This continues the attack that Gérard made in 11.2 (‘Chaos Revisited’) on the uncritical use of computers to run iterative algorithms. Gérard illustrates the problems found with IEEE floating point arithmetic by reference to different formulations of the logistic equation, this being quoted in most works on Chaos as a paradigm of the effects being discussed. However, as Donald McIntyre points out in 11.4 (‘The Perils of Subtraction’), older and much-loved algorithms such as that of Archimedes for calculating \(\pi\) are also sabotaged by standard computer number representation.

A Storm Brewing¶

In the current issue Gérard suggests that error checking, as supported by Mathematica, should be used and quoted whenever propagation of truncation errors is possible and that failure to do so may invalidate many respectable academic papers. Not unexpectedly, he reports that this suggestion is raising a storm among francophones. While I love a good fight and shall do my best to get to the battleground of this one, my own interest is less technical and mathematical and more philosophical. So in this article I am going to tilt at my own choice of dragon, albeit one from the same lair.

Look who Whistled up the Wind¶

ZENO of Elea (c. 490–c. 435 B.C.). Greek philosopher and mathematician. Propounded paradoxes to point out the problem of the relation of the discrete to the continuous……

Zeno’s paradox of Achilles and the tortoise. A tortoise has a head start on Achilles equal to the distance from a to b and both start running, Achilles after the tortoise. Although Achilles runs faster than the tortoise, he would never catch up with the tortoise, since, while Achilles goes from a to b, the tortoise goes from b to c, and while Achilles goes from b to c, the tortoise goes from c to d, etc., this process never ending. The explanation of the fallacy is that motion is measured by space intervals per unit of time, not by numbers of points. If Achilles takes time \(t_1, t_2, t_3, \cdots\) to go from a to b, b to c, c to d, ···, then Achilles will catch the tortoise in time equal to \(\sum\limits_{1}^{\infty} t_i\), if this sum is finite. If the tortoise travels 10 feet per second and Achilles 20 feet and the tortoise starts 10 feet in advance, Achilles will catch him at the end of the first second, since then \(t_1 = \dfrac{1}{2},\ t_2 = \dfrac{1}{4},\ \cdots,\ t_n = \dfrac{1}{2}^{n},\ \cdots\) and \(\dfrac{1}{2} + \dfrac{1}{4} + \dfrac{1}{8} + \cdots = 1\).

James & James "Mathematics Dictionary" Van Nostrand Reinhold NY ed 5 1992

The Dragon at the Eye of the Storm¶

Now whenever in a mathematical expression I see ‘··· = n’, the hairs on the back of my neck start to prickle. What, I ask myself, is finite about this procedure? I do not of course quarrel with the suggestion that in the example quoted the limit is 1 but how, in a finite time, could I compute this result using this algorithm? Also, in the real world, what counts as Achilles drawing level with the tortoise? Let us imagine the tortoise followed by Achilles striding along an (imaginary) time-line. Has Achilles drawn level when his leading edge touches the trailing edge of the tortoise? We know that in the real world they both have very fuzzy edges, the fuzziness depending on the resolution of our observing equipment.

To many readers all this will read like a philosophic quibble but remember that Chaos gurus aspire to model phenomena which are not amenable to Calculus. They are experts in discontinuity and the destabilising effects of small but precisely expressed differences. They ought therefore to be acutely aware of the minimum difference which is appropriate to the model they are investigating. In response to similar remarks, Frank Evans (long time APL user) said, ‘one must always be aware of the fineness of one’s grid’. So before accepting that an algorithm is a plausible model of any real-world system one must ask the crucial questions: ‘What is the resolution of your lens?’, ‘What is the wavelength of your light?’, ‘What is the grain size of your film?’, ‘What is the pixel count of your screen?’, ‘What are the limits of accuracy of your computer in relation to the algorithm you are computing?’ A statistician would say, ‘what is your population?’. So I ask myself ‘what do we populate our computers with?’ and the answer is ‘bits’. So let us declare a maximum population, say p. How might we organise p bits to emulate the logistic equation, as cited by Langlet, \(\mathbf{x_{n+1} = 4x_n(1-x_n)}\)?

We have to set a starting value between zero and p. Let us call this v. So the equivalent of \(\mathbf{1-x_n}\) is p-v. We can imagine the starting value as a set of bits v⍴1. Multiplication is then equivalent to replication 4/(p-v)/v⍴1 but this number of bits exceeds the maximum population so we must emulate division in order to restore the original ‘fineness of the grid’. The resolution of the grid is p and so we use this to set up a template to fit the result of the replication:

template←(+/4/(p-v)/v⍴1)⍴⌽p↑1

This gives a boolean vector of the same length as the replicated bits, all zeros except for one bit at each interval of p. The number of bits not equal to zero in this template will be the value of the next generation of the iterated equation. The result is always integer and is equivalent to taking the floor of a division by p. To take the ceiling one would reshape p↑1. To round the result the position of the bit set to 1 would be at the midpoint of p↑0. The other possible positions for the 1 do not have an obvious interpretation but might be interesting to speculate about.

To many readers this will seem like a convoluted way of doing a bit of simple integer arithmetic but I hope that it serves to illuminate the problem that arises when the logistic equation is iterated using floating point. This is that the ‘fineness of the grid’ is made finer with every iteration, in a doomed pursuit of precision. This cannot be consonant with the use of \(\mathbf{1-x_n}\) as an expression intended to impose a limit to population growth. This criticism applies wherever the phenomena for which such an equation might be a model are not continuous and those are the very systems with which Chaos theorists are most concerned.

APL: the Dragon Tamer¶

What are the characteristics of the logistic equation if it is constrained to work at given resolutions? I have used some APL functions to calculate the results for various maximum populations, p being set to some power of 2. Thus we can compare and contrast the effects of doubling the ‘fineness of the grid’ on the sets of iterated results. A consequence of limiting the population to p is that there can only be ⍳p starting values and a consequence of taking the floor of the division is that every result must also be taken from the maximum set of integers ⍳p. The set of all possible result vectors can thus be shown as a p p matrix. Three examples are given below.

      ∇try[⎕]∇
    ∇ r←try p;⎕io
[1]   ⎕io←1
[2]   r←(p,1)⍴⍳p
[3]   l1:→(p=¯1↑⍴r)↑0
[4]   r←r,⌊(÷p)×4×r[;¯1↑⍴r]×p-r[;¯1↑⍴r]
[5]   →l1
    ∇

      try 8
1 3 7 3 7 3 7 3
2 6 6 6 6 6 6 6
3 7 3 7 3 7 3 7
4 8 0 0 0 0 0 0
5 7 3 7 3 7 3 7
6 6 6 6 6 6 6 6
7 3 7 3 7 3 7 3
8 0 0 0 0 0 0 0
      try 16
 1  3  9 15  3  9 15  3  9 15  3  9 15  3  9 15
 2  7 15  3  9 15  3  9 15  3  9 15  3  9 15  3
 3  9 15  3  9 15  3  9 15  3  9 15  3  9 15  3
 4 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12
 5 13  9 15  3  9 15  3  9 15  3  9 15  3  9 15
 6 15  3  9 15  3  9 15  3  9 15  3  9 15  3  9
 7 15  3  9 15  3  9 15  3  9 15  3  9 15  3  9
 8 16  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 9 15  3  9 15  3  9 15  3  9 15  3  9 15  3  9
10 15  3  9 15  3  9 15  3  9 15  3  9 15  3  9
11 13  9 15  3  9 15  3  9 15  3  9 15  3  9 15
12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12
13  9 15  3  9 15  3  9 15  3  9 15  3  9 15  3
14  7 15  3  9 15  3  9 15  3  9 15  3  9 15  3
15  3  9 15  3  9 15  3  9 15  3  9 15  3  9 15
16  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
      32 16↑try 32
 1  3 10 27 16 32  0  0  0  0  0  0  0  0  0  0
 2  7 21 28 14 31  3 10 27 16 32  0  0  0  0  0
 3 10 27 16 32  0  0  0  0  0  0  0  0  0  0  0
 4 14 31  3 10 27 16 32  0  0  0  0  0  0  0  0
 5 16 32  0  0  0  0  0  0  0  0  0  0  0  0  0
 6 19 30  7 21 28 14 31  3 10 27 16 32  0  0  0
 7 21 28 14 31  3 10 27 16 32  0  0  0  0  0  0
 8 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24
 9 25 21 28 14 31  3 10 27 16 32  0  0  0  0  0
10 27 16 32  0  0  0  0  0  0  0  0  0  0  0  0
11 28 14 31  3 10 27 16 32  0  0  0  0  0  0  0
12 30  7 21 28 14 31  3 10 27 16 32  0  0  0  0
13 30  7 21 28 14 31  3 10 27 16 32  0  0  0  0
14 31  3 10 27 16 32  0  0  0  0  0  0  0  0  0
15 31  3 10 27 16 32  0  0  0  0  0  0  0  0  0
16 32  0  0  0  0  0  0  0  0  0  0  0  0  0  0
17 31  3 10 27 16 32  0  0  0  0  0  0  0  0  0
18 31  3 10 27 16 32  0  0  0  0  0  0  0  0  0
19 30  7 21 28 14 31  3 10 27 16 32  0  0  0  0
20 30  7 21 28 14 31  3 10 27 16 32  0  0  0  0
21 28 14 31  3 10 27 16 32  0  0  0  0  0  0  0
22 27 16 32  0  0  0  0  0  0  0  0  0  0  0  0
23 25 21 28 14 31  3 10 27 16 32  0  0  0  0  0
24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24
25 21 28 14 31  3 10 27 16 32  0  0  0  0  0  0
26 19 30  7 21 28 14 31  3 10 27 16 32  0  0  0
27 16 32  0  0  0  0  0  0  0  0  0  0  0  0  0
28 14 31  3 10 27 16 32  0  0  0  0  0  0  0  0
29 10 27 16 32  0  0  0  0  0  0  0  0  0  0  0
30  7 21 28 14 31  3 10 27 16 32  0  0  0  0  0
31  3 10 27 16 32  0  0  0  0  0  0  0  0  0  0
32  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0

As seen the set of result vectors for ⍳p starting values falls into three categories:

The iteration reaches the value p. All values after that are zero.
The iteration reaches some other value which stays the same for every subsequent iteration.
The iteration reaches one of a small subset of values and subsequently cycles through that set. Different starting values may enter the same cycle but at different points within it. This could be expressed as ‘having the same periodicity but a different phase’. Other starting values may enter a different cycle constituting a different periodicity, again with phase variation. Different cycles, of course, have no values in common.

I have explored populations up to 2 power 12 and am intrigued to find that the mix of categories 1,2,3 above is very variable and that there are many more values which appear only once in an iteration than the small subset which participate in the cycles and that the same cycle is common to many starting values.

As larger tables are unwieldy and take too long to calculate in full, I have provided a function to analyse the distribution of values, eliminating rows which lead only to zero or to a steady state and terminating the calculation when all results are merely permutations of the same small subset of numbers. To facilitate printing I show only the first column and the last 10 or 16 columns of the larger matrices.

      ∇limited[⎕]∇
    ∇ r←limited p;⎕io;v1;v2
[1]   ⍝ analyse results of each iteration &
[2]   ⍝ stop when only recycling a subset of ⍳p
[3]   ⎕io←1
[4]   ⍝ set r to all possible initial values
[5]   r←(p,1)⍴v1←⍳p
[6]   ⍝ v2 is next generation
[7]   l1:v2←⌊(4×v1×p-v1)÷p
[8]   ⍝ exclude rows where values will never again change
[9]   v2←(b←~(v1=v2)∨0=v2)/v2
[10]  r←(b⌿r),v2
[11]  ⍝ keep only first occurrence of each value
[12]  v1←(b←(v2⍳v2)=⍳⍴v2)/v2
[13]  r←b⌿r
[14]  →((⍴v1)=⍴v2)/0
[15]  →l1
    ∇

      limited 64
1  3 11 36 63  3 11 36 63  3
 3 11 36 63  3 11 36 63  3 11
 4 15 45 53 36 63  3 11 36 63
 8 28 63  3 11 36 63  3 11 36

      limited 128
  1   3  11  40 110  61 127   3  11  40 110
  2   7  26  82 117  40 110  61 127   3  11
  3  11  40 110  61 127   3  11  40 110  61
  7  26  82 117  40 110  61 127   3  11  40
  8  30  91 105  75 124  15  52 123  19  64
  9  33  97  93 101  85 114  49 120  30  91
 10  36 103  80 120  30  91 105  75 124  15
 11  40 110  61 127   3  11  40 110  61 127
 14  49 120  30  91 105  75 124  15  52 123
 24  78 121  26  82 117  40 110  61 127   3
 30  91 105  75 124  15  52 123  19  64 128
 31  93 101  85 114  49 120  30  91 105  75
 33  97  93 101  85 114  49 120  30  91 105
 35 101  85 114  49 120  30  91 105  75 124
 36 103  80 120  30  91 105  75 124  15  52
 48 120  30  91 105  75 124  15  52 123  19

      matrix←limited 256
      ⍴matrix
4 28

      matrix[;1,12+⍳16]
  1 140 253  11  42 140 253  11  42 140 253  11  42 140 253  11  42
  2 241  56 175 221 120 255   3  11  42 140 253  11  42 140 253  11
  3 253  11  42 140 253  11  42 140 253  11  42 140 253  11  42 140
  4  11  42 140 253  11  42 140 253  11  42 140 253  11  42 140 253

      matrix←limited 512

      ⍴matrix
14 29

      matrix[;1,13+⍳16]
  1 261 511   3  11  43 157 435 261 511   3  11  43 157 435 261 511
  2 211 496  62 217 500  46 167 450 217 500  46 167 450 217 500  46
  3 511   3  11  43 157 435 261 511   3  11  43 157 435 261 511   3
  4 144 414 316 483 109 343 452 211 496  62 217 500  46 167 450 217
  5  11  43 157 435 261 511   3  11  43 157 435 261 511   3  11  43
  6 157 435 261 511   3  11  43 157 435 261 511   3  11  43 157 435
  7 496  62 217 500  46 167 450 217 500  46 167 450 217 500  46 167
  8   3  11  43 157 435 261 511   3  11  43 157 435 261 511   3  11
 10 500  46 167 450 217 500  46 167 450 217 500  46 167 450 217 500
 12 450 217 500  46 167 450 217 500  46 167 450 217 500  46 167 450
 17 435 261 511   3  11  43 157 435 261 511   3  11  43 157 435 261
 19  43 157 435 261 511   3  11  43 157 435 261 511   3  11  43 157
 49 177 463 177 463 177 463 177 463 177 463 177 463 177 463 177 463
177 463 177 463 177 463 177 463 177 463 177 463 177 463 177 463 177

      matrix←limited 4096
      ⍴matrix
20 117
            matrix[;1,101+⍳16]
   1 3301 2562 3837  970 2961 3281 2611 3786 1146 3301 2562 3837  970 2961 3281 2611
   2 3786 1146 3301 2562 3837  970 2961 3281 2611 3786 1146 3301 2562 3837  970 2961
   3 2562 3837  970 2961 3281 2611 3786 1146 3301 2562 3837  970 2961 3281 2611 3786
   4 1146 3301 2562 3837  970 2961 3281 2611 3786 1146 3301 2562 3837  970 2961 3281
   5 2961 3281 2611 3786 1146 3301 2562 3837  970 2961 3281 2611 3786 1146 3301 2562
   9 3837  970 2961 3281 2611 3786 1146 3301 2562 3837  970 2961 3281 2611 3786 1146
  14 3767 1210 3410 2284 4041  217  822 2628 3767 1210 3410 2284 4041  217  822 2628
  19 3281 2611 3786 1146 3301 2562 3837  970 2961 3281 2611 3786 1146 3301 2562 3837
  25  970 2961 3281 2611 3786 1146 3301 2562 3837  970 2961 3281 2611 3786 1146 3301
  29 2611 3786 1146 3301 2562 3837  970 2961 3281 2611 3786 1146 3301 2562 3837  970
  55 1210 3410 2284 4041  217  822 2628 3767 1210 3410 2284 4041  217  822 2628 3767
  84 4041  217  822 2628 3767 1210 3410 2284 4041  217  822 2628 3767 1210 3410 2284
 145  217  822 2628 3767 1210 3410 2284 4041  217  822 2628 3767 1210 3410 2284 4041
 203 2503 3893  771 2503 3893  771 2503 3893  771 2503 3893  771 2503 3893  771 2503
 217 3410 2284 4041  217  822 2628 3767 1210 3410 2284 4041  217  822 2628 3767 1210
 519 2628 3767 1210 3410 2284 4041  217  822 2628 3767 1210 3410 2284 4041  217  822
 559  822 2628 3767 1210 3410 2284 4041  217  822 2628 3767 1210 3410 2284 4041  217
 771 3893  771 2503 3893  771 2503 3893  771 2503 3893  771 2503 3893  771 2503 3893
 822 2284 4041  217  822 2628 3767 1210 3410 2284 4041  217  822 2628 3767 1210 3410
1593  771 2503 3893  771 2503 3893  771 2503 3893  771 2503 3893  771 2503 3893  771

I hope that this analysis demonstrates why an iteration using floating point calculations will seem chaotic because the equation is applied to successively finer breakdowns of the grid but that where the grid is fixed (however large p may be) most starting values will lead in not too many iterations to a few cycles with periods short by comparison with ⍴⍳p. From these experiments we might come to expect that different populations will have different propensities to reach a small subset of values and that any assumption of equiprobability of outcomes would be grossly mistaken. Some values can only be a starting value as they can never be reached from any other value.