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