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'))
|
|
|
|