Mini-Lecture 6-31 - Linear differential equations (continuous time) - The diagonal case.ipynb

Kernel: SageMath 9.0

from itertools import product

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

Learning goals:

• Be able to describe the behavior, both short-term and long-term, for any linear differential equation (continuous-time) model whose matrix is diagonal.

• For such a model, be able to identify the dominant eigenvalue, and be able to explain what this means about the long-term behavior of the model.

Linear differential equation models, with a diagonal matrix:

$$\begin{bmatrix} X' \\ Y' \end{bmatrix} = \begin{pmatrix} \lambda_1 & 0 \\ 0 & \lambda_2 \end{pmatrix} \begin{bmatrix} X \\ Y \end{bmatrix}$$

or

$$\begin{cases} X' = \lambda_1 X \\ Y' = \lambda_2 Y \end{cases}$$

Solution:

$$X(t) = X(0) e^{\lambda_1 t} \qquad \text{and} \qquad Y(t) = Y(0) e^{\lambda_2 t}$$

def colorchoice(a):
    return "green" if a < 0 else "red" if a > 0 else "purple"

def diagonal_interactive(eigenvalues, initial_state=None):
    axismax = 10
    t_range = srange(0, 15.01, 0.05)
    xcolor = "red"
    ycolor = "green"
    def initchoice(a):
        if a < 0:
            return round((0.8 + uniform(-0.2, 0.2))*axismax, 1)
        if a > 0:
            return round((0.2 + uniform(-0.2, 0.2))*axismax, 1)
        return round(uniform(0, axismax), 1)
    if initial_state is None:
        initial_state = vector([initchoice(eigenvalues[0]), initchoice(eigenvalues[1])])
    else:
        initial_state = vector(initial_state)
    D = diagonal_matrix(RDF, eigenvalues)
    vectorfield(X, Y) = tuple(D*vector((X, Y)))
    solution = desolve_odeint(vectorfield, initial_state, t_range, [X,Y])

    axisfield = Graphics()
    wholefield = Graphics()
    xylist = srange(-0.95*axismax, axismax, 0.1*axismax)
    scale = 0.1*axismax/max([vectorfield(x,0).norm() for x in xylist] + 
                            [vectorfield(0,y).norm() for y in xylist])
    for x in xylist:
        axisfield += arrow2d((x,0), vector((x,0)) + scale*vectorfield(x,0), 
                             color=xcolor, width=1, arrowsize=2, zorder=4)
    for y in xylist:
        axisfield += arrow2d((0,y), vector((0,y)) + scale*vectorfield(0,y), 
                             color=ycolor, width=1, arrowsize=2, zorder=4)
    for x, y in product(xylist, xylist):
        wholefield += arrow2d((x,y), vector((x,y)) + scale*vectorfield(x,y), 
                              color="darkgray", width=1, arrowsize=2)

    dont_update = False
    def update(t, showspikes, showaxisfield, showwholefield, axismax):
        if dont_update:
            return
        p = point(solution[t], size=30, color="blue")
        p += list_plot(solution[:t+1], plotjoined=True, color="blue")
        if showspikes:
            x0, y0 = solution[t]
            p += line(((x0,0), (x0,y0), (0,y0)), color="lightgray")
            p += point((x0,0), size=30, color=xcolor)
            p += point((0,y0), size=30, color=ycolor)
            p += plot(vectorfield(x0, y0), start=(x0,y0), color="darkgray", arrowsize=2)
            p += plot(vectorfield(x0, 0), start=(x0,0), color=xcolor, thickness=3, arrowsize=2)
            p += plot(vectorfield(0, y0), start=(0,y0), color=ycolor, thickness=3, arrowsize=2)
        if showaxisfield:
            p += axisfield
        if showwholefield:
            p += wholefield
        p.show(xmin=-axismax, xmax=axismax, ymin=-axismax, ymax=axismax, 
               axes_labels=("$X$", "$Y$"), aspect_ratio=1, figsize=9.5)

    def change_initial_state(change):
        nonlocal dont_update, solution
        initial_state = change["new"].strip("()[]").split(",")
        try:
            initial_state = vector([float(x) for x in initial_state])
            solution = desolve_odeint(vectorfield, initial_state, t_range, [X,Y])
        except:
            return
        dont_update = True
        t.value = 1
        showspikes = False
        showaxisfield = False
        showwholefield = False
        dont_update = False
        t.value = 0

    layout = Layout(width="250px", margin="15px 0px 15px 0px")
    label = (r"$$\begin{bmatrix} X' \\ Y' \end{bmatrix} = %s "
             r"\begin{bmatrix} X \\ Y \end{bmatrix}$$")
    equation = HTMLMath(label % latex(D), layout=layout)
    t_values = {round(t,2): n for n, t in enumerate(t_range)}
    t  = SelectionSlider(options=t_values, layout=layout, description="$t$ for solution:")
    showspikes = Checkbox(value=False, layout=layout, description="Show XY lines")
    showaxisfield = Checkbox(value=False, layout=layout, description="Show vec fld on axes")
    showwholefield = Checkbox(value=False, layout=layout, description="Show whole vec fld")
    zoom = SelectionSlider(options=[1, 2, 5, 10, 20, 50, 100, 200, 500, 1000], value=10, 
                           layout=layout, description="Zoom:")
    state = Text(value=str(initial_state), continuous_update=False, 
                 layout=layout, description="Initial state:")
    state.observe(change_initial_state, "value")
    output = interactive_output(update, dict(t=t, showspikes=showspikes, showaxisfield=showaxisfield, 
                                             showwholefield=showwholefield, axismax=zoom))
    controls = VBox((equation, state, t, showspikes, showaxisfield, showwholefield, zoom))
    display(HBox((controls, output)))

