在python中限制浮点值以避免溢出

Question

我正在写一个程序，利用欧拉的方法来计算小行星的轨道是否会导致与地球的碰撞。在主循环的每次迭代结束时，都有一个if语句，它使用小行星和地球之间的距离来确定是否会发生碰撞。当我试着运行这个程序时，我得到一个溢出错误：数值结果超出范围，我假设这是由于我正在使用三角函数来转换极坐标和极坐标之外的事实，并且想知道如何限制这些函数返回的浮点值的大小以修复错误？在

编辑：以下是例外：

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')

代码如下：

^{pr2}$

提前谢谢！在

Answer 1

问题出在你的Dist函数上。计算两点之间的距离时，将平方距离计算为中间值。这个中间值可以在相当大的距离内溢出。Wikipedia对这个问题及其解决方案进行了很好的讨论。简而言之，以下dist功能的替换将解决您当前的问题：

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)

它只是一个数学上等价的函数，它避免了将平方距离作为中间值来计算。在修复了这个溢出错误之后，我又得到了两个可以使用类似技术修复的错误。我把你的重力函数改成这个：

^{pr2}$

以及v\r公式：

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

从您的原始版本：

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

但是，仍然存在一些问题。一旦我做了这些更改，我就可以避免溢出错误，但是当theta_a变得太大时，cos函数中会出现域错误。如果θa是我认为的，你可以通过添加一个mod 2*pi来解决最后一个问题，如下所示：

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

代替：

theta_a += r_a*theta_a*delta_t

下面是所有更改后的工作代码。我不确定这是否正确，但没有出现错误。在

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"