Regular Perturbation Solutions for Systems of 1st-Order ODE Initial-Value Problems

4100 days ago by medlock

# This cell is a bunch of code to find regular perturbation solutions # for systems of 1st-order ODE initial-value problems def regularPerturbationSeries(X, t, epsilon, order = 2): """Return regular perturbation series for a list of unknown functions.""" # Define constituent functions Xp = [[function('X_%d_%d' % (i, o), t) for o in range(order)] for i in range(len(X))] # Return perturbation series return [sum(epsilon^o * Xp[i][o] for o in range(order)) for i in range(len(X))] def reduceEquation(eq, epsilon, order, sol): """Pull out the coefficient of epsilon^order substituting in any info in sol.""" term = eq.diff(epsilon, order).subs(epsilon = 0) / factorial(order) for s in sol: for z in s: term = term.subs(z) return term def reduceEquations(eqs, epsilon, order, sol): """Pull out the coefficient of epsilon^order substituting in any info in sol.""" return [reduceEquation(eq, epsilon, order, sol) for eq in eqs] def solveReducedProblem(ODE, IC, Xp, t, epsilon, order, sol = []): """Solve the problem at the given order, given the solutions at lower orders in sol.""" reducedX = reduceEquations(Xp, epsilon, order, sol) reducedODE = reduceEquations(ODE(Xp, t), epsilon, order, sol) reducedIC = reduceEquations(IC(Xp, t), epsilon, order, sol) convertedIC = [0] + [I.rhs() for I in reducedIC] return desolve_system(reducedODE, reducedX, ics = convertedIC, ivar = t) def perturbationSolve(ODE, IC, X, t, epsilon, order = 2): """Return the perturbation solution of the system of 1st-order ODE initial-value problem.""" Xp = regularPerturbationSeries(X, t, epsilon) sol = [] for o in range(order): sol.append(solveReducedProblem(ODE, IC, Xp, t, epsilon, o, sol)) return [sum(epsilon^o * sol[o][i].rhs() for o in range(order)) for i in range(len(X))] 
       
# Set up the particular problem # Independent variable and parameter var('t epsilon') # Unknown functions u = function('u', t) q = function('q', t) # The system of 1st-order ODEs def ODE(X, t): (u, q) = X return [diff(u, t) == 1 - u * exp(epsilon * (q - 1)), diff(q, t) == u * exp(epsilon * (q - 1)) - q] # The initial conditions. They are all at t = 0. def IC(X, t): (u, q) = X return [u(t = 0) == 0, q(t = 0) == 0] 
       
# 1-term solution P1 = perturbationSolve(ODE, IC, (u, q), t, epsilon, 1) show(P1) # 2-term solution P2 = perturbationSolve(ODE, IC, (u, q), t, epsilon, 2) show(P2) 
       
\newcommand{\Bold}[1]{\mathbf{#1}}\left[-e^{\left(-t\right)} + 1, -t e^{\left(-t\right)} - e^{\left(-t\right)} + 1\right]
\newcommand{\Bold}[1]{\mathbf{#1}}\left[\frac{1}{2} \, {\left(t^{2} e^{\left(-t\right)} + 2 \, t e^{\left(-2 \, t\right)} + 2 \, t e^{\left(-t\right)} + 4 \, e^{\left(-2 \, t\right)} - 4 \, e^{\left(-t\right)}\right)} \epsilon - e^{\left(-t\right)} + 1, \frac{1}{6} \, {\left(t^{3} e^{\left(-t\right)} - 12 \, t e^{\left(-2 \, t\right)} - 18 \, t e^{\left(-t\right)} - 30 \, e^{\left(-2 \, t\right)} + 30 \, e^{\left(-t\right)}\right)} \epsilon - t e^{\left(-t\right)} - e^{\left(-t\right)} + 1\right]
\newcommand{\Bold}[1]{\mathbf{#1}}\left[-e^{\left(-t\right)} + 1, -t e^{\left(-t\right)} - e^{\left(-t\right)} + 1\right]
\newcommand{\Bold}[1]{\mathbf{#1}}\left[\frac{1}{2} \, {\left(t^{2} e^{\left(-t\right)} + 2 \, t e^{\left(-2 \, t\right)} + 2 \, t e^{\left(-t\right)} + 4 \, e^{\left(-2 \, t\right)} - 4 \, e^{\left(-t\right)}\right)} \epsilon - e^{\left(-t\right)} + 1, \frac{1}{6} \, {\left(t^{3} e^{\left(-t\right)} - 12 \, t e^{\left(-2 \, t\right)} - 18 \, t e^{\left(-t\right)} - 30 \, e^{\left(-2 \, t\right)} + 30 \, e^{\left(-t\right)}\right)} \epsilon - t e^{\left(-t\right)} - e^{\left(-t\right)} + 1\right]
# Set up numerical solution # ODE right-hand side def RHS(t, Y, p): (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 T = ode_solver(function = RHS, y_0 = Y0, algorithm = 'rk8pd') 
       
epsilon = 0.3 tMax = 10 # Pull apart perturbation solutions (uP1, qP1) = P1 (uP2, qP2) = P2 # Numerically solve T.ode_solve(params = [epsilon], t_span = [0, tMax], num_points = 1000) # Store solution in handy format uN = [(s[0], s[1][0]) for s in T.solution] qN = [(s[0], s[1][1]) for s in T.solution] # Plot u (plot(uP1(epsilon = epsilon), (t, 0, tMax), color = 'red') + plot(uP2(epsilon = epsilon), (t, 0, tMax), color = 'green') + line(uN, color = 'blue')).show(axes_labels = ('t', 'u')) # Plot q (plot(qP1(epsilon = epsilon), (t, 0, tMax), color = 'red') + plot(qP2(epsilon = epsilon), (t, 0, tMax), color = 'green') + line(qN, color = 'blue')).show(axes_labels = ('t', 'q'))