Search code examples
pythonnumpyprecisionnumericinfinity

Handling operations with infinities in python


I have a piece of code that does a simple calculation.

import numpy as np

#Constants
R = 8.314462
T = 298.15
e = -678.692
e_overkbT = e*1000/(R*T)

#Independent variable
mu = np.linspace(-2000,2000,1000)  
mu_overkbT = mu*1000/(R*T)


#Calculation
aa = (np.exp(mu_overkbT- e_overkbT))
theta = aa/(1+aa)

For negative values of 'mu', 'aa' is very small and thus the variable "theta" is very close to 0. For positive values of 'mu', 'aa' is very large. Thus for large numbers 'theta' approaches 1. (large number over large number + 1).

For large values of 'aa' python rounds 'theta' to be 1, which is fine. However, eventually for large enough numbers python will say 'aa' is 'inf'. Thus in the final step of calculating 'theta' I encounter a runtime error of dividing 'inf'/'inf'.

I need someway to handle this error such that it gives me '1' as the result for 'theta'. I can't reduce the range of the variable 'mu' and stop before the error, because this calculation is inside of a large function that changes the value of 'e', and thus this error does not always occur at the same spot.

Thanks.


Solution

  • Such overflow happens very often when using the exponential function on large terms. Other than the good very good comment noting that exp(x)/(1+exp(x)) = 1/(1+exp(-x)), another general approach in case you don't find easy transformations is to use the logarithm to make intermediary numbers more manageable, and then in the end to reverse this operation. This is especially true with products of many large (or very small) terms, which by applying the logarithm become a simple sum.