Recurrence Relations in Python

Question

I am working with the recurrence "x_n+1 = (13/3)x_n - (4/3)x_n-1". I am trying to write a Python script to print the first 50 values of x_n with the given values of x₀ = 1 and x₁ = 1/3. This is what my code currently looks like:

import math
def printRecurrence():
    x = [0]*51 #initialize list of x values
    x[0] = 1
    x[1] = 1/3
    for i in range(1, 51):
        x[i+1] = (13/3)*x[i] - (4/3)*x[i-1]
        print(x[i])

and the output I receive is:

0.3333333333333333
0.11111111111111094
0.03703703703703626
0.012345679012342514
0.004115226337435884
0.0013717421124321456
0.00045724737062478524
0.00015241578946454185
5.0805260179967644e-05
1.6935074827137338e-05
5.644977344304949e-06
1.8814687224716613e-06
6.263946716372672e-07
2.0575194713260943e-07
5.63988753916179e-08
-2.994080281313502e-08
-2.049419793790756e-07
-8.481608402251475e-07
-3.402107668470205e-06
-1.361158544307069e-05
-5.444739336201271e-05
-0.00021778992397796082
-0.0008711598127551465
-0.0034846392899683535
-0.013938557172856001
-0.05575422869575153
-0.22301691478444863
-0.8920676591382753
-3.5682706365532617
-14.2730825462131
-57.092330184852415
-228.36932073940963
-913.4772829576384
-3653.909131830553
-14615.63652732221
-58462.546109288836
-233850.18443715532
-935400.7377486213
-3741602.950994485
-14966411.80397794
-59865647.21591176
-239462588.86364704
-957850355.4545882
-3831401421.8183527
-15325605687.27341
-61302422749.09364
-245209690996.37457
-980838763985.4983
-3923355055941.993

which is only correct for the first 13 printed values. The proof I was provided was that x_n = 3^-n, which for the most part does not match the values of my script. Did I do something wrong in my calculations? I've been unable to see it.

Answer 1

This recurrence relation can be answered exactly through the use of the fractions module or with a variable level of precision using the decimal module, which demonstrates the very high level of precision required to find compute 50 iterations accurately.

Jean-François' point about floating point accumulation errors is correct. However the fraction module does not seem to be able to multiply multiple Fraction objects together, so all of the numeric values must be stated within the fraction object. Credit to him for the use of the correct module for this problem.

Exact answer

import math
import fractions

def printRecurrence():
    x = [0]*51 #initialize list of x values
    x[0] = fractions.Fraction(1,1)
    x[1] = fractions.Fraction(1,3)
    for i in range(1, 50):
        x[i+1] = fractions.Fraction(13*x[i]/3) - fractions.Fraction(4*x[i-1]/3)
        print(float(x[i+1]), 3**(-(i+1)), x[i+1]-(3)**(-(i+1)))

printRecurrence()

The print-out shows that the computation exactly matches the proof answer.

Inexact but instructive answer

The decimal module allows for a custom level of precision; to get to 50 iterations one needs a precision of 60 points or more.

import math
from decimal import *
getcontext().prec = 60

def printRecurrence():
    x = [0]*51 #initialize list of x values
    x[0] = Decimal(1)
    x[1] = Decimal(1)/Decimal(3)
    for i in range(1, 50):
        x[i+1] = (Decimal(13)/Decimal(3))*x[i] - (Decimal(4)/Decimal(3))*x[i-1]
        print(float(x[i+1]), 3**(-(i+1)), float(x[i+1]-Decimal(3)**(-(i+1))))

printRecurrence()

As with Jean-François' answer, I've printed the result, the computed value of 3**-n and the difference. One can play with the precision level getcontext().prec to see the effect on the result.

