# 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))]