Python computing error

Question

I’m using the API mpmath to compute the following sum

Let us consider the serie u0, u1, u2 defined by:

u0 = 3/2 = 1,5

u1 = 5/3 = 1,6666666…

un+1 = 2003 - 6002/un + 4000/un un-1

The serie converges on 2, but with rounding problem it seems to converge on 2000.

n   Calculated value    Rounded off exact value

2   1,800001            1,800000000
3   1,890000            1,888888889
4   3,116924            1,941176471
5   756,3870306         1,969696970
6   1996,761549         1,984615385
7   1999,996781         1,992248062
8   1999,999997         1,996108949
9   2000,000000         1,998050682
10  2000,000000         1,999024390

My code :

from mpmath import *
mp.dps = 50
u0=mpf(3/2.0)
u1=mpf(5/3.0)
u=[]
u.append(u0)
u.append(u1)
for i in range (2,11):
    un1=(2003-6002/u[i-1]+(mpf(4000)/mpf((u[i-1]*u[i-2]))))
    u.append(un1)
print u

my bad results :

[mpf('1.5'),   
 mpf('1.6666666666666667406815349750104360282421112060546875'),     
 mpf('1.8000000000000888711326751945268011597589466120961647'),  
 mpf('1.8888888889876302386905492787148253684796100079942617'),  
 mpf('1.9411765751351638992775070422559330255517747908588059'),    
 mpf('1.9698046831709839591526211645628191427874374792786951'),      
 mpf('2.093979191783975876606205176530675127058752077926479'),   
 mpf('106.44733511712489354422046139349654833300787666477228'),     
 mpf('1964.5606972399290690749220686397494349501387742896911'),   
 mpf('1999.9639916238009625032390578545797067344576357100626'),   
 mpf('1999.9999640260895343960004614025893194430187653900418')]

I tried to perform with some others functions (fdiv…) or to change the precision: same bad result

What’s wrong with this code ?

Question: How to change my code to find the value 2.0 ??? with the formula :

un+1 = 2003 - 6002/un + 4000/un un-1

thanks

Answer 1

The exact solution to your recurrence relation (with initial values u_0 = 3/2, u_1 = 5/3) is easily verified to be

u_n = (2^(n+1) + 1) / (2^n + 1).    (*)

The problem you're seeing is that although the solution is such that

lim_{n -> oo} u_n = 2,

this limit is a repelling fixed point of your recurrence relation. That is, any departure from the correct values of u_{n-1}, u{n-2}, for some n, will result in further values diverging from the correct limit. Consequently, unless your implementation of the recurrence relation correctly represents every u_n value exactly, it can be expected to exhibit eventual divergence from the correct limit, converging to the incorrect value of 2000 that just happens to be the only attracting fixed point of your recurrence relation.

---

(*) In fact, u_n = (2^(n+1) + 1) / (2^n + 1) is the solution to any recurrence relation of the form

u_n = C + (7 - 3C)/u_{n-1} + (2C - 6)/(u_{n-1} u_{n-2})

with the same initial values as given above, where C is an arbitrary constant. If I haven't made a mistake finding the roots of the characteristic polynomial, this will have the set of fixed points {1, 2, C - 3}\{0}. The limit 2 can be either a repelling fixed point or an attracting fixed point, depending on the value of C. E.g., for C = 2003 the set of fixed points is {1, 2, 2000} with 2 being a repellor, whereas for C = 3 the fixed points are {1, 2} with 2 being an attractor.

Answer 2

Well, as casevh said, I just added the mpf function in first initials terms in my code :

u0=mpf(3)/mpf(2)

u1=mpf(5)/mpf(3)

and the value converge for 16 steps to the correct value 2.0 before diverged again (see below).

So, even with a good python library for arbitrary-precision floating-point arithmetic and some basics operations the result can become totally false and it is not algorithmic, mathematical or recurrence problem as I read sometimes.

So it is necessary to remain watchful and critic !!! ( I’m very afraid about the mpmath.lerchphi(z, s, a) function ;-)

