Mini-Lecture 6-30 - Complex eigenvalues.ipynb

Kernel: SageMath 9.3

import numpy as np
import plotly.graph_objects

from ipywidgets import interactive_output, VBox, HBox, HTMLMath

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.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"

def round_complex(z, digits):
    if z.imag_part():
        return CDF(round(z.real_part(), digits) + round(z.imag_part(), digits) * I)
    return RDF(round(z, digits))

def latex_vector(v, round=None):
    if round is None:
        result = column_matrix(v)
    else:
        result = column_matrix(RDF, v).round(round)
    result = latex(result).replace("left(", "left[").replace("right)", "right]")
    return LatexExpr(result)

def latex_discrete_linear(M, state_vars=None):
    if state_vars is None:
        state_vars = ("X", "Y", "Z")[:M.nrows()]
    state_vector = latex_vector(var([X + "_t" for X in state_vars]))
    result = state_vector.replace("_{t}", "_{t+1}")
    return result + " = " + latex(M) + state_vector

def latex_eigenvalues(M, digits=4, show_abs=False):
    eigenvalues = []
    for l in sorted(M.eigenvalues(), key=lambda l: -abs(l)):
        l0 = round_complex(l, digits)
        if l.imag_part() < 0:
            continue
        if l.imag_part():
            text = r"{} \pm {}i".format(l0.real(), l0.imag())
            if show_abs:
                text += r" \text{ (abs. value %s) }" % round(abs(l0), digits)
            eigenvalues.append(text)
        else:
            eigenvalues.append(str(l0))
    return ", ".join(eigenvalues)

def bracket(start, end, offset, label=None, **options):
    options = options.copy()
    fontsize = options.pop("fontsize", 14)
    rotation = matrix(2, (0, 1, -1, 0))
    unit = (vector(end) - vector(start)).normalized()
    offset = offset*rotation*unit
    p1 = vector(start) + offset
    p2 = vector(end) + offset
    p  = line((p1 - 0.5*offset, p1 + 0.5*offset), **options)
    p += line((p2 - 0.5*offset, p2 + 0.5*offset), **options)
    p += line((p1, p2), linestyle="dashed", **options)
    if label:
        rotation = atan2(unit[1], unit[0]) * 180 / pi
        p += text(label, 0.5*p1 + 0.5*p2 + offset, fontsize=fontsize, #rotation=rotation, 
                  **options)
    return p

Learning goals:

• Know that eigenvalues can be complex numbers, and that these still matter in determining the behavior of a model.

• Know what complex eigenvalues tell us about the behavior of a discrete-time linear model.

• Know how to compute the absolute value of a complex number, and why this is significant.

Example:

$$\begin{bmatrix} X_{t+1} \\ Y_{t+1} \end{bmatrix} = \begin{pmatrix} 2 & -1 \\ 5 & 4 \end{pmatrix} \begin{bmatrix} X_t \\ Y_t \end{bmatrix}$$

Characteristic polynomial: $\qquad \lambda^2 - 6 \lambda + 13$

Roots (by the quadratic formula): $$\begin{aligned} \lambda &= \frac{6 \pm \sqrt{(-6)^2 - 4 \cdot 13}}{2} \\ &= \frac{6 \pm \sqrt{-16}}{2} = 3 \pm \frac{4 \sqrt{-1}}{2} = 3 \pm 2i \end{aligned}$$

Review:

Definitions: A complex number is a number of the form $a + bi$, where $a$ and $b$ can be any real numbers, and $i$ represents $\sqrt{-1}$. (A number of the form $bi$ by itself is called an imaginary number.)

In a complex number $a + bi$, the number $a$ is called the real part of the complex number, and the number $b$ is called the imaginary part.

Example: In the complex number $z = 2 - 7i$, the real part of $z$ is $2$, and the imaginary part of $z$ is $-7$. (Not $-7i$!)