0.111111111111 0.111111111111 -1e-60
0.037037037037 0.037037037037 -4e-60
0.0123456790123 0.0123456790123 -1.57e-59
0.00411522633745 0.00411522633745 -6.243e-59
0.00137174211248 0.00137174211248 -2.4961e-58
0.000457247370828 0.000457247370828 -9.984e-58
0.000152415790276 0.000152415790276 -3.993577e-57
5.08052634253e-05 5.08052634253e-05 -1.59742978e-56
1.69350878084e-05 1.69350878084e-05 -6.38971879e-56
5.64502926948e-06 5.64502926948e-06 -2.5558875048e-55
1.88167642316e-06 1.88167642316e-06 -1.02235500146e-54
6.27225474386e-07 6.27225474386e-07 -4.08942000567e-54
2.09075158129e-07 2.09075158129e-07 -1.63576800226e-53
6.96917193763e-08 6.96917193763e-08 -6.54307200905e-53
2.32305731254e-08 2.32305731254e-08 -2.61722880362e-52
7.74352437514e-09 7.74352437514e-09 -1.04689152145e-51
2.58117479171e-09 2.58117479171e-09 -4.18756608579e-51
8.60391597238e-10 8.60391597238e-10 -1.67502643432e-50
2.86797199079e-10 2.86797199079e-10 -6.70010573727e-50
9.55990663597e-11 9.55990663597e-11 -2.68004229491e-49
3.18663554532e-11 3.18663554532e-11 -1.07201691796e-48
1.06221184844e-11 1.06221184844e-11 -4.28806767185e-48
3.54070616147e-12 3.54070616147e-12 -1.71522706874e-47
1.18023538716e-12 1.18023538716e-12 -6.86090827496e-47
3.93411795719e-13 3.93411795719e-13 -2.74436330998e-46
1.3113726524e-13 1.3113726524e-13 -1.09774532399e-45
4.37124217466e-14 4.37124217466e-14 -4.39098129597e-45
1.45708072489e-14 1.45708072489e-14 -1.75639251839e-44
4.85693574962e-15 4.85693574962e-15 -7.02557007356e-44
1.61897858321e-15 1.61897858321e-15 -2.81022802942e-43
5.39659527735e-16 5.39659527735e-16 -1.12409121177e-42
1.79886509245e-16 1.79886509245e-16 -4.49636484708e-42
5.99621697484e-17 5.99621697484e-17 -1.79854593883e-41
1.99873899161e-17 1.99873899161e-17 -7.19418375532e-41
6.66246330538e-18 6.66246330538e-18 -2.87767350213e-40
2.22082110179e-18 2.22082110179e-18 -1.15106940085e-39
7.40273700597e-19 7.40273700597e-19 -4.60427760341e-39
2.46757900199e-19 2.46757900199e-19 -1.84171104136e-38
8.22526333997e-20 8.22526333997e-20 -7.36684416545e-38
2.74175444666e-20 2.74175444666e-20 -2.94673766618e-37
9.13918148886e-21 9.13918148886e-21 -1.17869506647e-36
3.04639382962e-21 3.04639382962e-21 -4.71478026589e-36
1.01546460987e-21 1.01546460987e-21 -1.88591210636e-35
3.38488203291e-22 3.38488203291e-22 -7.54364842542e-35
1.12829401097e-22 1.12829401097e-22 -3.01745937017e-34
3.76098003645e-23 3.76098003657e-23 -1.20698374807e-33
1.25366001171e-23 1.25366001219e-23 -4.82793499227e-33
4.17886668798e-24 4.1788667073e-24 -1.93117399691e-32
1.39295549185e-24 1.3929555691e-24 -7.72469598763e-32

Answer 2

You are experiencing floating point accumulation error.

Since current result depends on previous result, and float type is just an approximation, the more recurrences you're running, the more the accumulation error increases.

Since you're dealing with rational number, I'd advise the fraction module.

import fractions
def printRecurrence():
    x = [0]*51 #initialize list of x values
    x[0] = fractions.Fraction(1,1)
    x[1] = fractions.Fraction(1,3)
    for i in range(1, 50):
        x[i+1] = fractions.Fraction(13/3)*x[i] - fractions.Fraction(4/3)*x[i-1]
        print(float(x[i]),3**(-i),x[i]-3**(-i))

printRecurrence()

I've printed the result, the computed value of 3**-n and the difference:

0.3333333333333333 0.3333333333333333 0.0
0.11111111111111109 0.1111111111111111 -1.3877787807814457e-17
0.037037037037036924 0.037037037037037035 -1.1102230246251565e-16
0.01234567901234521 0.012345679012345678 -4.683753385137379e-16
0.004115226337446681 0.00411522633744856 -1.878705524482882e-15
0.001371742112475337 0.0013717421124828531 -7.516123140538511e-15
0.0004572473707975519 0.0004572473708276177 -3.006584781486965e-14
0.00015241579015560884 0.00015241579027587258 -1.2026374362518466e-13
5.080526294423582e-05 5.080526342529086e-05 -4.810550354871108e-13
1.6935085884210096e-05 1.6935087808430286e-05 -1.9242201893822884e-12
5.645021572595982e-06 5.645029269476762e-06 -7.696880780399043e-12
1.8816456356357935e-06 1.8816764231589208e-06 -3.0787523127313645e-11
6.27102324293796e-07 6.272254743863069e-07 -1.2315009251094865e-10
2.0858255775872449e-07 2.0907515812876897e-07 -4.926003700444828e-10
6.772131789607814e-08 6.969171937625632e-08 -1.9704014801781827e-09
1.534896720470598e-08 2.3230573125418773e-08 -7.881605920712794e-09
-2.378289930771161e-08 7.743524375139592e-09 -3.15264236828512e-08
-1.2352451993969162e-07 2.581174791713197e-09 -1.2610569473140483e-07
-5.035623873283815e-07 8.603915972377324e-10 -5.044227789256192e-07
-2.0174043185033973e-06 2.8679719907924413e-10 -2.0176911157024764e-06
-8.070668863743547e-06 9.559906635974805e-11 -8.070764462809906e-06
-3.228302598488417e-05 3.186635545324935e-11 -3.228305785123962e-05
...

After 20 recurrences, the difference increases, but it's already much better.

EDIT: David answer made me realize that the fraction to float conversion spoils it all.

The fraction module only performs a division between 2 exact computed integers, but conversion to float ruins the precision.