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 関数にあります。2 点間の距離を計算するときは、距離の 2 乗を中間値として計算します。その中間値は、適度に長い距離でオーバーフローする可能性があります。 ウィキペディアには、問題とその解決策に関する素晴らしい議論があります。つまり、 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)

これは、二乗距離を中間値として計算することを回避する数学的に同等の関数です。このオーバーフロー エラーを修正した後、同様の手法を使用して修正できるものをさらに 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 に答える