def discrete_linear_interactive3d(M, initial_state, **options):
    t_max = options.get("t_max", 100)
    show_abs = options.get("show_abs", False)
    axes_labels = options.get("axes_labels", ("X", "Y", "Z"))
    xmin = options.get("xmin", 0)
    ymin = options.get("ymin", 0)
    zmin = options.get("zmin", 0)
    xmax = options.get("xmax", 10)
    ymax = options.get("ymax", 10)
    zmax = options.get("zmax", 10)

    solution = [vector(initial_state)]
    for t in range(t_max):
        solution.append(M * solution[-1])
    solution = np.array(solution, dtype=float)

    fig = plotly.graph_objects.FigureWidget()
    fig.layout.margin = dict(t=10, b=10, r=0, l=0)
    fig.layout.showlegend = False
    fig.layout.scene.xaxis.title.text = axes_labels[0]
    fig.layout.scene.yaxis.title.text = axes_labels[1]
    fig.layout.scene.zaxis.title.text = axes_labels[2]
    fig.layout.scene.xaxis.range = (xmin, xmax)
    fig.layout.scene.yaxis.range = (ymin, ymax)
    fig.layout.scene.zaxis.range = (zmin, zmax)
    fig.layout.scene.aspectmode = "cube"
    fig.add_scatter3d(mode="markers+lines", marker_color="black", marker_size=2, 
                      line_color="blue")
    trajectory = fig.data[-1]
    eigenlines = []
    for evalue, evectors, mult in M.eigenvectors_right():
        if evalue.imag_part():
            continue
        evector = vector(RDF, evectors[0])
        end = min(np.array((xmax, ymax, zmax)) / evector) * evector
        start = min(np.array((xmin, ymin, zmin)) / evector) * evector
        x, y, z = np.array((start, end)).transpose()
        fig.add_scatter3d(x=x, y=y, z=z, mode="lines", line_color="red")
        eigenlines.append(fig.data[-1])

    def update(t, show_eigenlines):
        with fig.batch_update():
            trajectory.x = solution[:t+1,0]
            trajectory.y = solution[:t+1,1]
            trajectory.z = solution[:t+1,2]
            for eigenline in eigenlines:
                eigenline.visible = show_eigenlines

    def update_text(change=None):
        label = latex_discrete_linear(M)
        if show_eigenvals.value:
            label += r"\qquad \text{Eigenvalues: } " + latex_eigenvalues(M, show_abs=show_abs)
        text.value = "${}$".format(label)

    t = slider(0, t_max, 1, label="$t$")
    show_eigenlines = checkbox(False, label="Show eigenline(s)")
    show_eigenvals = checkbox(False, label="Show eigenvalues")
    output = interactive_output(update, dict(t=t, show_eigenlines=show_eigenlines))
    controls = HBox((t, show_eigenlines, show_eigenvals))
    text = HTMLMath()
    show_eigenvals.observe(update_text, "value")
    update_text()
    display(controls, text, output)

    return fig

# Lionfish example from lab:
M = matrix(RDF, [
    [0      , 0    , 35315    ], 
    [0.00003, 0.777,     0    ], 
    [0      , 0.071,     0.949], 
])
# Locust example from textbook:
M = matrix(RDF, [
    [0   , 0   , 1000], 
    [0.02, 0   ,    0], 
    [0   , 0.06,    0], 
])
M = matrix(RDF, [
    [0.05, 0   , 65   ], 
    [0.30, 0.03,  0   ], 
    [0   , 0.05,  0.05], 
])
M = matrix(RDF, [
    [0   , 0   , 1000], 
    [0.02, 0   ,    0], 
    [0   , 0.06,    0], 
])
for evalue, evectors, mult in M.eigenvectors_right():
    print(round_complex(evalue, 4), round(abs(evalue), 4))
    print("         ", evectors[0])

-0.5313 + 0.9203*I 1.0627
          (0.9998223729551143, -0.009408688754320866 - 0.016296326955085667*I, -0.000531234906140622 + 0.0009201258481896411*I)
-0.5313 - 0.9203*I 1.0627
          (0.9998223729551143, -0.009408688754320866 + 0.016296326955085667*I, -0.000531234906140622 - 0.0009201258481896411*I)
1.0627 1.0627
          (-0.9998223729551141, -0.018817377508641736, -0.0010624698122812446)

M = matrix(RDF, [
    [0   , 0   , 1000   ], 
    [0.02, 0   ,    0   ], 
    [0   , 0.05,    0.10], 
])
state = vector((500, 300, 15))
discrete_linear_interactive3d(M, state, xmax=80000, ymax=2000, zmax=100, 
                              axes_labels=("Eggs", "Hoppers", "Adults"))

Where complex numbers usually come from:

Complex numbers arise naturally as roots of polynomials, e.g. from the quadratic formula:

Theorem: The roots (solutions) of $Ax^2 + Bx + C = 0$ are $$x = \frac{-B \pm \sqrt{B^2 - 4AC}}{2A}$$

So if $B^2 - 4AC$ is negative, both roots will be complex numbers. Furthermore, they will have the form $$p + qi \qquad \text{and} \qquad p - qi .$$

Review: (maybe?)

Definition: Given a complex number $a + bi$, its complex conjugate is the number $a - bi$.

A complex number and its conjugate always have the same real part, but their imaginary parts are negatives of each other.

Theorem: Whenever a polynomial has a complex number as a root, its complex conjugate will also be a root. In other words, as roots of polynomials, complex numbers always occur in complex conjugate pairs.

