Precision loss numpy - mpmath

Question

Я использую numpy и mpmath в моей программе Python. Я использую numpy, потому что он обеспечивает легкий доступ ко многим операциям с линейной алгеброй. Но поскольку решатель numpy для линейных уравнений не точен, я использую mpmath для более точных операций. После того, как я вычислить решение системы:

solution = mpmath.lu_solve(A,b) 

я хочу решение, как массив. Поэтому я использую

array = np.zeros(m) 

, а затем сделать петлю для установки значения:

for i in range(m): 
    array[i] = solution[i] 

или

for i in range(m): 
    array.put([i],solution[i]) 

но с обоими способами я получаю снова численные неустойчивости, как:

solution[0] = 12.375 
array[0] = 12.37500000000000177636 

Is t вот способ избежать этих ошибок?

Answer 1

numpy ndarrays имеют однородный тип. Когда вы делаете array, по умолчанию dtype будет некоторый тип поплавка, который не имеет столько точность, как вы хотите:

>>> array = np.zeros(3) 
>>> array 
array([ 0., 0., 0.]) 
>>> array.dtype 
dtype('float64') 

Вы можете обойти эту проблему с помощью dtype=object:

>>> mp.mp.prec = 65 
>>> mp.mpf("12.37500000000000177636") 
mpf('12.37500000000000177636') 
>>> array = np.zeros(3, dtype=object) 
>>> array[0] = 12.375 
>>> array[1] = mp.mpf("12.37500000000000177636") 
>>> array 
array([12.375, mpf('12.37500000000000177636'), 0], dtype=object) 

но обратите внимание, что при этом вы получаете значительное повышение производительности.

Answer 2

Для полноты и для таких людей, как я, которые наткнулись на этот вопрос, потому что линейный решатель numpy недостаточно точен (возможно, он способен обрабатывать только 64-битные номера), есть также sympy.

API немного похож на numpy, но требует несколько настроек время от времени.

In [104]: A = Matrix([ 
[17928014155669123106522437234449354737723367262236489360399559922258155650097260907649387867023242961198972825743674594974017771680414642705007756271459833, 13639120912900071306285490050678803027789658562874829601953000463099941578381997997439951418291413106684405816668933580642992992427754602071359317117391198, 2921704428390104906296926198429197524950528698980675801502622843572749765891484935643316840553487900050392953088680445022408396921815210925936936841894852], 
[14748352608418286197003564525494635018601621545162877692512866070970871867506554630832144013042243382377181384934564249544095078709598306314885920519882886, 2008780320611667023380867301336185953729900109553256489538663036530355388609791926150229595099056264556936500639831205368501493630132784265435798020329993, 6522019637107271075642013448499575736343559556365957230686263307525076970365959710482607736811185215265865108154015798779344536283754814484220624537361159], 
[ 5150176345214410132408539250659057272148578629292610140888144535962281110335200803482349370429701981480772441369390017612518504366585966665444365683628345, 1682449005116141683379601801172780644784704357790687066410713584101285844416803438769144460036425678359908733332182367776587521824356306545308780262109501, 16960598957857158004200152340723768697140517883876375860074812414430009210110270596775612236591317858945274366804448872120458103728483749408926203642159476]]) 

In [105]: B = Matrix([ 
    .....: [13229751631544149067279482127723938103350882358472000559554652108051830355519740001369711685002280481793927699976894894174915494730969365689796995942384549941729746359], 
    .....: [ 6297029075285965452192058994038386805630769499325569222070251145048742874865001443012028449109256920653330699534131011838924934761256065730590598587324702855905568792], 
    .....: [ 2716399059127712137195258185543449422406957647413815998750448343033195453621025416723402896107144066438026581899370740803775830118300462801794954824323092548703851334]]) 

In [106]: A.solve(B) 
Out[106]: 
Matrix([ 
[358183301733], 
[498758543457], 
[ 1919512167]]) 

In [107]: