8

我正在使用 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 计算值 四舍五入的精确值

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 ??? 用公式：

联合国+1 = 2003 - 6002/联合国 + 4000/联合国联合国-1

谢谢

4

6 回答 6

13

使用十进制模块，您可以看到该系列也有一个收敛于 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。

于 2012-01-02T09:11:16.037 回答
3

您的（非线性）循环序列具有三个固定点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)任意精度进行计算。但这现在容易多了，因为所有的值都是整数。

于 2012-01-02T18:31:08.167 回答
2

您将初始值计算为 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))
于 2012-01-02T10:01:16.733 回答
2

如果你以无限精度计算，那么你会得到，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
于 2012-01-02T10:25:27.260 回答
1

好吧，正如 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

于 2012-01-03T07:21:16.193 回答
1

您的递归关系的精确解（初始值 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 的不正确值，这恰好是您的递推关系的唯一吸引固定点.

(*) 实际上，u_n = (2^(n+1) + 1) / (2^n + 1) 是形式的任何递归关系的解

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 是吸引子。

于 2012-01-05T05:30:58.553 回答