Mini-Lecture 6-32 - Linear differential equations (continuous time) - The general case.ipynb

Kernel: SageMath 9.5

from itertools import product

import numpy as np
import plotly.graph_objects
import plotly.figure_factory
from ipywidgets import IntSlider, SelectionSlider, Checkbox, Text, \
    interactive_output, VBox, HBox, HTMLMath, Layout

mydefault = plotly.graph_objects.layout.Template()
mydefault.layout.xaxis.showgrid = False
mydefault.layout.yaxis.showgrid = False
mydefault.layout.xaxis.showline = True
mydefault.layout.yaxis.showline = True
mydefault.layout.yaxis.linewidth = 2
mydefault.layout.xaxis.ticks = "outside"
mydefault.layout.yaxis.ticks = "outside"
mydefault.layout.hovermode = False
mydefault.layout.dragmode = "pan"
mydefault.layout.scene.hovermode = False
mydefault.layout.xaxis.showspikes = False
mydefault.layout.yaxis.showspikes = False
mydefault.layout.scene.xaxis.showspikes = False
mydefault.layout.scene.yaxis.showspikes = False
mydefault.layout.scene.zaxis.showspikes = False
plotly.io.templates["mydefault"] = mydefault
plotly.io.templates.default = "mydefault"

Learning goals:

• Be able to classify the type of equilibrium point (at the origin) for any linear differential equation.

• Be able to identify the dominant eigenvalue of any linear differential equation, and be able to explain what this mean's about the model's behavior.

The set-up:

We start with a linear differential equation model: $$\begin{bmatrix} X' \\ Y' \end{bmatrix} = M \begin{bmatrix} X \\ Y \end{bmatrix}$$

Use the eigenvectors of $M$ to form a basis of $\mathbb{R}^2$, and use this to define a new coordinate system, which we call $R,S$-coordinates, as usual. When we re-write the system of differential equations in this new coordinate system, the matrix turns into a diagonal matrix, with the eigenvalues of $M$ on its diagonal:

$$\begin{bmatrix} R' \\ S' \end{bmatrix} = \begin{pmatrix} \lambda_1 & 0 \\ 0 & \lambda_2 \end{pmatrix} \begin{bmatrix} R \\ S \end{bmatrix}$$

def plot_graph_paper(T, **options):
    color1, color2 = options.get("colors", ("red", "green"))
    v1_color = colors[color1].darker(0.25)
    v2_color = colors[color2].darker(0.25)
    grid1_color = colors[color1].lighter(0.75)
    grid2_color = colors[color2].lighter(0.75)
    axis1_label, axis2_label = options.get("axes_labels", ("$R$", "$S$"))
    show_axes = options.get("show_axes", True)
    gridlines1, gridlines2 = options.get("gridlines", (6, 6))
    v1 = T.column(0)
    v2 = T.column(1)
    xymax = min(g / a for g, a in zip((gridlines1, gridlines2), ~T*vector((1,1))) if a != 0)
    xmax, ymax = options.get("xyrange", (xymax, xymax))
    p = Graphics()
    if options.get("show_grid", True):
        for n in range(-gridlines1, gridlines1 + 1):
            if n == 0 and show_axes:
                continue
            p += line((v1*n - gridlines2*v2, v1*n + gridlines2*v2), color=grid2_color, thickness=0.5)
        for n in range(-gridlines2, gridlines2 + 1):
            if n == 0 and show_axes:
                continue
            p += line((v2*n - gridlines1*v1, v2*n + gridlines1*v1), color=grid1_color, thickness=0.5)
    if show_axes:
        p += line((-gridlines1*v1, gridlines1*v1), color=color1, thickness=1.5)
        p += line((-gridlines2*v2, gridlines2*v2), color=color2, thickness=1.5)
        p += text(axis1_label, 1.1*min(abs(xmax / v1[0]), abs(ymax / v1[1]))*v1, color=color1, fontsize=16)
        p += text(axis2_label, 1.1*min(abs(xmax / v2[0]), abs(ymax / v2[1]))*v2, color=color2, fontsize=16)
    if options.get("show_vectors", True):
        p += plot(v1, color=v1_color, thickness=2, zorder=5)
        p += plot(v2, color=v2_color, thickness=2, zorder=5)
    p.set_axes_range(xmin=-xmax, xmax=xmax, ymin=-ymax, ymax=ymax)
    return p

