Euler-Cromer ODE Python例程

Question

我使用歐拉克羅默方案來計算哈雷彗星的位置和速度。該腳本爲一個範圍內的每個初始速度值嘗試幾個值作爲時間步長（tau）。對於每個tau值，它通過歐拉克羅默慣例，並比較從軌道開始到第一個軌道週期結束的總機械能。如果兩個能量之間的百分比差異小於1％，則將tau的當前（最優）值添加到列表中。經過所有迭代後，tau和初始速度的值用半對數圖上的pyplot繪圖，以便可以解釋哈雷彗星的真實初始速度。然而，我最優的taus列表中的每個元素都是我的tau範圍的第一個元素（當前爲0.1）。腳本似乎比必要的更復雜，也許有點複雜。那就是： 進口matplotlib.pyplot如PLT 進口numpy的爲NP

## Set physical parameters and initial position and velocity of the comet 
GM = 4 * np.pi ** 2    # Grav const. times mass of sun (AU^3/yr^2) 
mass = 1.      # Mass of comet 
r0 = 35       # Initial aphelion position 
v0 = np.arange(0.5, 1.1, 0.1) # Create a range of initial velocities in AU/year 


## Function that takes position and velocity vectors and the initial total energy 
## and outputs the value of tau that provides a less than 1% error in the total energy 
def tau_test(vel): 
    tau = np.arange(.1, .009, -.001) 
    r = [r0, 0] 
    v = vel 
    optimal_tau = 0 
    for t in tau: 
     i = 0 
     i_max = 5 * 76/t 
     r_mag_initial = np.sqrt(r[0] ** 2 + r[1] ** 2)  # initial magnitude of the position vector 
     v_mag_initial = np.sqrt(v[0] ** 2 + v[1] ** 2)  # initial magnitude of the velocity vector 
     KE_initial = 0.5 * mass * v_mag_initial ** 2  # initial kinetic energy 
     PE_initial = -GM * mass/r_mag_initial    # initial potential energy 
     ME_initial = KE_initial + PE_initial    # initial total energy 
     ## Main looping function (using Euler-Cromer method for now) 
     while not i == i_max: 
      r_mag = np.sqrt(r[0] ** 2 + r[1] ** 2)  # magnitude of the position vector 
      v_mag_a = np.sqrt(v[0] ** 2 + v[1] ** 2) # current velocity magnitude 
      a = [-GM * r[0]/r_mag ** 3,    # acceleration vector 
       -GM * r[1]/r_mag ** 3]    # based on the current position 
      v = [v[0] + t * a[0], v[1] + t * a[1]]  # new velocity vector using Euler-Cromer method 
      r = [r[0] + t * v[0], r[1] + t * v[1]]  # new position vector using Euler-Cromer method 
      v_mag_b = np.sqrt(v[0] ** 2 + v[1] ** 2) # new velocity magnitude 
      if i > .75 * 76/t:      # Ensure that the comet is at least past the sun 
       if v_mag_b > v_mag_a:     # Break out of the while loop if the comet is starting 
        break        # to move back toward the sun 
      i += 1 
     v_mag = np.sqrt(v[0] ** 2 + v[1] ** 2)   # final magnitude of the velocity vector 
     r_mag = np.sqrt(r[0] ** 2 + r[1] ** 2)   # final magnitude of the position vector 
     KE = 0.5 * mass * v_mag ** 2     # final kinetic energy 
     PE = -GM * mass/r_mag       # final potential energy 
     ME = KE + PE         # final total energy 
     #print ME 
     if abs(((ME - ME_initial)/ME_initial) * 100) <= 1: # If the percent error between the initial and final 
      optimal_tau = t          # total energies is less than 1%, set t as the optimal 
      break            # tau and break out of the for loop 
    return optimal_tau 

## Loop through each initial velocity and test it against several values of tau 
taus = [] 
for u in v0: 
    v = [0, u] 
    #print ME_initial 
    taus.append(tau_test(v)) 

#print taus 
## Plot the values of tau and initial velocity on a semi-log graph 
## so that the true velocity of Halley's comet can be interpreted 
plt.semilogy(v0, taus) 
plt.grid(True) 
plt.title('Tau versus Initial Velocities') 

plt.show() 

任何幫助表示讚賞。

Answer 1

我決定回到我的原始腳本上，描繪彗星的軌跡及其能量隨着時間的流逝。我編輯它以找出每個軌道之後初始總能量與彗星總能量之間的百分比差異。我繪製了總能量百分比差異趨勢，並簡單地重新計算程序幾次，對0.3到1之間的每個初始速度使用不同的時間步長值。雖然這基本上是一種強力解決方案，但它起作用。