Question

The program construes equations describing the concentrations of reagents in time in consecutive reactions A->B->C->...

http://en.wikipedia.org/wiki/Rate_equation#Consecutive_reactions

The algorithm has been devised by myself and allows one to completely ignore the need to solve the differential equations to obtain the concentration functions. It makes use of the pattern that emerges when one closely examines these equations.

The code that works uses floating point numbers, but they come at a disadvantage. When the rate constants k_n for any two given n values are similar in value, the inexactness of floating point numbers kicks in and the results drown in the ocean of error amplification:

graph

I figured that substituting float() with decimal() should fix this issue, but to my dismay, the changes caused unexpected errors:

After substituting float() with decimal():

Traceback (most recent call last):
  File "conreact9.py", line 280, in <module>
    exec(comm)
  File "<string>", line 1, in <module>
  File "conreact9.py", line 256, in graphit
    p += plt.plot(t,eval(c_(i)),label= "c_" + str(i) + "(t)")
  File "<string>", line 1, in <module>
  File "C:\Python27\lib\decimal.py", line 658, in __new__
    raise TypeError("Cannot convert %r to Decimal" % value)
TypeError: Cannot convert array([0, -0.1, -0.2, -0.3, -0.4, -0.5, -0.6, -0.7, -0
.8, -0.9, -1.0, -1.1,
       -1.2, -1.3, -1.4, -1.5, -1.6, -1.7, -1.8, -1.9, -2.0, -2.1, -2.2,
       -2.3, -2.4, -2.5, -2.6, -2.7, -2.8, -2.9, -3.0, -3.1, -3.2, -3.3,
       -3.4, -3.5, -3.6, -3.7, -3.8, -3.9, -4.0, -4.1, -4.2, -4.3, -4.4,
       -4.5, -4.6, -4.7, -4.8, -4.9], dtype=object) to Decimal

And this error when I substitute NumPy exp() to decimal equivalent:

Traceback (most recent call last):
  File "conreact10.py", line 280, in <module>
    exec(comm)
  File "<string>", line 1, in <module>
  File "conreact10.py", line 256, in graphit
    p += plt.plot(t,eval(c_(i)),label= "c_" + str(i) + "(t)")
  File "<string>", line 1, in <module>
AttributeError: exp

Didn't find anything using Google or official documentation on the decimal module.

Here's the untouched code for reference:

http://pastebin.com/mfuwVq7B

And here's my attempt at modifying it:

http://pastebin.com/xde4gvHu

This may as well be something very basic I'm not able to see. I'm a student at the Faculty of Chemistry and I'm a novice when it comes to programming (no previous background). Thank you very much.

Was it helpful?

Solution

The problem is the substitution of numpy.exp() with Decimal.decimal.exp().

The former takes an array and returns an array in which e has been raised to each element's power:

>>> x = np.linspace(-2*np.pi, 2*np.pi, 100)
>>> xx = x + 1j * x[:, np.newaxis] # a + ib over complex plane
>>> out = np.exp(xx)

The latter works on single numbers:

>>> decimal.Decimal(1).exp()
Decimal('2.718281828459045235360287471')

(examples taken from the pages linked above).

Licensed under: CC-BY-SA with attribution
Not affiliated with StackOverflow
scroll top