def initchoice(a, vmax):
    if a < 0:
        return randint(round(0.5*vmax + 1), vmax)
    if a > 0:
        return randint(1, round(0.5*vmax - 1))
    return randint(1, vmax - 1)

def nondiagonal_interactive(T, eigenvalues, **options):
    initial_xy = options.get("initial_state", None)
    state_vars = options.get("state_vars", ("X", "Y"))
    rsmax = 10
    t_range = srange(0, 15.01, 0.05)
    lambda1, lambda2 = eigenvalues
    D = diagonal_matrix(RDF, eigenvalues)
    M = matrix(RDF, T*D*~T).round(4)
    P1 = T*diagonal_matrix((1,0))*~T
    P2 = T*diagonal_matrix((0,1))*~T
    vectorfield_rs(R, S) = tuple(D*vector((R, S)))
    vectorfield_xy(X, Y) = tuple(M*vector((X, Y)))
    scale = 0.1/max(abs(lambda1), abs(lambda2))
    if initial_xy is None:
        initial_rs = vector([initchoice(lambda1, rsmax), initchoice(lambda2, rsmax)])
        initial_xy = vector([round(x,3) for x in T * initial_rs])
    else:
        initial_rs = vector([round(x, 3) for x in ~T * vector(RDF, initial_xy)])
    solution_rs = solution_xy = []

    dont_update = False
    def update(xyt, rst, show_eigen, show_grid, showaxisfield_xy, showaxisfield_rs, 
               showwholefield_xy, showwholefield_rs, axis1max, axis2max):
        if dont_update or len(solution_rs) == 0:
            return
        rsmax = ceil(axis1max / abs(np.array(T).flatten()).min())
        p1  = plot_graph_paper(T, show_vectors=show_eigen and not showaxisfield_xy, 
                               show_axes=show_eigen, show_grid=show_grid, 
                               gridlines=(rsmax,rsmax), xyrange=(axis1max,axis1max))
        p1 += list_plot(solution_xy[:xyt+1], plotjoined=True, color="blue", zorder=5)
        p1 += point(solution_xy[xyt], size=30, color="black", zorder=5)
        p2  = list_plot(solution_rs[:rst+1], plotjoined=True, color="blue", zorder=5)
        p2 += point(solution_rs[rst], size=30, color="black", zorder=5)
        rslist = srange(-0.95*axis1max, axis1max, 0.1*axis1max)
        if showaxisfield_xy:
            for r in rslist:
                p1 += arrow2d(T*vector((r,0)), T*vector((r,0)) + scale*T*vectorfield_rs(r,0), 
                              color="darkred", width=3, arrowsize=2, zorder=4)
            for s in rslist:
                p1 += arrow2d(T*vector((0,s)), T*vector((0,s)) + scale*T*vectorfield_rs(0,s), 
                              color="darkgreen", width=3, arrowsize=2, zorder=4)
        if showwholefield_xy:
            for r, s in product(rslist, rslist):
                p1 += arrow2d(T*vector((r,s)), T*vector((r,s)) + scale*T*vectorfield_rs(r,s), 
                              color="orange", width=1, arrowsize=2)
        elif show_grid:
            mypnt = vector(solution_xy[xyt])
            proj1 = P1*mypnt
            proj2 = P2*mypnt
            p1 += line((proj1, mypnt, proj2), color="lightgray")
            if not showaxisfield_xy:
                p1 += arrow2d(mypnt, mypnt + vectorfield_xy(*mypnt), arrowsize=2, color="orange")
                p1 += arrow2d(proj1, proj1 + vectorfield_xy(*proj1), arrowsize=2, width=3, zorder=4, color="red")
                p1 += arrow2d(proj2, proj2 + vectorfield_xy(*proj2), arrowsize=2, width=3, zorder=4, color="green")
        rslist = srange(-0.95*axis2max, axis2max, 0.1*axis2max)
        if showaxisfield_rs:
            for r in rslist:
                p2 += arrow2d(vector((r,0)), vector((r,0)) + scale*vectorfield_rs(r,0), 
                              color="red", width=1, arrowsize=2, zorder=4)
            for s in rslist:
                p2 += arrow2d(vector((0,s)), vector((0,s)) + scale*vectorfield_rs(0,s), 
                              color="green", width=1, arrowsize=2, zorder=4)
        if showwholefield_rs:
            for r, s in product(rslist, rslist):
                p2 += arrow2d(vector((r,s)), vector((r,s)) + scale*vectorfield_rs(r,s), 
                              color="orange", width=1, arrowsize=2)
        elif show_grid:
            mypnt = solution_rs[rst]
            proj1 = mypnt * (1,0)
            proj2 = mypnt * (0,1)
            p2 += line((proj1, mypnt, proj2), color="lightgray")
            if not showaxisfield_rs:
                p2 += arrow2d(mypnt, mypnt + vectorfield_rs(*mypnt), arrowsize=2, color="orange")
                p2 += arrow2d(proj1, proj1 + vectorfield_rs(*proj1), arrowsize=2, zorder=4, color="red")
                p2 += arrow2d(proj2, proj2 + vectorfield_rs(*proj2), arrowsize=2, zorder=4, color="green")
        p1.set_axes_range(xmin=-axis1max, xmax=axis1max, ymin=-axis1max, ymax=axis1max)
        p2.set_axes_range(xmin=-axis2max, xmax=axis2max, ymin=-axis2max, ymax=axis2max)
        p1.axes_labels(["$%s$" % label for label in state_vars])
        p2.axes_labels(("$R$", "$S$"))
        p1.set_aspect_ratio(1)
        p2.set_aspect_ratio(1)
        both = multi_graphics([(p1, (0, 0.25, 0.57, 0.57)), (p2, (0.43, 0.25, 0.57, 0.57))])
        both.show(figsize=12)

    def change_initial_rs(change):
        nonlocal dont_update, solution_rs, solution_xy
        initial_rs = change["new"].strip("()[]").split(",")
        try:
            initial_rs = vector([float(x) for x in initial_rs])
        except:
            return
        initial_xy = T * initial_rs
        solution_rs = desolve_odeint(vectorfield_rs, initial_rs, t_range, [R, S])
        solution_xy = desolve_odeint(vectorfield_xy, initial_xy, t_range, [X, Y])
        dont_update = True
        rst.value = 0
        xyt.value = 1
        dont_update = False
        xyt.value = 0

    def change_initial_xy(change):
        initial_xy = change["new"].strip("()[]").split(",")
        try:
            initial_xy = vector([float(x) for x in initial_xy])
        except:
            return
        state_rs.value = str(vector([round(x, 3) for x in ~T * initial_xy]))

    p2b = str.maketrans("()", "[]")
    label = ((r"$ \begin{bmatrix} %s' \\ %s' \end{bmatrix} = " % state_vars) + 
             latex(M) + (r"\begin{bmatrix} %s \\ %s \end{bmatrix} \qquad " % state_vars) + 
             r"\vec{v}_1 = " + latex(T[:,0]).translate(p2b) + r" \qquad " + 
             (r"\text{with } \lambda_1 = %.2f \qquad" % lambda1) + 
             r"\vec{v}_2 = " + latex(T[:,1]).translate(p2b) + r" \qquad " + 
             (r"\text{with } \lambda_2 = %.2f $" % lambda2) )
    equation = HTMLMath(label)
    layout = Layout(width="250px", margin="0px 0px 0px 0px")
    t_values = {round(t,2): n for n, t in enumerate(t_range)}
    xyt  = SelectionSlider(options=t_values, layout=layout, description="$t$ for left sol'n:")
    rst  = SelectionSlider(options=t_values, layout=layout, description="$t$ for right sol'n:")
    showaxisfield_xy = Checkbox(value=False, layout=layout, description="Show vec fld on axes")
    showaxisfield_rs = Checkbox(value=False, layout=layout, description="Show vec fld on axes")
    showwholefield_xy = Checkbox(value=False, layout=layout, description="Show whole vec fld")
    showwholefield_rs = Checkbox(value=False, layout=layout, description="Show whole vec fld")
    show_eigen = Checkbox(value=False, layout=layout, description="Show eigenlines")
    show_grid = Checkbox(value=False, layout=layout, description="Show gridlines")
    zoom1 = SelectionSlider(options=[1, 2, 5, 10, 20, 50, 100, 200, 500, 1000], value=10, 
                            layout=layout, description="Zoom (left):")
    zoom2 = SelectionSlider(options=[1, 2, 5, 10, 20, 50, 100, 200, 500, 1000], value=10, 
                            layout=layout, description="Zoom (right):")
    state_rs = Text(continuous_update=False, layout=layout, description="Initial (R,S):")
    state_xy = Text(value=str(initial_xy), continuous_update=False, 
                    layout=layout, description="Initial (%s,%s):" % state_vars)
    state_rs.observe(change_initial_rs, "value")
    state_xy.observe(change_initial_xy, "value")
    state_rs.value = str(initial_rs)
    output = interactive_output(update, dict(
            xyt=xyt, rst=rst, show_eigen=show_eigen, show_grid=show_grid, 
            showaxisfield_xy=showaxisfield_xy, showaxisfield_rs=showaxisfield_rs, 
            showwholefield_xy=showwholefield_xy, showwholefield_rs=showwholefield_rs, 
            axis1max=zoom1, axis2max=zoom2))
    controls1 = HBox((state_xy, state_rs, show_eigen, show_grid))
    controls2 = HBox((showaxisfield_xy, showwholefield_xy, showaxisfield_rs, showwholefield_rs))
    controls3 = HBox((xyt, zoom1, rst, zoom2))
    display(VBox((equation, controls1, controls2, controls3, output)))

