Restricting floating point values to avoid overflow in python

Question

I'm writing a program that utilizes Euler's Method to calculate whether or not an asteroid's orbit will result in a collision with Earth. At the end of each iteration of the main loop, there is an if statement that uses the distance between the asteroid and earth to determine whether or not a collision would have occurred. When I try running the program I get a Overflow error: numerical result is out of range, I assume that this is due to the fact that I'm using trigonometric functions to convert to and out of polar coordinates and was wondering how I would restrict the size of the floating point values returned by these functions so as to fix the error?

EDIT: Here's the exception:

Traceback (most recent call last):
  File "/home/austinlee/eclipse/plugins/org.python.pydev_2.7.0.2013032300/pysrc/pydevd.py", line 1397, in <module>
    debugger.run(setup['file'], None, None)
  File "/home/austinlee/eclipse/plugins/org.python.pydev_2.7.0.2013032300/pysrc/pydevd.py", line 1090, in run
    pydev_imports.execfile(file, globals, locals) #execute the script
  File "/home/austinlee/workspace/test/src/orbit.py", line 72, in <module>
    if( Dist(x_a, x_e, y_a, y_e) < d_close):
  File "/home/austinlee/workspace/test/src/orbit.py", line 37, in Dist
    return sqrt((b-a)**2+(d-c)**2)
OverflowError: (34, 'Numerical result out of range')

Here's the code:

from math import *

##position of earth and ast. relative to sun, units in m/s/kg

r_earth = 1.4959787E11
x_e = r_earth
y_e = 0
v_ye = 29784.3405
v_xe = 0

x_a = 1.37801793E11
y_a = 2.31478719E11
v_ya = -14263.6905
v_xa = -8490.32975

##constants- masses and radius's of objects
G = 6.67384E-11
M_sun = 1.988500E30
M_earth = 5.9726E24
R_earth = 6371.0E3
M_ast = 1.30E11
R_ast = 250
t = 0 
delta_t = 10
t_max = 10000000
a_xe = 0
a_ye = 0
a_xa = 0
a_ya = 0


##Define Acceleration due to Gravity and Distance Formulas
def Grav(a,b):
    return (G*a)/(b**2)

def Dist(a,b,c,d):
    return sqrt((b-a)**2+(d-c)**2)

##Derived Constants
t_close = 0
d_close = Dist(x_a,x_e,y_a,y_e)
r_a = Dist(0,x_a,0,y_a)
theta_e = 0
theta_a = atan2(y_a,x_a)
v_angle = sqrt(v_xa**2+v_ya**2)/r_a
v_r1 = v_angle
v_r = sqrt(v_xa**2+v_ya**2)
T = 2* pi/(365*24*3600)
a_r = v_xa**2+v_ya**2
a_theta = (-Grav(M_sun, Dist(x_a,0,y_a,0))-Grav(M_earth,Dist(x_a,x_e,y_a,y_e)))**2-a_r**2

## Main Loop- assuming constant, circular orbit for earth (i.e M_ast is negligible)

for t in range(0, t_max):
    t += delta_t
    theta_e = T*t
    x_e = r_earth*cos( theta_e )
    y_e = r_earth*sin( theta_e )

## Convert asteroid parameters into polar coordinates and solve using Euler's Method

    a_r = v_xa**2+v_ya**2
    a_theta = (-Grav(M_sun, Dist(x_a,0,y_a,0))-Grav(M_earth,Dist(x_a,x_e,y_a,y_e)))**2-a_r**2 
    v_r1 = v_r
    v_r += (a_r + r_a*v_angle**2)*delta_t
    v_angle += (a_theta - 2*v_r1*v_angle)* delta_t
    theta_a += r_a*theta_a*delta_t
    r_a += v_r*delta_t
    x_a = r_a*cos( theta_a)
    y_a = r_a*sin( theta_a)

## Test for minimum distance  

    if( Dist(x_a, x_e, y_a, y_e) < d_close):
        d_close = Dist( x_a, x_e, y_a, y_e)
        t_close = t/(3600*24)
    continue
##Print Results:

print "d_close: ", d_close/1000, "km"
print "t_close: ", t_close
if( d_close < R_earth):
    print "Impact: Y"

else:
    print "Impact: N"

Thanks in advance!

