Mini-Lecture 6-35 - Linear dynamics in 3D.ipynb

Kernel: SageMath 9.2

from itertools import product

import numpy as np
import plotly.graph_objects

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

def plotly_line(points, **options):
    "Draw a line through the given points"
    useroptions = options
    options = dict(mode="lines", hoverinfo="none")
    options.update(useroptions)
    options.setdefault("line_color", options.pop("color", "blue"))
    x, y, z = np.array(points, dtype=float).transpose()
    return plotly.graph_objects.Scatter3d(x=x, y=y, z=z, **options)

def plotly_plane(v, w, axismax=10, **options):
    "Draw a plane through the origin spanned by v and w"
    useroptions = options
    options = dict(showscale=False, hoverinfo="none", contours_x_highlight=False, 
                   contours_y_highlight=False, contours_z_highlight=False)
    options.update(useroptions)
    color = options.pop("color", "lightblue")
    options.setdefault("colorscale", [[0, color], [1, color]])
    u1 = v.normalized()
    u2 = (w - (w*u1)*u1).normalized()
    u = np.array((-axismax, axismax), dtype=float)
    u, v = np.meshgrid(u, u)
    x, y, z = [float(a)*u + float(b)*v for a, b in zip(u1, u2)]
    return plotly.graph_objects.Surface(x=x, y=y, z=z, **options)

def plotly_text(text, location, **options):
    useroptions = options
    options = dict(showarrow=False)
    if useroptions.get("paper", False):
        options.update(xanchor="left", yanchor="bottom")
    options.update(useroptions)
    options.setdefault("font_color", options.pop("color", "black"))
    if options.pop("paper", False):
        options.update(xref="paper", yref="paper")
    x, y = np.array(location, dtype=float)
    return plotly.graph_objects.layout.Annotation(text=text, x=x, y=y, **options)

def plotly_add_text3d(fig, text, location, **options):
    useroptions = options
    options = dict(showarrow=False)
    options.update(useroptions)
    options.setdefault("font_color", options.pop("color", "black"))
    x, y, z = np.array(location, dtype=float)
    annotation = dict(text=text, x=x, y=y, z=z, **options)
    fig.layout.scene.annotations = fig.layout.scene.annotations + (annotation,)

def plotly_trajectory3d(trajectory, **options):
    useroptions = options
    options = dict(mode="lines", hoverinfo="none", legendgroup=hex(randint(0, 1 << 32))[2:])
    options.update(useroptions)
    options.setdefault("line_color", options.pop("color", "blue"))
    arrows = options.pop("arrows", 6)
    try:
        arrows = np.array(list(arrows), dtype=int)
    except:
        arclength = np.cumsum(np.linalg.norm(trajectory[1:] - trajectory[:-1], axis=1))
        arrows_so_far = np.rint(arrows*arclength / arclength[-1])
        arrows = np.where(arrows_so_far[1:] - arrows_so_far[:-1])[0]
    p1 = plotly.graph_objects.Scatter3d(x=trajectory[:,0], y=trajectory[:,1], z=trajectory[:,2], **options)

    del options["mode"]
    options["showlegend"] = False
    color = options.pop("line_color")
    options.setdefault("colorscale", [[0, color], [1, color]])
    options.setdefault("showscale", False)
    options.setdefault("sizeref", 0.33) # Adjust this for smaller/largers cones
    cone_pos = trajectory[arrows]
    cone_dir = trajectory[arrows+1] - cone_pos
    cone_dir /= np.linalg.norm(cone_dir, axis=1)[:,np.newaxis]
    p2 = plotly.graph_objects.Cone(x=cone_pos[:,0], y=cone_pos[:,1], z=cone_pos[:,2], 
                                   u=cone_dir[:,0], v=cone_dir[:,1], w=cone_dir[:,2], **options)
    return (p1, p2)