2 1.8000000000000000000000000000000000000000000000022 3 1.8888888888888888888888888888888888888888888913205 4 1.9411764705882352941176470588235294117647084569125 5 1.9696969696969696969696969696969696969723495083846 6 1.9846153846153846153846153846153846180779422496889 7 1.992248062015503875968992248062018218070968279944 8 1.9961089494163424124513618677070049064461141667961 9 1.998050682261208576998050684991268132991329645551 10 1.9990243902439024390243929766241359876402781522945 11 1.9995119570522205954151303455889283862002420414092 12 1.9997559189650964147435086295745928366095548127257 13 1.9998779445868451615169464386495752584786229236677 14 1.9999389685715481608370784691478769380770569091713 15 1.9999694860884747554701272066241108169217231319376 16 1.9999874767910784720428384947047783821702386000249 17 2.0027277350948824117795762659330557916802871427763 18 4.7316350177463946015607576536159982430500337286276 19 1156.6278675611076227796014310764287933259776352198 20 1998.5416721291457644804673979070312813731252347786 21 1999.998540608689366669273522363692463645090555294 22 1999.9999985406079725746311606572627439743947878652

Answer 3

Your (non-linear) recurrence sequence has three fixed points: 1, 2 and 2000. The values 1 and 2 are close to each other compared to 2000, which is usually an indication of unstable fixed points because they are "almost" double roots.

You need to do some maths in order to diverge less early. Let v(n) be a side sequence:

v(n) = (1+2^n)u(n)

The following holds true:

v(n+1) = (1+2^(n+1)) * (2003v(n)v(n-1) - 6002(1+2^n)v(n-1) + 4000(1+2^n)(1+2^n-1)) / (v(n)v(n-1))

You can then simply compute v(n) and deduce u(n) from u(n) = v(n)/(1+2^n):

#!/usr/bin/env python

from mpmath import *
mp.dps = 50
v0 = mpf(3)
v1 = mpf(5)
v=[]
v.append(v0)
v.append(v1)

u=[]
u.append(v[0]/2)
u.append(v[1]/3)

for i in range (2,25):
    vn1 = (1+2**i) * (2003*v[i-1]*v[i-2] \
                     - 6002*(1+2**(i-1))*v[i-2] \
                     + 4000*(1+2**(i-1))*(1+2**(i-2))) \
                   / (v[i-1]*v[i-2])
    v.append(vn1)
    u.append(vn1/(1+2**i))

print u

And the result:

