i am doing some scientific calculations using python.I got into trouble when solving linear equations,where some coefficients are very large ~ E13 , some are very small ~E-69.It gives me a inaccurate solutions.
The equations,as folowing in the code, are very simple,just like rate[0]=rate[1]=...=rate[6] rate[7]=0.I use slove to get the solutions,but rate[0]!=rate[1]=rate[2]=...Although the most are right,but the physical meaning is totally wrong,which is unacceptable.rate[0]~rate[6] must be equal.
I tried some ways to promote the accuracy.
#convert float to symbols
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_]
But they both don't work.
I tried the same code on matlab,using solve or vpasolveinsymbolic math tool box.The solutions are right.But for some reasons,the calutions must be done using python.
So my question is how to promote the accuracy?
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)))
The result for 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
The most serious is rate[0] is negative and rate[6],which is supposed to be zero, has the largest absorbate value.
And the right solutions in matlab
6.245878938898438e-28
6.245878938898395e-28
6.245878938898395e-28
6.245878938898395e-28
6.245878938898395e-28
6.245878938898395e-28
6.245878938898395e-28
0
To avoid any loss of precision in floating point arithmetics, recast the given coefficients as rationals. Assuming from sympy import Rational, this is done with
kf_ = [Rational(x) for x in kf_]
kb_ = [Rational(x) for x in kb_]
Then use the code you have to solve the system and compute the rates. To display floating point numbers instead of giant rationals, use evalf method:
print('\n'.join(str(r.evalf()) for r in rate_))
prints
6.24587893889840e-28
6.24587893889840e-28
6.24587893889840e-28
6.24587893889840e-28
6.24587893889840e-28
6.24587893889840e-28
6.24587893889840e-28
0
Note: SymPy documentation recommends using linsolve for linear systems, but you'd need to adapt the code to deal with the different return type of linsolve.
Also, a linear system with numeric coefficients can be solved directly by mpmath which allows setting an arbitrarily large precision of floating point computations.