v1 = vector((2, 1))
v2 = vector((3, -1))
T = column_matrix((v1, v2))
nondiagonal_interactive(T, (0.2, -0.6), initial_state=(9,-2))

M = matrix(RDF, [[0.36, 0.35], [0.24, 0.82]])
M0 = matrix(ZZ, 100*M)
show(M0)
for evalue, evectors, mult in M0.eigenvectors_right():
    print(evalue, evectors[0])

96 (1, 12/7)
22 (1, -2/5)

v1 = vector((1, 1))
v2 = vector((1, -2))
T = column_matrix((v1, v2))
nondiagonal_interactive(T, (1.5, 0.6), initial_state=(0,3))

v1 = vector((-1, 1))
v2 = vector((2, 1))
T = column_matrix((v1, v2))
nondiagonal_interactive(T, (-0.2, -0.8), initial_state=(9,9))

def plotly_streamline(M, x_range, y_range, **options):
    options = options.copy()
    plotpoints = options.pop("plotpoints", 100)
    try:
        plotpointsx, plotpointsy = plotpoints
    except:
        plotpointsx = plotpointsy = plotpoints
    color = options.pop("color", "red")
    options.setdefault("line_color", color)
    x, xmin, xmax = x_range
    y, ymin, ymax = y_range
    x = np.linspace(float(xmin), float(xmax), plotpointsx)
    y = np.linspace(float(ymin), float(ymax), plotpointsy)
    X, Y = np.meshgrid(x, y)
    (a,b), (c,d) = M.rows()
    u = a*X + b*Y
    v = c*X + d*Y
    streamline = plotly.figure_factory.create_streamline(x=x, y=y, u=u, v=v, **options)
    return streamline.data[0]