[mpf('1.5'),
mpf('1.6666666666666666666666666666666666666666666666666676'),
mpf('1.8000000000000000000000000000000000000000000000000005'),
mpf('1.8888888888888888888888888888888888888888888888888892'),
mpf('1.9411764705882352941176470588235294117647058823529413'),
mpf('1.969696969696969696969696969696969696969696969696969'),
mpf('1.9846153846153846153846153846153846153846153846153847'),
mpf('1.992248062015503875968992248062015503875968992248062'),
mpf('1.9961089494163424124513618677042801556420233463035019'),
mpf('1.9980506822612085769980506822612085769980506822612089'),
mpf('1.9990243902439024390243902439024390243902439024390251'),
mpf('1.9995119570522205954123962908735968765251342118106393'),
mpf('1.99975591896509641200878691725652916768367097876495'),
mpf('1.9998779445868424264616135725619431221774685707311133'),
mpf('1.9999389685688129386634116570033567287152883735123589'),
mpf('1.9999694833531691537733833806341359211449845890933504'),
mpf('1.9999847414437645909944001098616048949448403192089965'),
mpf('1.9999923706636759668276456631037666033431751771913355'),
...

Note that this will still diverge eventually. In order to really converge, you need to compute v(n) with arbitrary precision. But this is now a lot easier since all the values are integers.

Answer 4

Using the decimal module, you can see the series also has a solution converging at 2000:

from decimal import Decimal, getcontext
getcontext().prec = 100

u0=Decimal(3) / Decimal(2)
u1=Decimal(5) / Decimal(3)
u=[u0, u1]
for i in range(100):
    un1 = 2003 - 6002/u[-1] + 4000/(u[-1]*u[-2])
    u.append(un1)
    print un1

The recurrence relation has multiple fixed points (one at 2 and the other at 2000):

>>> u = [Decimal(2), Decimal(2)]
>>> 2003 - 6002/u[-1] + 4000/(u[-1]*u[-2])
Decimal('2')

>>> u = [Decimal(2000), Decimal(2000)]
>>> 2003 - 6002/u[-1] + 4000/(u[-1]*u[-2])
Decimal('2000.000')

The solution at 2 is an unstable fixed-point. The attractive fixed-point is at 2000.

The convergence gets very close to two and when the round-off causes the value to slightly exceed two, that difference gets amplified again and again until hitting 2000.

Answer 5

If you compute with infinite precision then you get 2 otherwise you get 2000:

import itertools
from fractions import Fraction

def series(u0=Fraction(3, 2), u1=Fraction(5, 3)):
    yield u0
    yield u1
    while u0 != u1:
        un = 2003 - 6002/u1 + 4000/(u1*u0)
        yield un
        u1, u0 = un, u1

for i, u in enumerate(itertools.islice(series(), 100)):
    err = (2-u)/2 # relative error
    print("%d\t%.2g" % (i, err))

Output

0   0.25
1   0.17
2   0.1
3   0.056
4   0.029
5   0.015
6   0.0077
7   0.0039
8   0.0019
9   0.00097
10  0.00049
11  0.00024
12  0.00012
13  6.1e-05
14  3.1e-05
15  1.5e-05
16  7.6e-06
17  3.8e-06
18  1.9e-06
19  9.5e-07
20  4.8e-07
21  2.4e-07
22  1.2e-07
23  6e-08
24  3e-08
25  1.5e-08
26  7.5e-09
27  3.7e-09
28  1.9e-09
29  9.3e-10
30  4.7e-10
31  2.3e-10
32  1.2e-10
33  5.8e-11
34  2.9e-11
35  1.5e-11
36  7.3e-12
37  3.6e-12
38  1.8e-12
39  9.1e-13
40  4.5e-13
41  2.3e-13
42  1.1e-13
43  5.7e-14
44  2.8e-14
45  1.4e-14
46  7.1e-15
47  3.6e-15
48  1.8e-15
49  8.9e-16
50  4.4e-16
51  2.2e-16
52  1.1e-16
53  5.6e-17
54  2.8e-17
55  1.4e-17
56  6.9e-18
57  3.5e-18
58  1.7e-18
59  8.7e-19
60  4.3e-19
61  2.2e-19
62  1.1e-19
63  5.4e-20
64  2.7e-20
65  1.4e-20
66  6.8e-21
67  3.4e-21
68  1.7e-21
69  8.5e-22
70  4.2e-22
71  2.1e-22
72  1.1e-22
73  5.3e-23
74  2.6e-23
75  1.3e-23
76  6.6e-24
77  3.3e-24
78  1.7e-24
79  8.3e-25
80  4.1e-25
81  2.1e-25
82  1e-25
83  5.2e-26
84  2.6e-26
85  1.3e-26
86  6.5e-27
87  3.2e-27
88  1.6e-27
89  8.1e-28
90  4e-28
91  2e-28
92  1e-28
93  5e-29
94  2.5e-29
95  1.3e-29
96  6.3e-30
97  3.2e-30
98  1.6e-30
99  7.9e-31

Answer 6

You calculate your initial values to 53-bits of precision and then assign that rounded value to the high-precision mpf variable. You should use u0=mpf(3)/mpf(2) and u1=mpf(5)/mpf(3). You'll stay close to 2 for a few more interations, but you'll still end up converging at 2000. This is due to rounding error. One alternative is to compute with fractions. I used gmpy and the following code converges to 2.

from __future__ import print_function
import gmpy

u = [gmpy.mpq(3,2), gmpy.mpq(5,3)]
for i in range(2,300):
    temp = (2003 - 6002/u[-1] + 4000/(u[-1]*u[-2]))
    u.append(temp)

for i in u: print(gmpy.mpf(i,300))