Regular Perturbation for System of 2 1st-Order ODEs

3720 days ago by medlock

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