Code
using NonlinearSolve
using DifferentialEquations
using SciMLBase
using StaticArrays
using ForwardDiff
using LinearAlgebra
using PlotlyJSThe following Julia packages will be needed for this document:
In the document below, many code cells have an extra
;after the code to suppress junk from being output into this document.
Phase Portrait Analysis and Plotting Functions
Fixed Points
There might be multiple fixed points for a given phase portrait, and the following function finds fixed points within a given region with NonlinearSolve.
Code
function find_fixed_points(system_of_eqs; guess_xs::AbstractRange, guess_ys::AbstractRange)
fixed_points = []
for (guess_x, guess_y) ∈ Base.product(guess_xs, guess_ys)
u0 = SA[guess_x, guess_y]
prob = NonlinearProblem(system_of_eqs, u0)
sol = solve(prob, NewtonRaphson())
if SciMLBase.successful_retcode(sol)
found = false
for fixed_point ∈ fixed_points
if isapprox(fixed_point[1], sol.u[1], atol = 1e-3) &&
isapprox(fixed_point[2], sol.u[2], atol = 1e-3)
found = true
break
end
end
if !found
push!(fixed_points, sol.u)
end
end
end
return fixed_points
end;Find Jacobians at Fixed Points
The following function finds the Jacobians for the given system of equations and fixed points.
Code
function find_jacobians(system_of_eqs, fps)
jacobians = []
for fp ∈ fps
jacobian = ForwardDiff.jacobian(system_of_eqs, fp)
push!(jacobians, jacobian)
end
return jacobians
end;Classify Fixed Points
This function classifies a fixed point depending trace and determinant of the Jacobian at that fixed point.
Code
function classify_jacobian(A::Matrix{Float64})
τ = tr(A)
Δ = det(A)
discriminant = tr(A)^2 - 4*det(A)
if Δ < 0.0
return "Saddle"
else
if isapprox(discriminant, 0.0)
return "Star, Degenerate"
elseif discriminant > 0.0
if isapprox(τ, 0.0)
return "Neutral Stable"
elseif τ < 0.0
return "Stable"
else
return "Unstable"
end
else
return "Spiral"
end
end
return "Unknown"
end;Nullcline Contours
Code
function nullcline_contours(f, g, xs, ys)
f_xy = [f([x, y]) for x ∈ xs, y ∈ ys]
g_xy = [g([x, y]) for x ∈ xs, y ∈ ys]
return f_xy, g_xy
end;Slope Field Calculation
This function calculates the slope field defined by f and g over the points defined by xs and ys.
Code
function slope_field(f, g, xs, ys)
scaler = 1 / length(xs)
start_xys = Base.product(xs, ys)
end_xys = [
(start_xy[1] + f(start_xy)*scaler, start_xy[2] + g(start_xy)*scaler) for
start_xy ∈ start_xys
]
return start_xys, end_xys
end;Trajectory Calculation
This function calculates trajectories of the trajectory_eqs! from the initial conditions given in u0 over the spans given by tspans.
Code
function calculate_trajectories(trajectory_eqs!, u0s, tspans)
trajectories = []
for (u0, tspan) ∈ zip(u0s, tspans)
trajectory_prob = ODEProblem(trajectory_eqs!, u0, tspan)
trajectory_sol = solve(trajectory_prob, RK4(), dt = 0.01)
push!(trajectories, trajectory_sol.u)
end
return trajectories
end;FINALLY: The Plotting Function
The final function makes a plot with Plotly
Code
function final_plot(;
fps,
As,
contour_xs,
contour_ys,
contour_f_xy,
contour_g_xy,
slope_start_xys,
slope_end_xys,
trajectories,
)
annotations = []
for (fp, A) ∈ zip(fps, As)
classification = classify_jacobian(A)
annotation = attr(x = fp[1], y = fp[2], text = "<b>$classification</b>")
push!(annotations, annotation)
end
traces::Vector{GenericTrace} = []
trace_fxy = contour(
x = contour_xs,
y = contour_ys,
z = contour_f_xy',
contours_start = 0,
contours_end = 0,
contours_coloring = "lines",
colorscale = [[0, "gold"], [1.0, "white"]],
line = attr(width = 2),
name = "f(x,y)",
showlegend = false,
)
push!(traces, trace_fxy)
trace_gxy = contour(
x = contour_xs,
y = contour_ys,
z = contour_g_xy',
contours_start = 0,
contours_end = 0,
contours_coloring = "lines",
colorscale = [[0, "darkorange"], [1.0, "white"]],
line = attr(width = 2),
name = "g(x,y)",
showlegend = false,
)
push!(traces, trace_gxy)
for (i, (start_xy, end_xy)) ∈ enumerate(zip(slope_start_xys, slope_end_xys))
showlegend = i == 1
trace_slope = scatter(
x = [start_xy[1], end_xy[1]],
y = [start_xy[2], end_xy[2]],
mode = "lines",
line = attr(color = "lawngreen"),
name = "slope",
showlegend = showlegend,
)
push!(traces, trace_slope)
end
for (i, trajectory) ∈ enumerate(trajectories)
showlegend = i == 1
trace_start = scatter(
x = [trajectory[1][1]],
y = [trajectory[1][2]],
mode = "markers",
marker = attr(color = "blue", size = 10),
name = "start",
showlegend = showlegend,
)
trace_trajectory = scatter(
x = [x for (x, _) ∈ trajectory],
y = [y for (_, y) ∈ trajectory],
mode = "lines",
line = attr(color = "black"),
name = "trajectory",
showlegend = showlegend,
)
trace_end = scatter(
x = [trajectory[end][1]],
y = [trajectory[end][2]],
mode = "markers",
marker = attr(color = "red", size = 10),
name = "end",
showlegend = showlegend,
)
push!(traces, trace_start)
push!(traces, trace_trajectory)
push!(traces, trace_end)
end
for (i, (fp, A)) ∈ enumerate(zip(fps, As))
classification = classify_jacobian(A)
size = 15
color = "darkorchid"
symbol = classification == "Saddle" ? "circle-open" : "circle"
showlegend = i == 1
trace_fp = scatter(
x = [fp[1]],
y = [fp[2]],
mode = "markers",
marker = attr(color = color, symbol = symbol, size = size),
showlegend = showlegend,
name = "Fixed Points",
)
push!(traces, trace_fp)
end
plot_bgcolor = "white"
paper_bgcolor = "white"
border_width = 1
gridwidth = 1
border_color = "black"
gridcolor = "lightgray"
layout = Layout(
width = 550,
height = 500,
plot_bgcolor = plot_bgcolor,
paper_bgcolor = paper_bgcolor,
xaxis = attr(
showline = true,
linewidth = border_width,
linecolor = border_color,
mirror = true,
showgrid = true,
gridcolor = gridcolor,
gridwidth = gridwidth,
),
yaxis = attr(
showline = true,
linewidth = border_width,
linecolor = border_color,
mirror = true,
showgrid = true,
gridcolor = gridcolor,
gridwidth = gridwidth,
),
annotations = annotations,
)
return plot(traces, layout)
end;Example 6.3.1
Here is a phase portrait for the following system of equations:
\[ \dot x = -x + x^3 \] \[ \dot y = -2y \]
Code
function example_6_3_1()
# Find fixed points
fixed_points_system_of_eqs(u, p) = SA[-u[1]+u[1]^3, -2*u[2]]
fps = find_fixed_points(
fixed_points_system_of_eqs;
guess_xs = range(-2.0, 2.0, 10),
guess_ys = range(-2.0, 2.0, 10),
)
println(fps)
# Find Jacobians
forward_diff_system_of_eqs(u) = [-u[1] + u[1]^3, -2*u[2]]
As = find_jacobians(forward_diff_system_of_eqs, fps)
println(As)
# Find contours to plot nullclines
min_x, max_x = -1.5, 1.5
min_y, max_y = -1.0, 1.0
f(u::Union{Vector{Float64},Tuple{Float64,Float64}}) = -u[1] + u[1]^3
g(u::Union{Vector{Float64},Tuple{Float64,Float64}}) = -2*u[2]
contour_xs = range(min_x, max_x, 100)
contour_ys = range(min_y, max_y, 100)
contour_f_xy, contour_g_xy = nullcline_contours(f, g, contour_xs, contour_ys)
start_xys, end_xys = slope_field(f, g, range(min_x, max_x, 10), range(min_y, max_y, 10))
# Compute trajectories
function trajectory_eqs!(du, u, p, t)
du[1] = -u[1] + u[1]^3
du[2] = -2*u[2]
end
u0s = [
[-0.833, -0.556],
[-1.167, -0.556],
[-1.167, 0.556],
[-0.833, 0.556],
[0.833, -0.556],
[1.167, -0.556],
[1.167, 0.556],
[0.833, 0.556],
[0.0, 0.8],
[0.0, -0.8],
[-0.4, 0.0],
[0.4, 0.0],
]
tspans = [
(-0.5, 0.5),
(-0.1, 0.2),
(-0.1, 0.2),
(-0.5, 0.5),
(-0.5, 0.5),
(-0.1, 0.2),
(-0.1, 0.2),
(-0.5, 0.5),
(-0.5, 0.125),
(-0.5, 0.125),
(-0.5, 0.5),
(-0.5, 0.5),
]
trajectories = calculate_trajectories(trajectory_eqs!, u0s, tspans)
# Create final plot
return final_plot(;
fps = fps,
As = As,
contour_xs = contour_xs,
contour_ys = contour_ys,
contour_f_xy = contour_f_xy,
contour_g_xy = contour_g_xy,
slope_start_xys = start_xys,
slope_end_xys = end_xys,
trajectories = trajectories,
)
end
example_6_3_1()Lotka-Volterra Competition
On Page 172 of Strogatz, the following system of equations that model two species vying for the same food resource, where is x(t) is the population of rabbits and y(t) the population of sheep.
\[ \dot x = x(3-x-2y) \] \[ \dot y = y(2-x-y) \]
Code
function lotka_volterra_competition()
# Find fixed points
system_of_eqs_01(u, p) = SA[u[1]*(3-u[1]-2*u[2]), u[2]*(2-u[1]-u[2])]
fps = find_fixed_points(
system_of_eqs_01;
guess_xs = range(0.0, 3.25, 10),
guess_ys = range(0.0, 2.25, 10),
)
# println(fps)
# Find Jacobians
system_of_eqs_02(u) = [u[1]*(3-u[1]-2*u[2]), u[2]*(2-u[1]-u[2])]
As = find_jacobians(system_of_eqs_02, fps)
# println(As)
# Find contours to plot nullclines
min_x, max_x = 0.0, 3.0
min_y, max_y = 0.0, 2.0
f(u::Union{Vector{Float64},Tuple{Float64,Float64}}) = u[1]*(3-u[1]-2*u[2])
g(u::Union{Vector{Float64},Tuple{Float64,Float64}}) = u[2]*(2-u[1]-u[2])
contour_xs = range(min_x, max_x, 100)
contour_ys = range(min_y, max_y, 100)
contour_f_xy, contour_g_xy = nullcline_contours(f, g, contour_xs, contour_ys)
start_xys, end_xys = slope_field(f, g, range(min_x, max_x, 10), range(min_y, max_y, 10))
# Compute trajectories
function trajectory_eqs!(du, u, p, t)
du[1] = u[1]*(3-u[1]-2*u[2])
du[2] = u[2]*(2-u[1]-u[2])
end
u0s = [
[0.1, 0.222],
[0.333, 0.222],
[0.333, 0.333],
[0.0, 0.222],
[0.133, 0.156],
[2.67, 2.0],
[2.333, 2.0],
[3.0, 2.0],
[3.0, 1.78],
[3.0, 1.56],
[0.333, 0.444],
[0.667, 0.222],
[1.0, 0.222],
[0.333, 0.667],
[0.333, 0.889],
[0.333, 1.111],
[2.0, 2.0],
[1.667, 2.0],
[1.333, 2.0],
[1.0, 2.0],
[0.667, 2.0],
[0.333, 2.0],
[0.333, 0.0],
[3.0, 1.333],
[3.0, 1.111],
[3.0, 0.889],
[3.0, 0.667],
[3.0, 0.444],
]
tspans = [
(-1.0, 2.0),
(-1.0, 2.0),
(-1.0, 2.0),
(-1.0, 1.5),
(-1.0, 2.0),
(-1.0, 0.0),
(-1.0, 0.0),
(-1.0, 0.0),
(-1.0, 0.0),
(-1.0, 0.0),
(-1.0, 2.0),
(-1.0, 2.0),
(-1.0, 2.0),
(-1.0, 4.0),
(-1.0, 4.0),
(-1.0, 4.0),
(-1.0, 0.0),
(-1.0, 0.0),
(-1.0, 0.0),
(-1.0, 0.0),
(-1.0, 0.0),
(-1.0, 0.0),
(-1.0, 0.75),
(-1.0, 0.0),
(-1.0, 0.0),
(-1.0, 0.0),
(-1.0, 0.0),
(-1.0, 0.0),
]
trajectories = calculate_trajectories(trajectory_eqs!, u0s, tspans)
# Create final plot
return final_plot(;
fps = fps,
As = As,
contour_xs = contour_xs,
contour_ys = contour_ys,
contour_f_xy = contour_f_xy,
contour_g_xy = contour_g_xy,
slope_start_xys = start_xys,
slope_end_xys = end_xys,
trajectories = trajectories,
)
end
lotka_volterra_competition()Example 6.5.2
Code
function example_6_5_2()
# Find fixed points
system_of_eqs_01(u, p) = SA[u[2], u[1]-u[1]^3]
fps = find_fixed_points(
system_of_eqs_01;
guess_xs = range(-1.0, 1.0, 3),
guess_ys = range(-0.1, 0.1, 3),
)
println(fps)
# Find Jacobians
system_of_eqs_02(u) = [u[2], u[1]-u[1]^3]
As = find_jacobians(system_of_eqs_02, fps)
println(As)
# Find contours to plot nullclines
min_x, max_x = -1.5, 1.5
min_y, max_y = -1.0, 1.0
f(u::Union{Vector{Float64},Tuple{Float64,Float64}}) = u[2]
g(u::Union{Vector{Float64},Tuple{Float64,Float64}}) = u[1]-u[1]^3
contour_xs = range(min_x, max_x, 100)
contour_ys = range(min_y, max_y, 100)
contour_f_xy, contour_g_xy = nullcline_contours(f, g, contour_xs, contour_ys)
start_xys, end_xys = slope_field(f, g, range(min_x, max_x, 10), range(min_y, max_y, 10))
# Compute trajectories
function trajectory_eqs!(du, u, p, t)
du[1] = u[2]
du[2] = u[1]-u[1]^3
end
u0s = [
[1.1, 1.1],
[0.167, 0.333],
[0.5, 0.111],
[-0.5, 0.111],
]
tspans = [
(0.0, 7.5),
(0.0, 11.75),
(0.0, 5.0),
(0.0, 5.0),
]
trajectories = calculate_trajectories(trajectory_eqs!, u0s, tspans)
# Create final plot
return final_plot(;
fps = fps,
As = As,
contour_xs = contour_xs,
contour_ys = contour_ys,
contour_f_xy = contour_f_xy,
contour_g_xy = contour_g_xy,
slope_start_xys = start_xys,
slope_end_xys = end_xys,
trajectories = trajectories,
)
end
example_6_5_2()