def rotation_matrix(scale, theta=None, field=RDF):
    if theta is None:
        a, b = scale.real_part(), scale.imag_part()
        return matrix(field, 2, (a, -b, b, a))
    a, b = cos(theta), sin(theta)
    return scale * matrix(field, 2, (a, -b, b, a))

def discrete_linear_interactive(M, initial_state, **options):
    options = options.copy()
    t_max = options.pop("t_max", 100)
    show_abs = options.pop("show_abs", False)
    options.setdefault("axes_labels", ("$X$", "$Y$"))
    options.setdefault("xmin", -10)
    options.setdefault("ymin", -10)
    options.setdefault("xmax",  10)
    options.setdefault("ymax",  10)
    options.setdefault("figsize", 7.5)

    solution = [vector(initial_state)]
    for t in range(t_max):
        solution.append(M * solution[-1])

    def update(t, plotjoined, show_vectors):
        p = list_plot(solution[:t+1], color="black", size=30)
        if plotjoined:
            p += list_plot(solution[:t+1], plotjoined=True, color="gray")
        if show_vectors:
            for v in solution[:t+1]:
                p += plot(v, color="gray", arrowsize=2)
        p.show(**options)

    def update_text(change=None):
        label = latex_discrete_linear(M)
        label = latex_discrete_linear(M)
        if show_eigenvals.value:
            label += r"\qquad \text{Eigenvalues: } " + latex_eigenvalues(M, show_abs=show_abs)
        text.value = "${}$".format(label)

    t = slider(0, t_max, 1, label="$t$")
    plotjoined = checkbox(False, label="Connect dots")
    show_vectors = checkbox(False, label="Show vectors")
    show_eigenvals = checkbox(False, label="Show eigenvalues")
    output = interactive_output(update, dict(t=t, plotjoined=plotjoined, show_vectors=show_vectors))
    controls = HBox((t, plotjoined, show_vectors, show_eigenvals))
    text = HTMLMath()
    show_eigenvals.observe(update_text, "value")
    update_text()
    display(controls, text, output)

M = rotation_matrix(0.8 - 0.6j)
discrete_linear_interactive(M, (2,9), aspect_ratio=1)

M = rotation_matrix(0.96 + 0.28j)
discrete_linear_interactive(M, (7,2), aspect_ratio=1)

M = rotation_matrix(0.912 - 0.266j)
discrete_linear_interactive(M, (-4,7), aspect_ratio=1)

M = rotation_matrix(1.008 + 0.294j)
discrete_linear_interactive(M, (1,1), aspect_ratio=1)

def discrete_linear_magnitude_interactive(M, initial_state, **options):
    options = options.copy()
    t_max = options.pop("t_max", 60)
    show_abs = options.pop("show_abs", False)
    options.setdefault("axes_labels", ("$X$", "$Y$"))
    options.setdefault("xmin", -10)
    options.setdefault("ymin", -10)
    options.setdefault("xmax",  10)
    options.setdefault("ymax",  10)
    options.setdefault("aspect_ratio", 1)

    solution = [vector(initial_state)]
    for t in range(t_max):
        solution.append(M * solution[-1])
    ymax = max([v.norm() for v in solution])

    def update(t, show_vectors, show_trend):
        p1 = list_plot(solution[:t+1], color="black", size=30, **options)
        p2 = Graphics()
        if show_vectors:
            for t0, v in enumerate(solution[:t+1]):
                p1 += plot(v, color="gray", arrowsize=2)
                p2 += plot(vector((0, v.norm())), start=(t0,0), color="gray", arrowsize=2)
        if show_trend:
            c = solution[0].norm()
            r = abs(M.eigenvalues()[0])
            f(t1) = c*r^t1
            p2 += plot(f, (t1, 0, t_max), color="red")
        p2.set_axes_range(xmax=t_max, ymax=ymax)
        p2.axes_labels(("$t$", "Length of vectors"))
        both = multi_graphics([(p1, (0,0.27,0.46,0.46)), (p2, (0.54,0.27,0.46,0.46))])
        both.show(figsize=12)


    def update_text(change=None):
        label = latex_discrete_linear(M)
        label = latex_discrete_linear(M)
        if show_eigenvals.value:
            label += r"\qquad \text{Eigenvalues: } " + latex_eigenvalues(M, show_abs=show_abs)
        text.value = "${}$".format(label)

    t = slider(0, t_max, 1, label="$t$")
    show_vectors = checkbox(False, label="Show vectors")
    show_trend = checkbox(False, label="Show trend")
    show_eigenvals = checkbox(False, label="Show eigenvalues")
    output = interactive_output(update, dict(t=t, show_vectors=show_vectors, show_trend=show_trend))
    controls = HBox((t, show_vectors, show_trend, show_eigenvals))
    text = HTMLMath()
    show_eigenvals.observe(update_text, "value")
    update_text()
    display(controls, text, output)

