Series Stability
Consider the following two series:
Series A
xn+1=4xn(1-xn)
Series B
xn+1=2xn(1-xn)
If we start each series with x1=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: 1101002 and 11001002. 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