Heun Project

import numpy as np import matplotlib.pyplot as plt

p = (4*9.806) * 1 v = 0 xi = 0 #initial position xm = 2.12/1.057e-16 #Half the total distance xf = 4.24/1.057e-16 #Total Distance F = 4 *9.806 #force c= 2.99e8 m = 1 gamma = 1 / ((1-(v/c)**2)**0.5) ke = m*c**2*(gamma - 1) dt = 604800.0 #1 week in seconds t = 0 tau = 0.0 s = np.array([xi,p,tau]) def derivs_accel(s,t): #this is setting up the derivatives to use in HeunStep x = s[0] p = s[1] tau = s[2] c= 2.99e8 m = 1 v = (p * c) / (p**2 + c**2)**0.5 gamma = 1 / ((1-(v/c)**2)**0.5) if s[0] < xm: dpdt = (4*9.806) else: dpdt = -4*9.806 dxdt = v dtaudt = 1 / gamma return np.array([dxdt,dpdt,dtaudt]) def HeunStep(s, t, derivs, dt): #explained in Algorithm section f1=derivs(s,t) f2=derivs(s+f1*dt,t+dt) return s + 0.5*(f1+f2)*dt tlist = [t] #Lists created to store the data for each point xlist = [xi] plist = [p] taulist = [tau] glist = [gamma] kelist = [ke] vlist = [v] while s[1] > 0: s = HeunStep(s, t, derivs_accel, dt) #calling the Heun Step v = (s[1] * c) / (s[1]**2 + c**2)**0.5 #calculatin velocity gamma = 1 / ((1-(v/c)**2)**0.5) #calculating gamma ke = c**2 * (gamma - 1) #calculating kinetic energy t+= dt xlist.append(s[0]) #appending and adding the data to the lists plist.append(s[1]) taulist.append(s[2]) vlist.append(v) glist.append(gamma) kelist.append(ke) tlist.append(t) tlist = [i / 31556952 for i in tlist] #converting time to years taulist = [i / 31556952 for i in taulist]

#All of the plots are set up like this plt.xlabel("Time(years)") plt.ylabel("Position(m)") plt.title("Motion of the Spaceship") plt.plot(tlist, xlist,'r-') plt.grid()

print("The journey took", tlist[-1],"years") print("The dialated time was", taulist[-1], 'years') print("Max Kinetic Energy:", max(kelist), "J") print("Work at halfway point =", F*xm, "J") uncertainty = abs(F*xm - max(kelist)) print("Uncertainty is: +/-", format(uncertainty, ".3E"), "J") percenterror = (max(kelist) - (F*xm)) / (F*xm) * 100 print("Percent error is", abs(percenterror))