def latex_vector(v, round=None):
    if isinstance(v[0], str):
        return (r"\left[ \begin{array}{c} " + r" \\ ".join(v) 
                + r" \end{array} \right]")
    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_linear_discrete(M, state_vars=None):
    if state_vars is None:
        state_vars = ("X", "Y", "Z")[:M.nrows()]
    current_state = latex_vector([fr"{X}_t" for X in state_vars])
    next_state = latex_vector([fr"{X}_{{t+1}}" for X in state_vars])
    return fr"{next_state} = {latex(M)} {current_state}"

def latex_linear_diffeq(M, state_vars=None):
    if state_vars is None:
        state_vars = ("X", "Y", "Z")[:M.nrows()]
    current_state = latex_vector([str(X) for X in state_vars])
    change_vector = latex_vector([fr"{X}'" for X in state_vars])
    return fr"{change_vector} = {latex(M)} {current_state}"

def latex_complex(z, digits=4):
    a, b = float(round(z.real(), digits)), float(round(z.imag(), digits))
    return fr"{a:g} \pm {b:g} i"

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 = fr"{l0.real()} \pm {l0.imag()} i"
            if show_abs:
                text += fr" \text{{ (abs. value {round(abs(l0), digits)}) }}"
            eigenvalues.append(text)
        else:
            eigenvalues.append(str(l0))
    return ", ".join(eigenvalues)

QQI = GaussianIntegers().fraction_field()

def block_diagonal(eigenvalues):
    blocks = []
    for evalue in eigenvalues:
        a, b = evalue.real(), evalue.imag()
        if b:
            blocks.append(matrix(((a, b), (-b, a))))
        else:
            blocks.append(matrix(((a,),)))
    return block_diagonal_matrix(*blocks)

def undiagonalize(eigenvalues, basis):
    D = block_diagonal(eigenvalues)
    T = column_matrix(basis)
    return T*D*~T

def diagonalize(M):
    basis = []
    eigenvalues = []
    for evalue, evectors, mult in M.eigenvectors_right():
        if evalue.imag() < 0:
            continue
        if evalue.imag():
            for evector in evectors:
                eigenvalues.append(evalue)
                u = vector([a.real() for a in evector])
                v = vector([a.imag() for a in evector])
                basis.append(u * u.denominator())
                basis.append(v * v.denominator())
        else:
            for v in evectors:
                eigenvalues.append(evalue)
                basis.append(v * v.denominator())
    return eigenvalues, basis

def colorchoice(evalue):
    if abs(evalue) < 1:
        color = "green"
    elif abs(evalue) > 1:
        color = "coral"
    else:
        color = "orange"
    return color

def colorchoice_continuous(value):
    return colorchoice(e^value)

def linear_discrete_3d_interactive(eigenvalues, basis, **options):
    pass

