Python计算错误

Question

我正在使用 API mpmath 计算以下总和

让我们考虑由以下定义的系列 u0、u1、u2:

u0 = 3/2 = 1,5

u1 = 5/3 = 1,6666666…

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

序列收敛于 2，但由于四舍五入问题似乎收敛于 2000。

n   Calculated value    Rounded off exact value2   1,800001            1,8000000003   1,890000            1,8888888894   3,116924            1,9411764715   756,3870306         1,9696969706   1996,761549         1,9846153857   1999,996781         1,9922480628   1999,999997         1,9961089499   2000,000000         1,99805068210  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

我的坏结果:

[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')]

我尝试执行其他一些函数(fdiv ...)或更改精度:同样的错误结果

这段代码有什么问题？

问题:如何更改我的代码以找到值 2.0 ???使用公式:

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

谢谢

Answer 1

使用十进制模块，可以看到级数也有收敛于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

递推关系有多个不动点(一个在2，另一个在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')

2 处的解是一个不稳定的不动点。 attractive fixed-point为 2000。

收敛非常接近 2，当舍入导致值略微超过 2 时，差异会一次又一次地放大，直到达到 2000。