Python计算错误

发布于 2024-12-24 22:00:03 字数 1717 浏览 2 评论 0原文

我正在使用 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 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

谢谢

原文

I’m using the API mpmath to compute the following sum

Let us consider the serie u0, u1, u2 defined by:

u0 = 3/2 = 1,5

u1 = 5/3 = 1,6666666…

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

The serie converges on 2, but with rounding problem it seems to converge on 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

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

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 perform with some others functions (fdiv…) or to change the precision: same bad result

What’s 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

分享到QQ

分享到微博

如果你对这篇内容有疑问，欢迎到本站社区发帖提问参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。

发布评论

需要登录才能够评论，你可以免费注册一个本站的账号。

み格子的夏天 2024-12-31 22:00:03

使用十进制模块，您可以看到级数在 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。

Using the decimal module, you can see the series also has a solution converging at 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 multiple fixed points (one at 2 and the other at 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 2 is an unstable fixed-point. The attractive fixed-point is at 2000.

The convergence gets very close to two and when the round-off causes the value to slightly exceed two, that difference gets amplified again and again until hitting 2000.

回复收藏 0 原文

薄荷港 2024-12-31 22:00:03

您的（非线性）递归序列具有三个固定点：1、2 和 2000。与 2000 相比，值 1 和 2 彼此接近，这通常表明不动点不稳定，因为它们“几乎”是双根。

你需要做一些数学计算，才能避免过早出现分歧。令 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) 来自 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

结果：

[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)。但现在这要容易得多，因为所有值都是整数。

Your (non-linear) recurrence sequence has three fixed points: 1, 2 and 2000. The values 1 and 2 are close to each other compared to 2000, which is usually an indication of unstable fixed points because they are "almost" double roots.

You need to do some maths in order to diverge less early. Let v(n) be a side sequence:

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

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

You can then simply compute 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'),
...

Note that this will still diverge eventually. In order to really converge, you need to compute v(n) with arbitrary precision. But this is now a lot easier since all the values are integers.

回复收藏 0 原文

执笏见 2024-12-31 22:00:03

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

You calculate your initial values to 53-bits of precision and then assign that rounded value to the high-precision mpf variable. You should use u0=mpf(3)/mpf(2) and u1=mpf(5)/mpf(3). You'll stay close to 2 for a few more interations, but you'll still end up converging at 2000. This is due to rounding error. One alternative is to compute with fractions. 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))

回复收藏 0 原文

遮了一弯 2024-12-31 22:00:03

如果您以无限精度计算，那么您将得到 2 否则您将得到 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))

输出

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

If you compute with infinite precision then 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))

Output

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

回复收藏 0 原文

柠檬心 2024-12-31 22:00:03

好吧，正如 casevh 所说，我只是在代码中以首字母缩写形式添加了 mpf 函数：

u0=mpf(3)/mpf(2)

u1=mpf(5)/mpf(3)

并且值收敛 16 个步骤再次偏离之前的正确值 2.0（见下文）。

因此，即使有一个很好的用于任意精度浮点运算和一些基本运算的Python库，结果也可能变得完全错误，并且它不是我有时读到的算法、数学或递归问题。

所以有必要保持警惕和批评！（我非常担心 mpmath.lerchphi(z, s, a) 函数；-)

2 1.8000000000000000000000000000000000000000000000022 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.5416721291457644804673979070312813731252347786 21 1999.998540608689366669273522363692463645090555294 22 1999.9999985406079725746311606572627439743947878652

回复收藏 0 原文

第几種人 2024-12-31 22:00:03

递归关系的精确解（初始值 u_0 = 3/2，u_1 = 5/3）很容易被验证为

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

您看到的问题是，尽管解是这样的，但

lim_{n -> oo} u_n = 2,

这个限制是令人排斥的递归关系的固定点。也就是说，对于某些 n，任何偏离 u_{n-1}、u{n-2} 正确值的情况都将导致值进一步偏离正确限制。因此，除非您的递归关系实现正确地表示每个 u_n 值准确，否则预计它最终会偏离正确的限制，收敛到恰好是唯一的错误值 2000吸引你的递归关系的固定点。

(*) In fact, u_n = (2^(n+1) + 1) / (2^n + 1) is the solution to any recurrence relation of the form

u_n = C + (7 - 3C)/u_{n-1} + (2C - 6)/(u_{n-1} u_{n-2})

具有与上面给出的相同的初始值，其中C是任意常数。如果我在查找特征多项式的根时没有犯错误，则这将具有一组固定点 {1, 2, C - 3}\{0}。极限 2 可以是排斥固定点或吸引固定点，具体取决于 CEg 的值，对于 C = 2003，固定点集是 {1, 2, 2000}，其中 2 是排斥点，而对于 C = 3 不动点是 {1, 2}，其中 2 是吸引子。

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

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

The problem you're seeing is that although the solution is such that

lim_{n -> oo} u_n = 2,

this limit is a repelling fixed point of your recurrence relation. That is, any departure from the correct values of u_{n-1}, u{n-2}, for some n, will result in further values diverging from the correct limit. Consequently, unless your implementation of the recurrence relation correctly represents every u_n value exactly, it can be expected to exhibit eventual divergence from the correct limit, converging to the incorrect value of 2000 that just happens to be the only attracting fixed point of your recurrence relation.

(*) In fact, u_n = (2^(n+1) + 1) / (2^n + 1) is the 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 given above, where C is an arbitrary constant. If I haven't made a mistake finding the roots of the characteristic polynomial, this will have the set of fixed points {1, 2, C - 3}\{0}. The limit 2 can be either a repelling fixed point or an attracting fixed point, depending on the value of C. E.g., for C = 2003 the set of fixed points is {1, 2, 2000} with 2 being a repellor, whereas for C = 3 the fixed points are {1, 2} with 2 being an attractor.

回复收藏 0 原文

~没有更多了~