M = rotation_matrix(0.912 - 0.266j)
discrete_linear_magnitude_interactive(M, (-4,7))

M = rotation_matrix(1.008 + 0.294j)
discrete_linear_magnitude_interactive(M, (1,1), aspect_ratio=1)

@interact(x=slider(-10, 10, 1, label="Number:"), show_abs=checkbox(False, "Show absolute value"))
def real_line(x, show_abs):
    p = text("0", (0,-0.6), color="black", zorder=4)
    p += polygon(((-2,0.05), (1,0.05), (1,1.4), (-2,1.4)), color="white", zorder=3)
    p += polygon(((-2,-0.1), (1,-0.1), (1,-1.4), (-2,-1.4)), color="white", zorder=3)
    p += point((x,0), size=40, color="blue", zorder=5)
    if show_abs and x > 0:
        p += bracket((0,0), (x,0), 0.5, label=str(abs(x)), color="blue", zorder=5)
    elif show_abs and x < 0:
        p += bracket((x,0), (0,0), 0.5, label=str(abs(x)), color="blue", zorder=5)
    p.show(xmin=-10, xmax=10, ymin=-1.3, ymax=1.3, aspect_ratio=1, axes_labels=("Real", ""))

@interact(x=slider(-10, 10, 1, label="Real part:"), y=slider(-10, 10, 1, label="Imag. part:"), 
          show_abs=checkbox(False, "Show absolute value"), 
          show_legs=checkbox(False, "Show real and imag. parts"))
def complex_plane(x, y, show_abs, show_legs):
    p = point((x,y), size=40, color="purple")
    p += text("${} + {}i$".format(x, y), (x,y), color="purple", fontsize=14, 
              vertical_alignment="bottom", horizontal_alignment="left")
    if show_abs:
        p += line(((0,0), (x,y)), linestyle="dashed", color="purple")
        p += text("${}$".format(latex(sqrt(x^2 + y^2))), (x/2, y/2), color="purple", fontsize=14, 
                  vertical_alignment="bottom", horizontal_alignment="right")
    if show_legs:
        p += line(((0,0), (x,0)), linestyle="dashed", color="blue", thickness=2)
        p += text("${}$".format(x), (x/2, 0), color="blue", fontsize=14, 
                  vertical_alignment="top", horizontal_alignment="center")
        p += line(((x,0), (x,y)), linestyle="dashed", color="red")
        p += text("${}$".format(y), (x, y/2), color="red", fontsize=14, 
                  vertical_alignment="center", horizontal_alignment="left")
    p.show(xmin=-10, xmax=10, ymin=-10, ymax=10, aspect_ratio=1, axes_labels=("Real", "Imag"))

The absolute value of a complex number:

Definition: For a complex number $z = a + bi$, the absolute value of $z$ is defined to be $$|z| = \sqrt{a^2 + b^2}$$

Note that for its complex conjugate: $$|a - bi| = \sqrt{a^2 + (-b)^2} = \sqrt{a^2 + b^2}$$

M = rotation_matrix(0.96 + 0.28j)
discrete_linear_interactive(M, (8,3), aspect_ratio=1, show_abs=True)

M = rotation_matrix(0.912 - 0.266j)
discrete_linear_interactive(M, (-4,7), aspect_ratio=1, show_abs=True)

M = rotation_matrix(1.008 + 0.294j)
discrete_linear_interactive(M, (1,1), aspect_ratio=1, show_abs=True)

M = matrix(RDF, [
    [0   , 0   , 1000   ], 
    [0.02, 0   ,    0   ], 
    [0   , 0.05,    0.10], 
])
state = vector((500, 300, 15))
discrete_linear_interactive3d(M, state, xmax=80000, ymax=2000, zmax=100, 
                              axes_labels=("Eggs", "Hoppers", "Adults"), show_abs=True)

Conclusions:

1. For a discrete-time linear model, eigenvalues that are complex numbers are associated with rotating (or spiraling in or out) behavior.

2. The absolute value of a complex number $a + bi$ is defined as $\sqrt{a^2 + b^2}$.

3. More specifically:

   • A pair of complex eigenvalues with absoute value $> 1$ means the solution will spiral outward (rotating combined with exponential growth).

   • A pair of complex eigenvalues with absoute value $< 1$ means the solution will spiral inward (rotating combined with exponential growth).

   • A pair of complex eigenvalues with absoute value exactly equal to $1$ means the solution will oscillate.
     (neutral equilibrium point/center, neither exponential growth nor decay. This is a degenerate case, i.e. it rarely occurs in realistic models.)