var('t')
t t |
# Perturbation solutions
u0(t) = 1 - exp(-t)
q0(t) = 1 - (1 + t) * exp(-t)
u1(t) = (2 + t) * exp(-2 * t) + (1/2 * t^2 + t - 2) * exp(-t)
q1(t) = -(5 + 2 * t) * exp(-2 * t) + (1/6 * t^3 - 3 * t + 5) * exp(-t)
|
|
# Numerical solution
def RHS(t, Y, p):
'ODE right-hand side.'
(u, q) = Y
epsilon = p[0]
du = 1 - u * exp(epsilon * (q - 1))
dq = u * exp(epsilon * (q - 1)) - q
return (du, dq)
# Initial condition
Y0 = (0, 0)
# Set up solver
solver = ode_solver(function = RHS, y_0 = Y0, algorithm = 'rk8pd')
|
|
epsilon = 0.3
tMax = 10
# Numerically solve
solver.ode_solve(params = [epsilon], t_span = [0, tMax], num_points = 1000)
# Plot u
(plot(u0, (t, 0, tMax), color = 'red')
+ plot(u0 + epsilon * u1, (t, 0, tMax), color = 'green')
+ line(((s[0], s[1][0]) for s in solver.solution), color = 'blue')).show(axes_labels = ('t', 'u'))
# Plot 1
(plot(q0, (t, 0, tMax), color = 'red')
+ plot(q0 + epsilon * q1, (t, 0, tMax), color = 'green')
+ line(((s[0], s[1][1]) for s in solver.solution), color = 'blue')).show(axes_labels = ('t', 'q'))
|
|
|
|