diagonal_interactive((-0.2, 0.5))

diagonal_interactive((0.4, -0.6))

$X$ eigenvalue	$Y$ eigenvalue	Eq. pt. at $(0,0)$	Long term behavior
-0.2	0.5	Saddle point	Exp. growth along $Y$ axis
0.4	-0.6	Saddle point	Exp. growth along $X$ axis

diagonal_interactive((-0.3, -0.6))

diagonal_interactive((-0.25, -0.1))

$X$ eigenvalue	$Y$ eigenvalue	Eq. pt. at $(0,0)$	Long term behavior
-0.2	0.5	Saddle point	Exp. growth along $Y$ axis
0.4	-0.6	Saddle point	Exp. growth along $X$ axis
-0.3	-0.6	Sink	Exp. decay along $X$ axis
-0.25	-0.1	Sink	Exp. decay along $Y$ axis

diagonal_interactive((0.5, 1.5))

diagonal_interactive((1.1, 0.8))

diagonal_interactive((-1.2, 0.3))

$X$ eigenvalue	$Y$ eigenvalue	Eq. pt. at $(0,0)$	Long term behavior
-0.2	0.5	Saddle point	Exp. growth along $Y$ axis
0.4	-0.6	Saddle point	Exp. growth along $X$ axis
-0.3	-0.6	Sink	Exp. decay along $X$ axis
-0.25	-0.1	Sink	Exp. decay along $Y$ axis
0.5	1.5	Source	Exp. growth parallel to $Y$ axis
1.1	0.8	Source	Exp. growth parallel to $X$ axis
-1.2	0.3	Saddle	Exp. growth along $Y$ axis

Conclusions:

If multiple (two or more) variables in a differential equation model are exponentially growing/decaying independently of each other (i.e. the model is linear with a diagonal matrix), then the resulting behavior is as follows:

1. Overall behavior:

   • If all eigenvalues are positive, then any starting state will grow away from the origin. That is, the equilibrium point at the origin is a source (unstable).

   • If all eigenvalues are negative, then any starting state will grow toward the origin. That is, the equilibrium point at the origin is a sink (stable).

   • If some eigenvalues are negative and some others are positive, then a starting state will move inward along some axes, but grow outward along others. That is, the equilibrium point at the origin is a saddle point (unstable).

2. The long-term behavior will be in the direction of whichever axis has the dominant eigenvalue. In this setting, the dominant eigenvalue is the one that is greatest, as a real number. (It's not the one whose absolute value is greatest)