Ошибка вычисления 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 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

Используя десятичный модуль, вы можете увидеть, что ряд также имеет решение, сходящееся к 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.

Сходимость становится очень близкой к двум, и когда округление приводит к тому, что значение немного превышает два, эта разница снова и снова усиливается, пока не достигнет 2000.

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) из 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) с произвольной точностью. Но теперь это намного проще, так как все значения целые числа .

Answer 3

Если вы вычислите с бесконечной точностью, тогда вы получите 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

Answer 4

Вы вычисляете свои начальные значения с точностью до 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 5

Точное решение вашего рекуррентного отношения (с начальными значениями 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 является аттрактором.

Answer 6

Ну, как сказал casevh, я просто добавил в свой код функцию mpf в терминах первых инициалов:

u0 = mpf (3) / mpf (2)

u1 = mpf (5) / mpf (3)

и значение сходится за 16 шагов к правильному значению 2.0, прежде чем снова расходится (см. Ниже).

Таким образом, даже с хорошей библиотекой python для произвольноТочная арифметика с плавающей точкой и некоторые базовые операции: результат может стать полностью ложным, и это не алгоритмическая, не математическая или не повторяющаяся проблема, как я иногда читаю.

Так что нужно быть бдительным и критичным !!!(Я очень боюсь о mpmath.lerchphi (г, з, а) функции; -)

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 111,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