在求解方程系数范围为 E13 到 E-18 的线性系统时提高精度

Question

我正在使用 python.I 进行一些科学计算，在求解线性方程时遇到了麻烦，其中一些系数非常大 ~ E13，一些非常小 ~E-69.It 给我一个不准确的结果解决方案。

方程式，如下代码，很简单，就像rate[0]=rate[1]=...=rate[6] rate[7]=0。我用slove求解，但是rate[0]!=rate[1]=rate[2]=...虽然大部分都是对的，但物理意义完全错误，即unacceptable.rate[0]~rate[6]必须相等

我尝试了一些方法来提高准确性。

#convert float to symbols
kf_ = [symbols(str(k)) for k in kf_]
kb_= [symbols(str(k)) for k in kb_]

或

#convert float to decimal
kf_ = [Decimal(str(k)) for k in kf_]
kb_ = [Decimal(str(k)) for k in kb_]

但它们都不起作用。我在 matlab 上尝试了相同的代码，使用 solve 或 vpasolve 在 symbolic math tool box 上。由于某些原因，解决方案是 right.But，必须使用 python。所以我的问题是如何提高准确性？

from sympy import symbols , solve
from decimal import Decimal

#coeffcients   Some are very large ~ E13 , some are very small ~E-69
kf_= [804970533.1611289,
1.5474657692374676e-13,
64055206.72863516,
43027484.879109934,
239.58564380236825,
43027484.879109934,
0.6887097015872349,
43027484.879109934]

kb_=[51156022807606.22,
4.46863981338889e-18,
9.17599257631182,
8.862701377525092e-43,
4.203415587017237e-20,
2180151.4516747626,
5.590961781720337e-69,
0.011036598954049947]

#convert float to symbols , it takes quite a long time
#kf_ = [symbols(str(k)) for k in kf_]
#kb_= [symbols(str(k)) for k in kb_]
#or
#convert float to decimal
#kf_ = [Decimal(str(k)) for k in kf_]
#kb_ = [Decimal(str(k)) for k in kb_]

# define unkown
theta = list(symbols('theta1:%s' % (8 + 1)))

#define some expressions
rate=[0]*8
rate[0] = kf_[0] * theta[0] - kb_[0] * theta[1]
rate[1] = kf_[1] * theta[1] - kb_[1] * theta[2]
rate[2] = kf_[2] * theta[2] - kb_[2] * theta[3]
rate[3] = kf_[3] * theta[3] - kb_[3] * theta[4]
rate[4] = kf_[4] * theta[4] - kb_[4] * theta[5]
rate[5] = kf_[5] * theta[5] - kb_[5] * theta[6]
rate[6] = kf_[6] * theta[6] - kb_[6] * theta[0]
rate[7] = kf_[7] * theta[0] - kb_[7] * theta[7]

print('\n'.join(str(r) for r in rate))

#euqations
fun=[0]*8
fun[0]=sum(theta)-1
fun[1]=rate[0]-rate[1]# The coefficients kb[0] (~E13 )and kf[1] (~E-13) are merged
fun[2] = rate[1] - rate[2]
fun[3] = rate[2] - rate[3]
fun[4] = rate[3] - rate[4]
fun[5] = rate[4] - rate[5]
fun[6] = rate[5] - rate[6]
fun[7] = rate[7]

#solve
solThetaT=solve(fun,theta)
#print(solThetaT)
theta_=[solThetaT[t]  for t in theta]
#print(theta_)
rate_=[0]*8
for i in range(len(rate)):
    rate_[i]=rate[i].subs(solThetaT)

print('\n'.join(str(r) for r in rate_))

#when convert float to symbols
#for r in rate_:
    #print(eval(str(r)))

rate[0]~rate[7]的结果:

-1.11022302462516e-16
6.24587893889839e-28
6.24587893889839e-28
6.24587893889840e-28
6.24587893889840e-28
6.24587893222751e-28
6.24587895329296e-28
-3.81639164714898e-17

最严重的是rate[0]为负，rate[6]，本应为零，吸收率最大。

以及 matlab 中的正确解决方案

 6.245878938898438e-28
 6.245878938898395e-28
 6.245878938898395e-28
 6.245878938898395e-28
 6.245878938898395e-28
 6.245878938898395e-28
 6.245878938898395e-28
                 0

Answer 1

为避免浮点运算中的任何精度损失，请将给定系数重新转换为有理数。假设 from sympy import Rational，这是用

完成的

kf_ = [Rational(x) for x in kf_]
kb_ = [Rational(x) for x in kb_]

然后使用您必须的代码来求解系统并计算比率。要显示浮点数而不是巨型有理数，请使用 evalf 方法：

print('\n'.join(str(r.evalf()) for r in rate_))

打印

6.24587893889840e-28
6.24587893889840e-28
6.24587893889840e-28
6.24587893889840e-28
6.24587893889840e-28
6.24587893889840e-28
6.24587893889840e-28
0

注意：SymPy 文档建议对线性系统使用 linsolve，但您需要调整代码以处理 linsolve 的不同 return 类型。

此外，具有数值系数的线性系统可以直接通过 mpmath 求解，这允许设置任意大的浮点计算精度。

在求解方程系数范围为 E13 到 E-18 的线性系统时提高精度

Improve accuracy when solving a linear system where equations have coefficients ranging from E13 to E-18

python

sympy

floating-accuracy

equation-solving