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 を見つけるにはどうすればよいですか ??? 式で：

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

ありがとう

4

6 に答える 6

13

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

再帰関係には複数の固定点があります (1 つは 2 で、もう 1 つは 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

(非線形の) 繰り返しシーケンスには、 と の 3 つの固定点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)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

無限の精度で計算すると、次のようになり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 に答える
2

初期値を 53 ビットの精度で計算し、その丸められた値を高精度の mpf 変数に割り当てます。u0=mpf(3)/mpf(2) および u1=mpf(5)/mpf(3) を使用する必要があります。さらに数回の反復で 2 近くにとどまりますが、最終的には 2000 に収束します。これは丸め誤差によるものです。1 つの代替方法は、分数で計算することです。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 に答える
1

漸化式（初期値u_0 = 3/2、u_1 = 5/3）の正確な解は、次のように簡単に検証できます。

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

あなたが見ている問題は、解決策は

lim_{n -> oo} u_n = 2,

この制限は、漸化式の反発する不動点です。つまり、u_ {n-1}、u {n-2}の正しい値から逸脱すると、一部のnについて、さらに値が正しい制限から逸脱することになります。したがって、漸化式の実装がすべての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 に答える
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 に答える