2

我在我的 Python 程序中使用 numpy 和 mpmath。我使用 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

有没有办法避免这些错误？

4

2 回答 2

2

numpyndarrays 具有同质类型。当你 make 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)

但请注意，这样做会严重影响性能。

于 2014-01-16T15:30:08.887 回答
1

为了完整性，以及像我这样偶然发现这个问题的人，因为 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]:
于 2014-08-12T11:44:31.110 回答