# Non-constant Gravity

## 3901 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 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, s) for s in solver.solution), color = 'blue')).show(axes_labels = ('t', 'x'))  error0 = [(s, x0(t = s) - s) for s in solver.solution] error1 = [(s, (x0 + epsilon * x1)(t = s) - s) 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'))    