Euler-Cromer ODE Python Routine
Question: 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

    
    
    ## 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()
    

Any help is appreciated.


Answer: Instead of banging my head against the wall, I decided to go back to my
original scrip that plots the trajectory of the comet and its energy over
time. I edited it to find the percent difference between the initial total
energy and the total energy of the comet after each orbit. I plotted the total
energy percent difference trend and simply reran the program several times,
using different time step values for each initial velocity between 0.3 and 1.
While this was essentially a brute-force solution, it worked.

