Python calculation error

Question

Im using the mpmath API to calculate the next sum

Consider the series u0, u1, u2, defined as follows:

u0 = 3/2 = 1,5 u1 = 5/3 = 1,6666666… un+1 = 2003 - 6002/un + 4000/un un-1 

The series converges at 2, but with the rounding problem it seems to converge at 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,99676781 1,992248062
 8 1999,999997 1,996108949
 9 2000.000000 1.998050682
 10 2000.000000 1.999024390

My code is:

 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 execute some other functions (fdiv ...) or change the accuracy: the same bad result

What is 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

Using the decimal module, you can see that the series also has a solution converging in 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 several fixed points (one in 2, and the other in 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 point 2 is an unstable fixed point. Attractive fixed point located in 2000.

The convergence is very close to two, and when rounding causes the value to slightly exceed two, this difference amplifies again and again until it hits 2000.

Answer 2

Your (non-linear) repetition sequence has three fixed points: 1 , 2 and 2000 . Values ​​1 and 2 are close to each other compared to 2000, which usually indicates unstable fixed points, since they are “almost” double roots.

You need a little math to diverge less early. Let v(n) be a side sequence:

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

The following is 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)) 

Then you can simply calculate 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'), ... 

Please note that this will still be rejected . To really converge, you need to calculate v(n) with arbitrary precision. But now it is much simpler, since all values ​​are integers .

Answer 3

You calculate your initial values ​​to 53 bits of precision, and then assign this rounded value to the high-precision variable mpf. You should use u0 = mpf (3) / mpf (2) and u1 = mpf (5) / mpf (3). You will stay close to 2 for a few extra interactions, but you will end the rapprochement anyway in 2000. This is due to a rounding error. One option is fractional computation. 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)) 

Answer 4

If you calculate with infinite precision, 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)) 

Exit

 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 5

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

and 0 = MPF (3) / MPF (2)

u1 = MPF (5) / MPF (3)

and the value converges 16 steps to the correct value of 2.0 before diverging again (see below).

Thus, even with a good python library for arbitrary precision floating point arithmetic and some basic operations, the result can become completely false and this is not an algorithmic, mathematical or repetitive problem, as I sometimes read.

So you need to remain vigilant and critic !!! (Im very afraid of mpmath.lerchphi (z, s, a); -)

2 1.800000000000000000000000000000000000000000000000000 3 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.54167212 91457644804673979070312813731252347786 21 1999.9985406086893 66669273522363692463645090555294 22 1999.9999985406079725746311606572627439743947878652

Answer 6

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

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

The problem you see is that although the solution is such that

 lim_{n -> oo} u_n = 2, 

this limit is the repulsive fixed point of your recurrence relation. That is, any deviation from the correct values ​​u_ {n-1}, u {n-2}, for some n, will lead to further values ​​deviating from the correct limit. Therefore, if your implementation of the recurrence relation correctly does not represent any u_n value exactly, you can expect it to have a possible deviation from the correct limit, converging to the incorrect value 2000, which is simply the only attractive fixed point of your repetition relation,

---

(*) Indeed, u_n = (2 ^ (n + 1) + 1) / (2 ^ n + 1) is a 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 indicated above, where C is an arbitrary constant. If I were not mistaken in finding the roots of the characteristic polynomial, then this will have a set of fixed points {1, 2, C - 3} \ {0}. Limit 2 can be either a repelling fixed point or an attractive fixed point, depending on the value of CEg, for C = 2003 the set of fixed points is {1, 2, 2000}, where 2 is a repellent, whereas for C = 3 fixed points are {1, 2}, where 2 is an attractor.