Julia: Strogatz Example 7.1.2

This document builds on Strogatz’s Example 7.1.2 of the van der Pol equation.

The following Julia packages will be needed for this document:

using NonlinearSolve
using DifferentialEquations
using SciMLBase
using StaticArrays
using ForwardDiff
using LinearAlgebra
using PlotlyJS

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.

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.

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.

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

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.

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.

function calculate_trajectories(trajectory_eqs!, u0s, tspans, ps)
    trajectories = []
    for (u0, tspan, p) ∈ zip(u0s, tspans, ps)
        trajectory_prob = ODEProblem(trajectory_eqs!, u0, tspan, p)
        trajectory_sol = solve(trajectory_prob, RK4(), dt = 0.01)
        push!(trajectories, (trajectory_sol.u, trajectory_sol.t))
    end
    return trajectories
end;

FINALLY: The Plotting Function

The final function makes a plot with Plotly.

function final_plot(;
    title,
    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} = []
    x_line_traces::Vector{GenericTrace} = []
    y_line_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, trajectory_ts)) ∈ 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,
        )
        x_line_trace = scatter(
            x = trajectory_ts, 
            y = [x for (x, _) ∈ trajectory],
            mode = "lines",
            line = attr(color = "darkslateblue"),
            name = "X Component",
        )
        y_line_trace = scatter(
            x = trajectory_ts, 
            y = [y for (_, y) ∈ trajectory],
            mode = "lines",
            line = attr(color = "crimson"),
            name = "Y Component",
        )
        push!(traces, trace_start)
        push!(traces, trace_trajectory)
        push!(traces, trace_end)
        push!(x_line_traces, x_line_trace)
        push!(y_line_traces, y_line_trace)
    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"
    fig = make_subplots(
        rows = 2,
        cols = 2,
        vertical_spacing = 0.1,
        subplot_titles = ["Phase Portrait" "Y Trajectory"],
    )
    for trace ∈ traces
        add_trace!(fig, trace, row = 1, col = 1)
    end
    for x_line_trace ∈ x_line_traces
        add_trace!(fig, x_line_trace, row = 1, col = 2)
    end
    for y_line_trace ∈ y_line_traces
        add_trace!(fig, y_line_trace, row = 2, col = 2)
    end
    relayout!(
        fig,
        title = title,
        width = 800,
        height = 700,
        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(
            # domain = [phase_portrait_min_y, phase_portrait_max_y],
            showline = true,
            linewidth = border_width,
            linecolor = border_color,
            mirror = true,
            showgrid = true,
            gridcolor = gridcolor,
            gridwidth = gridwidth,
        ),
        xaxis2 = attr(
            title = "<b>t</b>",
            showline = true,
            linewidth = border_width,
            linecolor = border_color,
            mirror = true,
            showgrid = true,
            gridcolor = gridcolor,
            gridwidth = gridwidth,
        ),
        yaxis2 = attr(
            title = "<b>x</b>",
            showline = true,
            linewidth = border_width,
            linecolor = border_color,
            mirror = true,
            showgrid = true,
            gridcolor = gridcolor,
            gridwidth = gridwidth,
        ),
        xaxis4 = attr(
            title = "<b>t</b>",
            showline = true,
            linewidth = border_width,
            linecolor = border_color,
            mirror = true,
            showgrid = true,
            gridcolor = gridcolor,
            gridwidth = gridwidth,
        ),
        yaxis4 = attr(
            title = "<b>y</b>",
            showline = true,
            linewidth = border_width,
            linecolor = border_color,
            mirror = true,
            showgrid = true,
            gridcolor = gridcolor,
            gridwidth = gridwidth,
        ),
        annotations = annotations,
    )
    return fig
end;

Function to Assemble the Plots

function ex_7_1_2(μ)
    # Find fixed points
    eqs_01(u, p) = SA[u[2], u[2]*(1-u[1]^2)-u[1]]
    fps = find_fixed_points(
        eqs_01;
        guess_xs = range(-3.5, 3.5, 10),
        guess_ys = range(-3.5, 3.5, 10),
    )
    println(fps)

    # Find Jacobians
    eqs_02(u) = [u[2], u[2]*(1-u[1]^2)-u[1]]
    As = find_jacobians(eqs_02, fps)
    println(As)

    # Find contours to plot nullclines and slope field
    min_x, max_x = -3.0, 3.0
    min_y, max_y = -3.0, 3.0
    f(u::Union{Vector{Float64},Tuple{Float64,Float64}}) = u[2]
    g(u::Union{Vector{Float64},Tuple{Float64,Float64}}) = u[2]*(1-u[1]^2)-u[1]
    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] = p[1]*(1-u[1]^2)*u[2]-u[1]
    end

    u0s = [[-0.444, 0.444]]
    tspans = [(0.0, 20.0)]
    ps = [[μ]]

    trajectories = calculate_trajectories(trajectory_eqs!, u0s, tspans, ps)

    # Create final plot
    return final_plot(;
        title = "<b>μ=$μ</b>",
        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 7.1.2

This example explores the van der Pol equation:

\[ \ddot{x} + \mu(x^2-1)\dot{x} + x \]

To solve this second-order equation with an RK4 algorithm for visualization, I split it into a system of first-order equations:

\[ {dx_1 \over dt} = x_2 \]

\[ {dx_2 \over dt} = -\mu(x^2-1)x_2-x \]

Then I plotted the phase portrait with the trajectories on it alongside the x and y components versus t of the trajectories separately.

Below are phase portraits for differing values of mu.

ex_7_1_2(0.5)

Figure 1

ex_7_1_2(1.0)

Figure 2

ex_7_1_2(1.5)

Figure 3

ex_7_1_2(2.0)

Figure 4

ex_7_1_2(3.0)

Figure 5