# Series Stability

Consider the following two series:

### Series A

x_{n+1}=4x_{n}(1-x_{n})

### Series B

x_{n+1}=2x_{n}(1-x_{n})

If we start each series with x_{1}=0.52, what is the
52nd term?

This sounds like a fairly straight-forward exercise.

seriesA.py

#!/usr/bin/python3 x=0.52 for i in range (1,53): print(i,x) x=4*x*(1-x)

seriesB.py

#!/usr/bin/python3 x=0.52 for i in range (1,53): print(i,x) x=2*x*(1-x)

Having made the above files executable with

chmod +x seriesA.py seriesB.py

we can run them as "`./seriesA.py`

" and
"`./seriesB.py`

". The output should be something like:

Series A

1 0.52 2 0.9984 3 0.0063897600000001826 4 0.025395723868570322 5 0.09900312431104659 6 0.35680602275079215 [...] 50 0.34260639033414503 51 0.9009090065454101 52 0.3570878738830894

Series B

1 0.52 2 0.4992 3 0.49999872 4 0.4999999999967232 5 0.5 6 0.5 [...] 50 0.5 51 0.5 52 0.5

Job done?

## We have a problem

The correct answer for series B is indeed 0.5 (or, arguably, 0.49999..., but the number of 9's doubles on each iteration, so it is really quite close to 0.5).

The correct answer for series A is approximately 0.0001, and certainly nowhere near 0.357.

#!/usr/bin/python3 import decimal decimal.getcontext().prec=40 x=decimal.Decimal(52)/decimal.Decimal(100) for i in range (1,53): print(i,x) x=4*x*(1-x)

The above routine uses python's standard high-precision maths package, decimal. The line:

decimal.getcontext().prec=40

sets the number of decimal digits of precision to use, and the line

x=decimal.Decimal(52)/decimal.Decimal(100)

carefully sets `x`

to 0.52 to that precision. One should
not fall into the trap of trying

x=decimal.Decimal(0.52)

for that evaluates 0.52 as a standard double precision number, and
then converts that to a high-precision 40 decimal digit number. The problem
is that a standard double precision number is stored as a binary
fraction, and whilst 0.52 can be stored exactly with just two
decimal digits, as a binary fraction it recurs forever, just as one
seventh would expressed as a decimal number (or, indeed, a binary
one). The small integers 52 and 100 can be stored exactly using a
small number of binary digits: 110100_{2} and
1100100_{2}. The binary representation of 0.52 is

```
0.1000 0101 0001 1110 1011 1000 0101 0001 1110 1011 1000
```

where the overline indicates the infinitely-repeating sequence.

How do we know that our new code is correct? We can run it for different values of precision, and observe how the final term changes. Precision vs 52nd term are tabulated below:

14 0.24979965732246 15 0.00008050219746 16 0.00006594684318 17 0.00006629687479 18 0.00008699579715 20 0.00008543817532 25 0.00008542316204 30 0.00008542316220 35 0.00008542316220

And note that had we (incorrectly)
used `x=decimal.Decimal(0.52)`

we would conclude that the
answer is about 0.000947690, wrong by more than a factor of ten!

## A Warning

In general computers are marvellous for maths. But, just sometimes, they can be dangerously wrong. To me it is not immediately obvious that the two series presented here are so very different: both are simple quadratic expressions. Of course one can analyse their stability, but who would bother?

## Further Thoughts

Thirty decimal digits of precision were sufficient to reach the 52nd term of this series. What if we had wanted the hundredth, or the thousandth, term?

Standard double precision arithmetic is fast. How much slower is
this special high-precision arithmetic? Try removing the print
statement from the loop, and just printing `i`

and `x`

after the final iteration. Increase the loop
count to a million, and try timing the original double precision
code, the high precision code with a precision of 40 decimal places,
and the high precision code with a precision of 100 decimal
places.

It may be useful to have the Pi time itself, by using commands such as:

$ time ./seriesA.py