I am working with the recurrence "xn+1 = (13/3)xn - (4/3)xn-1". I am trying to write a Python script to print the first 50 values of xn with the given values of x0 = 1 and x1 = 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 xn = 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.
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