Non-constant Gravity

3724 days ago by medlock

var('t') 
       
t
t
# Perturbation solutions x0(t) = 1/2 * t * (2 - t) x1(t) = 1/12 * t^3 * (4 - t) 
       
# Numerical solution def RHS(t, Y, p): 'ODE right-hand side.' (x, xp) = Y epsilon = p[0] dx = xp dxp = - 1 / (1 + epsilon * x)^2 return (dx, dxp) # Initial condition Y0 = (0, 1) # Set up solver solver = ode_solver(function = RHS, y_0 = Y0, algorithm = 'rk8pd') 
       
epsilon = 0.1 # Initial velocity = 18000 * sqrt(epsilon) mi/hr! tMax0 = 2 tMax1 = 2.14 # Numerically solve solver.ode_solve(params = [epsilon], t_span = [0, tMax1], num_points = 1000) # Plot (plot(x0, (t, 0, tMax0), color = 'red') + plot(x0 + epsilon * x1, (t, 0, tMax1), color = 'green') + line(((s[0], s[1][0]) for s in solver.solution), color = 'blue')).show(axes_labels = ('t', 'x')) 
       
error0 = [(s[0], x0(t = s[0]) - s[1][0]) for s in solver.solution] error1 = [(s[0], (x0 + epsilon * x1)(t = s[0]) - s[1][0]) for s in solver.solution] (line(error0, color = 'red') + line(error1, color = 'green')).show(axes_labels = ('t', 'error')) line(error1, color = 'green').show(axes_labels = ('t', 'error'))