def linear_diffeq_3d_interactive(eigenvalues, basis, **options):
    colorchoice = colorchoice_continuous
    axismax = options.get("axismax", 10)
    state_vars = options.get("state_vars", "XYZ")
    other_vars = options.get("other_vars", "QRS")
    field = options.get("field", QQ)
    dominant_rate = max([l.real() for l in eigenvalues])
    t_max = float(log(2*axismax) / abs(dominant_rate))
    t_range = np.linspace(0, t_max, 1000, dtype=float)
    D = block_diagonal(eigenvalues)
    T = column_matrix(basis)
    M = matrix(field, T*D*~T)
    vectorfield(X, Y, Z) = tuple(M*vector((X, Y, Z)))

    fig = plotly.graph_objects.Figure()
    fig.layout.margin = dict(t=10, b=10, r=10, l=10)
    fig.layout.legend.y = 0.9
    fig.layout.scene.domain.x = (0.3, 1)
    fig.layout.scene.xaxis.title.text = state_vars[0]
    fig.layout.scene.yaxis.title.text = state_vars[1]
    fig.layout.scene.zaxis.title.text = state_vars[2]
    fig.layout.scene.xaxis.range = (-axismax, axismax)
    fig.layout.scene.yaxis.range = (-axismax, axismax)
    fig.layout.scene.zaxis.range = (-axismax, axismax)
    fig.layout.scene.aspectmode = "cube"
    fig.layout.scene.hovermode = "closest"

    text = {0.8: f"${latex_linear_diffeq(M)}$", 0.72: r"$\text{Eigenvalues:}$"}
    y = 0.72
    axis_info = {}
    initial_states = []
    i = 0
    for evalue in eigenvalues:
        color = colorchoice(evalue)
        if evalue.imag():
            # Draw the eigenplane:
            v, w = basis[i:i+2]
            scale = axismax / max(v.norm(), w.norm())
            fig.add_trace(plotly_plane(v, w, axismax=2*axismax, color="light" + color, 
                                       opacity=0.5, showlegend=False))
            fig.add_trace(plotly_line(((0,0,0), 0.9*scale*v), color=color, 
                                      line_width=4, showlegend=False))
            fig.add_trace(plotly_line(((0,0,0), 0.9*scale*w), color=color, 
                                      line_width=4, showlegend=False))
            plotly_add_text3d(fig, f"${other_vars[i]}$", 0.9*scale*v, color=color, 
                              font_size=16, yshift=10)
            plotly_add_text3d(fig, f"${other_vars[i+1]}$", 0.9*scale*w, color=color, 
                              font_size=16, yshift=10)
            # Choose some initial states:
            u = scale * (v + w)
            if evalue.real() > 0:
                axis_info[f"{other_vars[i:i+2]} plane"] = (color, [0.1*u])
                initial_states.append((0.1 * axismax,))
                initial_states.append(initial_states[-1])
            elif evalue.real() < 0:
                axis_info[f"{other_vars[i:i+2]} plane"] = (color, [u])
                initial_states.append((0.6 * axismax,))
                initial_states.append(initial_states[-1])
            else:
                axis_info[f"{other_vars[i:i+2]} plane"] = (color, [0.3*u, 0.6*u, 0.9*u])
                initial_states.append((0.2 * axismax, 0.4 * axismax, 0.6 * axismax))
                initial_states.append((0.5*axismax,))
            # Add the text for the eigenvalue and eigenvectors:
            y -= 0.15
            text[y] = (fr"$\lambda = {latex_complex(evalue)}, \text{{ with }} " + 
                       fr"{latex_vector(v)} \text{{ and }} {latex_vector(w)}$")
            i += 2
        else:
            # Draw the eigenline:
            v = axismax * basis[i].normalized()
            fig.add_trace(plotly_line((-2*v, 2*v), color=color, 
                                      line_width=4, showlegend=False))
            plotly_add_text3d(fig, f"${other_vars[i]}$", 0.9*v, color=color, 
                              font_size=16, yshift=12)
            # Choose some initial states:
            if evalue > 0:
                axis_info[f"{other_vars[i]} axis"] = (color, [-0.1*v, 0.1*v])
                initial_states.append((0.1 * axismax, -0.1 * axismax))
            elif evalue < 0:
                axis_info[f"{other_vars[i]} axis"] = (color, [v, v])
                initial_states.append((0.6 * axismax, -0.6 * axismax))
            else:
                initial_states.append((0, 0.2 * axismax, 0.4 * axismax, 0.6 * axismax, 
                                       -0.2 * axismax, -0.4 * axismax, -0.6 * axismax))
            # Add the text for the eigenvalue and eigenvector:
            y -= 0.15
            text[y] = (fr"$\lambda = {float(round(evalue, 4)):g}, " + 
                       fr"\text{{ with }} {latex_vector(basis[i])}$")
            i += 1
    text[y - 0.06] = r"$\text{Diagonalization:}$"
    text[y - 0.24] = f"${latex_linear_diffeq(D, state_vars=('Q', 'R', 'S'))}$"
    for y, label in text.items():
        fig.add_annotation(plotly_text(label, (0, y), paper=True, font_size=12))

    for name, (color, states) in axis_info.items():
        for i, initial_state in enumerate(states):
            solution = desolve_odeint(vectorfield, initial_state, t_range, [X, Y, Z])
            fig.add_traces(plotly_trajectory3d(solution, color=color, 
                                               name=name + "X"*i, showlegend=not i, 
                                               visible="legendonly", legendgroup=name))
    T0 = column_matrix([v.normalized() for v in basis])
    initial_states = [T0*vector(state) for state in product(*initial_states)]
    names = [f"Trajectory {i + 1}" for i in range(len(initial_states))]
    for initial_state, name in zip(initial_states, names):
        solution = desolve_odeint(vectorfield, initial_state, t_range, [X, Y, Z])
        fig.add_traces(plotly_trajectory3d(solution, name=name, visible="legendonly"))
    points = np.array(initial_states, dtype=float)
    fig.add_scatter3d(x=points[:,0], y=points[:,1], z=points[:,2], 
                      mode="markers", marker_color="black", marker_size=2, 
                      text=names, hoverinfo="text", showlegend=False)
    return fig