def plotly_vector_field(M, x_range, y_range, **options):
    options = options.copy()
    plotpoints = options.pop("plotpoints", 20)
    try:
        plotpointsx, plotpointsy = plotpoints
    except:
        plotpointsx = plotpointsy = plotpoints
    color = options.pop("color", "red")
    options.setdefault("line_color", color)
    x, xmin, xmax = x_range
    y, ymin, ymax = y_range
    diagonal = sqrt( ((xmax - xmin)^2/plotpointsx^2 + (ymax - ymin)^2/plotpointsy^2) )
    x = np.linspace(float(xmin), float(xmax), plotpointsx)
    y = np.linspace(float(ymin), float(ymax), plotpointsy)
    x, y = np.meshgrid(x, y)
    (a,b), (c,d) = M.rows()
    u = a*x + b*y
    v = c*x + d*y
    scale = diagonal / np.linalg.norm(np.array((u.flatten(), v.flatten())), axis=0).max()
    options.setdefault("scale", float(scale))
    quiver = plotly.figure_factory.create_quiver(x=x, y=y, u=u, v=v, **options)
    return quiver.data[0]

def complex_interactive():
    state_vars = var("X, Y")
    t_range = srange(0, 50.01, 0.1)

    fig = plotly.graph_objects.FigureWidget()
    fig.layout.margin = dict(b=10, t=10, l=10, r=10)
    fig.layout.xaxis.domain = (0, 0.65)
    fig.layout.yaxis.domain = (0, 1)
    fig.layout.xaxis.range = (-10, 10)
    fig.layout.yaxis.range = (-10, 10)
    fig.layout.xaxis.title.text = "$X$"
    fig.layout.yaxis.title.text = "$Y$"
    fig.layout.xaxis.showline = False
    fig.layout.yaxis.showline = False
    fig.layout.showlegend = False
    fig.add_scatter(mode="lines", line_color="limegreen")
    field = fig.data[-1]
    fig.add_scatter(mode="lines", line_color="red", line_smoothing=1.3)
    streamlines = fig.data[-1]
    fig.add_annotation(showarrow=False, font_size=16, x=0.68, y=0.6, 
                       xref="paper", yref="paper", xanchor="left", yanchor="middle")
    fig.add_annotation(showarrow=False, font_size=16, x=0.68, y=0.4, 
                       xref="paper", yref="paper", xanchor="left", yanchor="middle")
    equation, eigenvalues = fig.layout.annotations

    def update(a, b, other, ccw):
        T = matrix([[1, other], [0, (-1)^ccw]])
        R = matrix([[a, b], [-b, a]])
        M = matrix(RDF, T*R*~T).round(4)
        if a > 0:
            initial_state = (1, 0)
        elif a < 0: 
            initial_state = (10, 0)
        else:
            initial_state = (6, 0)
        vectorfield(X, Y) = tuple(M*vector((X, Y)))
        solution = desolve_odeint(vectorfield, initial_state, t_range, [X, Y])
        trace = plotly_vector_field(M, (X, -10, 10), (Y, -10, 10))
        label1 = ((r"$ \begin{bmatrix} %s' \\ %s' \end{bmatrix} = " % state_vars) + 
                 latex(M) + (r"\begin{bmatrix} %s \\ %s \end{bmatrix} $" % state_vars))
        label2 = r"$ \text{Eigenvalues: } %s \pm %s i $" % (a, b)
        with fig.batch_update():
            field.x = trace.x
            field.y = trace.y
            streamlines.x = solution[:,0]
            streamlines.y = solution[:,1]
            equation.text = label1
            eigenvalues.text = label2

    a_slider = slider([round(x,1) for x in srange(-2, 2.05, 0.1)], default=0, label="Real part:")
    b_slider = slider([round(x,1) for x in srange(0.1, 2.05, 0.1)], default=2, label="Imag part:")
    T_slider = slider(-2, 2, 1, default=0, label="Roundness?")
    ccw = checkbox(False, label="CCW?")
    output = interactive_output(update, dict(a=a_slider, b=b_slider, other=T_slider, ccw=ccw))
    display(VBox((HBox((a_slider, b_slider, T_slider, ccw)), output)))
    return fig

