def ODERHS(t, Y, params):
(N, P) = Y
(alpha, beta, gamma) = params
return (dNdt(N = N, P = P, alpha = alpha, beta = beta, gamma = gamma),
dPdt(N = N, P = P, alpha = alpha, beta = beta, gamma = gamma))
def computeSolution(alpha, beta, gamma, IC):
T = ode_solver(function = ODERHS, algorithm = 'rk8pd')
T.ode_solve(y_0 = IC, t_span=[0, 100], params = (alpha, beta, gamma),
num_points = 200)
(t, Y) = zip(*T.solution)
(N, P) = zip(*Y)
return (t, N, P)
def numericalSolution(alpha, beta, gamma, IC):
(t, Nn, Pn) = computeSolution(alpha, beta, gamma, IC)
return line(zip(t, Nn), color = 'red') + line(zip(t, Pn), color = 'green')
def numericalPhaseTrajectory(alpha, beta, gamma, IC):
(t, Nn, Pn) = computeSolution(alpha, beta, gamma, IC)
return line(zip(Nn, Pn))
ICs = [(n, 0.99) for n in (0.01, 0.1, 0.25, 0.5, 0.75, 0.9, 0.99)] + [(0.99, p) for p in (0.01, 0.1, 0.25, 0.5, 0.75, 0.9, 0.99)]
def numericalPhaseTrajectories(alpha, beta, gamma):
return phasePortrait(alpha, beta, gamma) \
+ sum(numericalPhaseTrajectory(alpha, beta, gamma,
(1.1 * gamma * N0,
1.2 * (gamma + 1)^2 / 4 / gamma * P0))
for (N0, P0) in ICs)