Learning goal:

• Be able to describe the dynamics of any linear differential equation, in 3 dimensions and even in higher dimensions.

Example 1: Growth and decay

eigenvalues = [1, -1, -2]
v1, v2, v3 = matrix([(0, 1, 1), (1, 0, -1), (2, -2, 1)]).rows()
linear_diffeq_3d_interactive(eigenvalues, (v1, v2, v3), field=RDF)

Example 2: All decaying

eigenvalues = [-1, -2, -4]
v1, v2, v3 = matrix([(0, 1, 1), (1, 0, -1), (1, -1, 0)]).rows()
linear_diffeq_3d_interactive(eigenvalues, (v1, v2, v3), field=RDF)

Example 3: All growing

eigenvalues = [1/2, 1, 3]
v1, v2, v3 = matrix([(1, 2, 2), (2, 1, -2), (2, -2, 1)]).rows()
linear_diffeq_3d_interactive(eigenvalues, (v1, v2, v3))

Example 4: One decaying, two growing

eigenvalues = [1, 2, -1/2]
v1, v2, v3 = matrix([(0, 1, 1), (2, 1, -2), (1, -1, 0)]).rows()
linear_diffeq_3d_interactive(eigenvalues, (v1, v2, v3), field=RDF)

Example 5: Spiralling in and growing?

eigenvalues = [1, -2 + 15*I]
v1, v2, v3 = matrix([(0, 3, 1), (-1, 1, 0), (1, -1, -2)]).rows()
linear_diffeq_3d_interactive(eigenvalues, (v1, v2, v3))

Example 6: Spiralling out and decaying?

eigenvalues = [1 + 12*I, -5]
v1, v2, v3 = matrix([(0, -1, 1), (1, 1, 0), (-3, 1, 2)]).rows()
linear_diffeq_3d_interactive(eigenvalues, (v1, v2, v3))

Example 7: Spiralling in and decaying

eigenvalues = [-5, -1 + 18*I]
v1, v2, v3 = matrix([(3, 0, 1), (0, 1, 1), (0, 3, -1)]).rows()
linear_diffeq_3d_interactive(eigenvalues, (v1, v2, v3))

Conclusions:

• In 3 dimensions, linear models work essentially the same way that they do in 2 dimensions:

  • Each eigenvalue has a corresponding eigenvector line, which we treat as an axis in the $Q,R,S$ coordinate system.

  • The eigenvalue determines whether the corresponding variable ($Q$, $R$, or $S$) will exponentially grow or decay.

  • For a complex pair of eigenvalues, we can use the “almost-diagonalizing” trick. Then the corresponding pair of axes span a plane in the state space, and the behavior in that plane is either spiraling in or spiraling out, determined by the eigenvalue.

• All of this works the same way for differential equations (continuous time) and discrete-time models. The only difference is how you interpret exponential growth and decay.

  • For differential equations, $\text{Re}(\lambda) > 0$ means growth, $\text{Re}(\lambda) < 0$ means decay.

  • For discrete-time models, $\lvert \lambda \rvert > 1$ means growth, $\lvert \lambda \rvert < 1$ means decay.