title: Euler-Cromer ODE Python Routine
body: <p>I am using an Euler-Cromer scheme to calculate the position and velocity of Halley's comet. The script tries several values for a time-step (tau) for each value of initial velocity in a range. For each tau value, it runs through the Euler-Cromer routine and compares the total mechanical energy from the beginning of the orbit to the end of the first orbit's cycle. If the percent difference between the two energies is less than 1%, the current (optimal) value of tau is added to a list. After all the iterations, the values of tau and initial velocities are plotted using pyplot on a semi-log graph so that the true initial velocity of Halley's comet can be interpreted. However, each element of my optimal taus list is the first element of my tau range (currently 0.1). The script seems more complex than necessary and perhaps a bit convoluted. Here it is:
      import matplotlib.pyplot as plt
      import numpy as np</p>
  
  <pre><code>## 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 &gt; .75 * 76 / t:                        # Ensure that the comet is at least past the sun
                  if v_mag_b &gt; 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) &lt;= 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()
  </code></pre>
  
  <p>Any help is appreciated.</p>