complex_interactive()

One explanation for the above conclusion:

Recall that the solution to $R' = \lambda R$ is:     $R(t) = R(0) e^{\lambda t}$.

What on earth is $e^{\lambda t}$ if $\lambda$ is a complex number?

Theorem: (Euler's Formula) $$e^{i x} = \cos(x) + i \sin(x)$$

As a result of this, if $\lambda = a + bi$, then $$\begin{aligned} e^{\lambda t} &= e^{(a + bi)t} \\ &= e^{at + ibt} \\ &= e^{at} \cdot e^{ibt} \\ &= e^{at} \big( \cos(bt) + i\sin(bt) \big) \end{aligned}$$

Conclusions:

• Just as we did with discrete-time linear models, for a linear differential equation model, once you look at the state space using the $R,S$ coordinates*, the system behaves exactly like a diagonal (decoupled) linear model. That is: $$\begin{bmatrix} X' \\ Y' \end{bmatrix} = \begin{pmatrix} \dotsm & \dotsm \\ \dotsm & \dotsm \end{pmatrix} \begin{bmatrix} X \\ Y \end{bmatrix} \implies \begin{bmatrix} R' \\ S' \end{bmatrix} = \begin{pmatrix} \lambda_1 & 0 \\ 0 & \lambda_2 \end{pmatrix} \begin{bmatrix} R' \\ S' \end{bmatrix}$$

• Therefore, everything we just learned about linear differential equations with diagonal matrices (in particular, the type of equilibrium point at the origin, and the concept of the dominant eigenvalue) now carries over to any linear differential equation whose matrix is diagonalizable.

• For linear differential equations, complex eigenvalues once again cause rotating/spiraling behavior. Spiraling inward/outward depends on whether the eigenvalues are to the right or left of $0$. That is:

  • If the real part of the eigenvalue is negative, trajectories will spiral inward.

  • If the real part of the eigenvalue is positive, trajectories will spiral outward.