当解方程的系数在E13到E18之间的线性系统时，提高精度

Question

我在用python做一些科学计算。我在解线性方程组时遇到了麻烦，有些系数非常大~E13，有些系数非常小~E-69。它给了我一个不准确的解。你知道吗

下面代码中的方程非常简单，就像rate[0]=rate[1]=...=rate[6] rate[7]=0，我用slove来求解，但是rate[0]!=rate[1]=rate[2]=...虽然大部分是对的，但是物理意义是完全错误的，这是错误的不合格率[0]~速率[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上尝试了相同的代码，在symbolic math tool box中使用了solve或vpasolve是的，但是由于某些原因，必须使用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)))

速率[0]~速率[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

最严重的是比率[0]为负，比率[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返回类型。你知道吗

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