Divison by zero due to insufficient precision using exp() Python

Question

I'm using a bissection algorithm to find the maximum of a function.

Here is my code:

from math import exp

h = 6.62606876e-34
c = 2.99792458e8
k = 1.3806504e-23
T = 7000
epsilon = 1.e-10

def fct(lam) :
    return 2*h**2*c**3*exp(h*c/(lam*k*T))/(k*T*lam**7*(exp(h*c/(lam*k*T)) - 1)) - 10*c**2*h/(lam**6*(exp(h*c/(lam*k*T)) - 1))

lam1 = 100.e-9
lam2 = 1000.e-9

delta = lam2 - lam1

f1 = fct(lam1)
f2 = fct(lam2)

while delta > epsilon :
    lamm = 0.5*(lam1 + lam2)
    fm = fct(lamm)
    f2=fct(lam2)
    if fm*f2 > 0 :
        lam2 = lamm
    else :
        lam1 = lamm
    delta = lam2 - lam1
    print('racine de la fonction: {}+/-{}'.format(lamm,delta))

The problem is that because of the finite precision of floats, I get a division by zero error when the line fm = fct(lamm) is evaluated.

How can I fix this issue?

Huh I can't reproduce you're error, everything appears to be working... I'm running 2.7 though so maybe that has something to do with it. — John, Jan 23 '16 at 20:55

score 1 · Answer 1 · edited May 23 '17 at 10:27

1

I think your best bet for doing math in python is using one of the many math processing libraries. https://stackoverflow.com/a/13442742/2534876 suggests mpmath for your floating point issues. bigfloat could possibly be sufficient. numpy also has stuff like float128 for this purpose.

Dig around and find something you like.

edited May 23 '17 at 10:27

Community

1
1

answered Jan 23 '16 at 20:52

JDong

2,304
3
24
42

score 1 · Answer 2 · answered Jan 23 '16 at 21:02

I don't have a great answer. One technique would be to recale your units, say set c=1 if you're dealing with large quantities, do your computation then express your answer back in the unscaled system. You'd have to do something analogous if you're dealing with small units.

Another solution which I've just discovered is the mpmath module, which claims to be able to do computations up to 1000 units of precision! Their website is http://mpmath.org, I hope this works out for you.

score 0 · Answer 3 · answered Jan 23 '16 at 23:10

One possibility is to transform your code into using exact arithmetics for all but the transcendental functions (exp in your example), for example using https://docs.python.org/2/library/fractions.html. Certainly, inputs and outputs of transcendental functions need to be converted.

Divison by zero due to insufficient precision using exp() Python

3 Answers3