0

我正在编写一个程序，该程序利用欧拉方法来计算小行星的轨道是否会导致与地球的碰撞。在主循环的每次迭代结束时，都有一个 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')

这是代码：

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"

提前致谢！

4

1 回答 1

2

问题出在您的 Dist 函数中。当您计算两点之间的距离时，您将计算平方距离作为中间值。该中间值可能会溢出适中的大距离。 维基百科对这个问题及其解决方案进行了很好的讨论。简而言之，以下 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)

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

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

和 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 函数中出现域错误。如果 theta_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"
于 2013-11-12T08:41:11.483 回答