Python计算错误

Question

我正在使用一个叫做mpmath的API来计算以下的总和。

我们考虑一个序列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 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

我的代码是：

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

你先计算出初始值，精确到53位，然后把这个四舍五入后的值赋给高精度的mpf变量。你应该用u0=mpf(3)/mpf(2)和u1=mpf(5)/mpf(3)。在接下来的几次迭代中，你的结果会一直接近2，但最后还是会收敛到2000。这是因为有舍入误差。一个替代的方法是用分数来计算。我使用了gmpy，下面的代码可以收敛到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 2

你的这个（非线性）递归序列有三个固定点：1、2和2000。其中，1和2之间的距离比2000要近，这通常意味着这些固定点不太稳定，因为它们“几乎”是双重根。

为了让结果不那么早就发散，你需要做一些数学运算。设v(n)为一个辅助序列：

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

接下来，这个关系是成立的：

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))

然后你可以简单地计算v(n)，并从u(n) = v(n)/(1+2^n)中推导出u(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

最后的结果是：

[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'),
...

请注意，这个方法最终还是会发散。要真正让它收敛，你需要用任意精度来计算v(n)。不过现在这变得简单多了，因为所有的值都是整数。

Answer 3

使用 decimal 模块，你会发现这个系列的解也会收敛到 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 的解是一个不稳定的固定点。而在 2000 的解是一个 吸引型固定点。

当结果接近 2 时，如果因为四舍五入的原因稍微超过了 2，这个差距会不断放大，最终会达到 2000。