x_eq(t, x0, a) = x0 * exp(a * t)
y_eq(t, y0) = y0 * exp(-t);Phase Portrait Inspiration
Chapter 5 in Strogatz, 3rd edition has many great phase portraits! In this post, I seek to make similar phase portraits in Julia to enhance my understanding of the material.
In function definitions below, the trailing
;prevents extraneous output in this document.
Fig. 5.1.5-Inspired
There are three phase portraits from the discussion of Figure 5.1.5 that I wanted to display. I set initial conditions at cartesian-coordinate conversions of points around the unit circle to create these plots. As shown in the text, the system being plotted is:
\[ x(t) = x_0 e^{at} \] \[ y(t) = y_0 e^{-t} \]
Which are conveniently expressed in Julia code:
The function to draw a selection of the phase portraits found in Figure 5.1.5 is the following:
function portrait(a::Float64, r::Float64, width::Int64 = 500, height::Int64 = 500)
traces::Vector{GenericTrace} = []
angles = [0.0, π/4, π/2, π, 3π/4, 5π/4, 3π/2, 7π/4]
x0s = [r * cos(θ) for θ ∈ angles]
y0s = [r * sin(θ) for θ ∈ angles]
ts = range(0.0, 2.0, length = 10)
for (i, (x0, y0)) ∈ enumerate(zip(x0s, y0s))
xs = x_eq.(ts, x0, a)
ys = y_eq.(ts, y0)
showlegend = i == 1
trace_start = scatter(
x = [x0],
y = [y0],
mode = "markers",
marker = attr(color = "blue", size = 10),
name = "start",
showlegend = showlegend,
)
trace_line = scatter(
x = xs,
y = ys,
mode = "lines",
line = attr(color = "black"),
name = "path",
showlegend = showlegend,
)
trace_end = scatter(
x = [xs[end]],
y = [ys[end]],
mode = "markers",
marker = attr(color = "red", size = 10),
name = "end",
showlegend = showlegend,
)
push!(traces, trace_start)
push!(traces, trace_line)
push!(traces, trace_end)
end
title = "<b>a = $a</b>"
plot_bgcolor = "white"
paper_bgcolor = "white"
border_width = 1
gridwidth = 1
border_color = "black"
gridcolor = "lightgray"
layout = Layout(
plot_bgcolor = plot_bgcolor,
paper_bgcolor = paper_bgcolor,
title = title,
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,
),
width = width,
height = height,
)
plot(traces, layout)
end;The text shows plots for varying values of a, which are reflected in the plots below. The blue dots are at the beginnings of the paths and the red dots are at the end of the paths to give a sense of directionality to the plots.
The plots are made with Plotly, so they can be hovered and zoomed, too!
portrait(-2.0, 0.1, 550, 500)portrait(-1.0, 1.0, 550, 500)portrait(1.0, 1.0, 550, 500)Fig 5.2.2-Inspired
Leading up to Figure 5.2.2, we are taken through solving the following:
\[ \begin{pmatrix} \dot x \\ \dot y \end{pmatrix} = \begin{pmatrix} 1 & 1 \\ 4 & -2 \end{pmatrix} \begin{pmatrix} x \\ y \end{pmatrix} \]
As shown in the book, the solution to this system is:
\[ x(t) = e^{2t} + e^{-3t} \] \[ y(t) = e^{2t} - 4e^{-3t} \]
Like the book, I will call the 2x2 matrices in the equation above A. The text’s solution solves for initial conditions (x0, y0) = (2, 3), but in the examples below I solve for different initial conditions to get many different curves.
The following functions takes a given A matrix and vector of initial conditions, computes eigensolutions, and solves for c1 and c2 and returns the appropriate x(t) and y(t) functions to plot the phase portraits.
function solve_for_ics(A::Matrix{Float64}, ics::Vector{Float64})
eig = eigen(A)
λ = eig.values
v = eig.vectors
c = eig.vectors \ ics
x_eq(t::Float64) = c[1]*v[1, 1]*exp(λ[1]*t) + c[2]*v[1, 2]*exp(λ[2]*t)
y_eq(t::Float64) = c[1]*v[2, 1]*exp(λ[1]*t) + c[2]*v[2, 2]*exp(λ[2]*t)
x_eq, y_eq
end;The following plot is very similar to portrait() above, except it is designed to accomodate the A matrix and call solve_for_ics():
function portrait02(
A::Matrix{Float64},
r::Float64,
ts::Vector{Float64},
width::Int64 = 500,
height::Int64 = 500,
)
angles = [0.0, π/4, π/2, π, 3π/4, 5π/4, 3π/2, 7π/4]
x0s = [r * cos(θ) for θ ∈ angles]
y0s = [r * sin(θ) for θ ∈ angles]
traces::Vector{GenericTrace} = []
for (i, (x0, y0)) ∈ enumerate(zip(x0s, y0s))
x_eq, y_eq = solve_for_ics(A, [x0, y0])
xs = x_eq.(ts)
ys = y_eq.(ts)
showlegend = i == 1
trace_start = scatter(
x = [xs[1]],
y = [ys[1]],
mode = "markers",
marker = attr(color = "blue", size = 10),
name = "start",
showlegend = showlegend,
)
trace_line = scatter(
x = xs,
y = ys,
mode = "lines",
marker = attr(color = "black"),
name = "path",
showlegend = showlegend,
)
trace_end = scatter(
x = [xs[end]],
y = [ys[end]],
mode = "markers",
marker = attr(color = "red", size = 10),
name = "stop",
showlegend = showlegend,
)
push!(traces, trace_start)
push!(traces, trace_line)
push!(traces, trace_end)
end
title = "<b>A = $(string(A))</b>"
plot_bgcolor = "white"
paper_bgcolor = "white"
border_width = 1
gridwidth = 1
border_color = "black"
gridcolor = "lightgray"
layout = Layout(
plot_bgcolor = plot_bgcolor,
paper_bgcolor = paper_bgcolor,
title = title,
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,
),
width = width,
height = height,
)
plot(traces, layout)
end;As before, the blue dots show beginings of paths and the red dots show the end.
portrait02([1.0 1.0; 4.0 -2.0], 1.0, collect(range(-0.75, 0.75, 100)), 550, 500)portrait02([2.0 2.0; 3.0 -3.0], 1.0, collect(range(-0.5, 0.5, 100)), 550, 500)Problem and Solution 5.2.2-Inspired
Problem 5.2.2 seeks to solve
\[ \begin{pmatrix} \dot x \\ \dot y \end{pmatrix} = \begin{pmatrix} 1 & -1 \\ 1 & 1 \end{pmatrix} \begin{pmatrix} x \\ y \end{pmatrix} \]
What makes this interesting is that it results in complex eigenvalues. In the back of the book, the solution is given as:
\[ \mathbf x(t) = c_1 e^t \begin{pmatrix} \cos t \\ \sin t \end{pmatrix} + c_2e^t \begin{pmatrix} -\sin t \\ \cos t \end{pmatrix} \]
However, no further guidance was given on how to find c1 and c2, so in orer to make the nifty spiral plots promised by such a system, I needed to find another approach. Enter DifferentialEquations.jl, the Julia massive Julia package for solving DEs. With this package in hand, the following code made quick work of the problem:
function solve_system(
A::Matrix{Float64},
x0::Vector{Float64},
tspan::Tuple{Float64,Float64},
)
f(x, p, t) = A * x
prob = ODEProblem(f, x0, tspan)
solve(prob, Tsit5())
endsolve_system (generic function with 1 method)The following function, very similar to those above, can plot such a system.
function portrait03(
A::Matrix{Float64},
tspan::Tuple{Float64,Float64},
rs::Vector{Float64},
angles::Vector{Float64},
width::Int64,
height::Int64,
)
x0s = [[r * cos(θ), r * sin(θ)] for θ ∈ angles, r ∈ rs]
traces::Vector{GenericTrace} = []
for (i, x0) ∈ enumerate(x0s)
sol = solve_system(A, x0, tspan)
xs = sol[1, :]
ys = sol[2, :]
showlegend = i == 1
trace_path = scatter(
x = xs,
y = ys,
mode = "lines",
line = attr(color = "black"),
name = "path",
showlegend = showlegend,
)
trace_start = scatter(
x = [xs[1]],
y = [ys[1]],
mode = "markers",
line = attr(color = "blue"),
name = "start",
showlegend = showlegend,
)
trace_end = scatter(
x = [xs[end]],
y = [ys[end]],
mode = "markers",
line = attr(color = "red"),
name = "end",
showlegend = showlegend,
)
push!(traces, trace_path)
push!(traces, trace_start)
push!(traces, trace_end)
end
title = "<b>A = $A</b>"
plot_bgcolor = "white"
paper_bgcolor = "white"
border_width = 1
gridwidth = 1
border_color = "black"
gridcolor = "lightgray"
layout = Layout(
title = title,
plot_bgcolor = plot_bgcolor,
paper_bgcolor = paper_bgcolor,
width = width,
height = height,
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,
),
)
plot(traces, layout)
endportrait03 (generic function with 1 method)Here is the plot of a solution to Problem 5.2.2
portrait03(
[1.0 -1.0; 1.0 1.0],
(0.0, 2.0),
[0.5],
[0.0, π/4, π/2, π, 3π/4, 5π/4, 3π/2, 7π/4],
550,
500,
)And a a plot of the solution to following system:
\[ \begin{pmatrix} \dot x \\ \dot y \end{pmatrix} = \begin{pmatrix} 3 & -3 \\ 0 & 2 \end{pmatrix} \begin{pmatrix} x \\ y \end{pmatrix} \]
portrait03(
[3.0 -3.0; 2.0 2.0],
(0.0, 2.0),
[0.5],
[0.0, π/4, π/2, π, 3π/4, 5π/4, 3π/2, 7π/4],
550,
500,
)Figure 5.3.2-inspired
Figure 5.3.2 represents the following system:
\[ \begin{pmatrix} \dot x \\ \dot y \end{pmatrix} = \begin{pmatrix} a & b \\ b & a \end{pmatrix} \begin{pmatrix} x \\ y \end{pmatrix} \]
Assuming:
\[ a < 0, b > 0 \]
To solve the systems, I use the solve_for_ics() function I use earlier since the eigenvalues are real-valued. To plot these phase portraits I use the following function, similar to those above, but with a few tweaks to make plotting these phase portraits easier:
function portrait(
A::Matrix{Float64},
angles::Vector{Float64},
rs::Vector{Float64},
ts::Vector{Float64},
width::Int64 = 550,
height::Int64 = 500,
)
x0s = [r * cos(θ) for θ ∈ angles, r ∈ rs]
y0s = [r * sin(θ) for θ ∈ angles, r ∈ rs]
traces::Vector{GenericTrace} = []
# for (i, (x0, y0)) ∈ enumerate(Base.product(x0s, y0s))
for (i, (x0, y0)) ∈ enumerate(zip(x0s, y0s))
x_eq, y_eq = solve_for_ics(A, [x0, y0])
xs = x_eq.(ts)
ys = y_eq.(ts)
showlegend = i == 1
trace_start = scatter(
x = [xs[1]],
y = [ys[1]],
mode = "markers",
marker = attr(color = "blue", size = 10),
name = "start",
showlegend = showlegend,
)
trace_line = scatter(
x = xs,
y = ys,
mode = "lines",
marker = attr(color = "black"),
name = "path",
showlegend = showlegend,
)
trace_end = scatter(
x = [xs[end]],
y = [ys[end]],
mode = "markers",
marker = attr(color = "red", size = 10),
name = "stop",
showlegend = showlegend,
)
push!(traces, trace_start)
push!(traces, trace_line)
push!(traces, trace_end)
end
title = "<b>A = $(string(A))</b>"
plot_bgcolor = "white"
paper_bgcolor = "white"
border_width = 1
gridwidth = 1
border_color = "black"
gridcolor = "lightgray"
layout = Layout(
plot_bgcolor = plot_bgcolor,
paper_bgcolor = paper_bgcolor,
title = title,
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,
),
width = width,
height = height,
)
plot(traces, layout)
end;(Actually, in the book this represents a Romeo and Juliet story, but I am switching back to x and y for the purposes of this demo.)
The first phase portrait has a = -2 and b = 1 to match the condition a2 > b2:
portrait(
[-2.0 1.0; 1.0 -2.0],
[0.0, π/2, π/4, 3π/4, π, 5π/4, 3π/2, 7π/4],
[1.0],
collect(range(-1.0, 1.0, 100)),
)The second phase portrait has a = -1 and b = 2 to match the condition a2 < b2:
portrait(
[-1.0 2.0; 2.0 -1.0],
[0.0, π/4, π/3, 2π/3, 3π/4, 5π/6, 7π/6, 5π/4, 4π/3, 5π/3, 7π/4, 11π/6],
[1.0],
collect(range(-0.5, 0.5, 100)),
)