According to the [documentation](http://docs.python.org/2/library/math.html#constants) for the `math` module python doesn't do anything special about the results. If you don't want to get an `OverflowError` just do `try: result = #action-that-raises-error except OverflowError: result = float('+inf') # or a big enough value`. — Bakuriu, Nov 12 '13 at 08:02

score 2 · Accepted Answer · answered Nov 12 '13 at 08:41

The problem is in your Dist function. When you calculate the distance between two points, you calculate the squared distance as an intermediate. That intermediate value can overflow for moderately large distances. Wikipedia has a nice discussion of the problem and its solution. In short, the following replacement for the dist function will solve your immediate problem:

def Dist(a,b,c,d):
    x = float(b - a)
    y = float(d - c)
    u = min(x, y)
    v = max(x, y)
    return abs(v) * sqrt(1 + (u/v)**2)

It's just a mathematically equivalent function that avoids calculating the squared distance as an intermediate. After fixing this overflow error, I got two more that can be fixed using similar techniques. I changed your Grav function to this:

def Grav(a,b):
    return (G*a/b)/(b)

and v_r formula to:

v_r += (a_r/v_angle + r_a*v_angle)*delta_t*v_angle

from your original:

v_r += (a_r + r_a*v_angle**2)*delta_t

However, there are still problems. Once I make those changes I am able to avoid the overflow errors, but eventually get a domain error in the cos function when theta_a gets too large. If theta_a is what I think it is, you can fix this last problem by adding a mod 2*pi, like this:

theta_a += r_a*theta_a*delta_t % (2*pi)

in place of:

theta_a += r_a*theta_a*delta_t

Below is the working code after all changes. I'm not sure if it's right, but there are no errors raised.

from math import *

##position of earth and ast. relative to sun, units in m/s/kg

r_earth = 1.4959787E11
x_e = r_earth
y_e = 0
v_ye = 29784.3405
v_xe = 0

x_a = 1.37801793E11
y_a = 2.31478719E11
v_ya = -14263.6905
v_xa = -8490.32975

##constants- masses and radius's of objects
G = 6.67384E-11
M_sun = 1.988500E30
M_earth = 5.9726E24
R_earth = 6371.0E3
M_ast = 1.30E11
R_ast = 250
t = 0 
delta_t = 10
t_max = 10000000
a_xe = 0
a_ye = 0
a_xa = 0
a_ya = 0


##Define Acceleration due to Gravity and Distance Formulas
def Grav(a,b):
    return (G*a/b)/(b) #Changed by jcrudy

def Dist(a,b,c,d): #Changed by jcrudy
    x = float(b - a)
    y = float(d - c)
    u = min(x, y)
    v = max(x, y)
    return abs(v) * sqrt(1 + (u/v)**2)

#    return sqrt((b-a)**2+(d-c)**2)

##Derived Constants
t_close = 0
d_close = Dist(x_a,x_e,y_a,y_e)
r_a = Dist(0,x_a,0,y_a)
theta_e = 0
theta_a = atan2(y_a,x_a)
v_angle = sqrt(v_xa**2+v_ya**2)/r_a
v_r1 = v_angle
v_r = sqrt(v_xa**2+v_ya**2)
T = 2* pi/(365*24*3600)
a_r = v_xa**2+v_ya**2
a_theta = (-Grav(M_sun, Dist(x_a,0,y_a,0))-Grav(M_earth,Dist(x_a,x_e,y_a,y_e)))**2-a_r**2

## Main Loop- assuming constant, circular orbit for earth (i.e M_ast is negligible)

for t in range(0, t_max):
    t += delta_t
    theta_e = T*t
    x_e = r_earth*cos( theta_e )
    y_e = r_earth*sin( theta_e )

## Convert asteroid parameters into polar coordinates and solve using Euler's Method

    a_r = v_xa**2+v_ya**2
    a_theta = (-Grav(M_sun, Dist(x_a,0,y_a,0))-Grav(M_earth,Dist(x_a,x_e,y_a,y_e)))**2-a_r**2 
    v_r1 = v_r
    v_r += (a_r/v_angle + r_a*v_angle)*delta_t*v_angle # Changed by jcrudy
    v_angle += (a_theta - 2*v_r1*v_angle)* delta_t
    theta_a += r_a*theta_a*delta_t % (2*pi) # Changed by jcrudy
    r_a += v_r*delta_t
    x_a = r_a*cos( theta_a)
    y_a = r_a*sin( theta_a)

## Test for minimum distance  

    if( Dist(x_a, x_e, y_a, y_e) < d_close):
        d_close = Dist( x_a, x_e, y_a, y_e)
        t_close = t/(3600*24)
    continue
##Print Results:

print "d_close: ", d_close/1000, "km"
print "t_close: ", t_close
if( d_close < R_earth):
    print "Impact: Y"

else:
    print "Impact: N"

Restricting floating point values to avoid overflow in python

1 Answers1