Small_Oscilations.ipynb

Kernel: SageMath 9.2

Small oscillations

In this notebook we consider a model of a two dimensional linear physical system which conserves energy. There are very many physical systems which can be described by this mathematics. Our linear $2$-dimensional system is governed by Newton's law $F = ma$. We assume that the force doesn't depend on the velocity (for example, a particle in some potential). In our notation, this means $\mathbf{x}'' = \frac{1}{m} F(\mathbf{x}).$ Our assumption that the system is linear means that $F(\mathbf{x}) = A\mathbf{x}$ for some matrix $A$. We can choose units where $m = 1$. Thus our system is $$\mathbf{x}'' = A \mathbf{x}.$$

Just as we convert an $n$-th order one-dimensional ODE to a $n$-dimensional system of ODEs, we can convert this 2nd order 2-dimensional system into a 1st order 4-dimensional system as follows: let $$\mathbf{x} = \begin{pmatrix} x_1\\x_2\end{pmatrix}.$$ Then define $p_1 = x_1'$ and $p_2 = x_2'$. It follows that $x_1'' = p_1'$ and $x_2'' = p_2'$. In vector notation, we have $\mathbf{p} = \mathbf{x}'$, so that $\mathbf{p}' = \mathbf{x}''$. Our system can then be written in the form of a 1st order $4\times 4$ linear system, which is given by a block matrix as follows:

$$\begin{pmatrix} \mathbf{p} \\ \mathbf{x} \end{pmatrix}' = \left( \begin{array}{c|c} 0 & A \\ \hline I & 0 \end{array}\right) \begin{pmatrix} \mathbf{p} \\ \mathbf{x} \end{pmatrix}$$

where our matrix $B$ is a $4\times 4$ matrix made up of $2\times 2$ blocks.

Let's use our last assumption in the problem, that the system conserves energy. Mathematically the corresponds to the fact that the eigenvalues $r + i \omega$ of $B$ cannot have non-zero real part, since if this were the case there would be $e^{rt}$ factors in the solutions, which either contract to $0$ or escape to $\infty$, with exponentially increasing or decreasing velocity. It follows that the eigenvalues of $B$ are purely imaginary. Since $B$ is a real matrix, its eigenvalues must be of the form $$i\omega_1, -i\omega_1, i\omega_2, -i\omega_2$$ so that our general solution (let us ignore the issue of Jordan blocks...) is $$\begin{equation}\begin{pmatrix} \mathbf{p} \\ \mathbf{x} \end{pmatrix} = c_1 \cos(\omega_1 t) \xi_1 + c_2 \sin(\omega_1 t) \xi_2 + c_3 \cos(\omega_2 t) \xi_3 + c_4 \sin(\omega_2 t) \xi_4 \end{equation}$$ for vectors $\xi_1, \ldots, \xi_4\in \mathbb{R}^4$.

It's clear that we have two decoupled oscillatory systems, one in the $\xi_1$ and $\xi_2$ directions with frequency $\omega_1$, and another in the $\xi_3$ and $\xi_4$ directions with frequency $\omega_2$. In hindsight we can say more: our original matrix $A$ must have had eigenvalues $-\omega_1^2$ and $-\omega_2^2$ because of the possible choices for the second derivative of our solution (1). Thus our original system decouples into a pair of oscillators $x_1'' = -\omega_1^2 x_1$ and $x_2'' = -\omega_2^2 x_2$, after a linear change of coordinates.

Since we cannot easily visualize the 4-dimensional phase space of the system, let's only plot the $2$-dimensional coordinates $\mathbf{x}$. All of this analysis shows that, after a linear change of $x$ coordinates to diagonalize $A$, we know the solutions are of the form $$\begin{align*} x_1(t) &= c_1 \sin(\omega_1 t + \delta)\\ x_2(t) &= c_2 \sin(\omega_2 t + \epsilon). \end{align*}$$ Shifting the time coordinate $t \mapsto t - \epsilon / \omega_2$ we can assume that $\epsilon = 0$, while subsequently rescaling $t\mapsto t/\omega_1$, we can assume $\omega_1 = 1$. So there are really only 2 free parameters, $\omega_2$ and $\delta$. We will now plot these solution curves. A few remarks:

1. These curves are known as Lissajous curves.

2. Note that in the 4-dimensional system the curves do not intersect (this would violate the uniqueness half of the Existence/Uniqueness theorem), but when we project to two dimensions, we get intersecting solution curves.

3. Along the same lines, as we vary $\delta$ we can alternatively interpret this literally as a rotation in the $4$-dimensional phases space, or a $3$-dimensional subspace thereof. It is a good exercise to make this precise mathematically. It explains the phenomenon we notice as we vary $\delta$.

4. Notice also the following: if $\omega_2/\omega_1$ is a rational number $a/b \in \mathbb{Q}$ then the curve eventually returns to itself and repeats. On the other hand, if $\omega_2/\omega_1 \not \in \mathbb{Q}$ then the curve will eventually fill the entire rectangle. Why is this?

@interact
def lissajous(omega_2 = input_box(default = 1.5, label = "$\\omega_2$"), 
             delta = slider(vmax = 6.28, vmin = 0, step_size = 0.01,default = pi/4, label = "$\\delta$"),
             tmax = slider(vmax = 500, vmin = 10, step_size = 5, default = 20, label = "max $t$"),
             tolerance = slider(vmax = 0.01, vmin = 0.00001, step_size = 0.0001, default = 0.01, label = "Plot tolerance"),
             show_axes = checkbox(default = True, label = "Show axes?"),
             t_values = input_box(default = [2,4,6,8,10], label = "$t$ values"),
             show_t = checkbox(default = False, label = "Plot $t$ values?"),
             ):
    t = var('t')
    if show_t:
        maxt = max(t_values)
        points = sum(point((sin(t0 + delta), sin(omega_2*t0)), size = 22, color = Color(1,1.0-t0/maxt,0)) for t0 in t_values)
    else: 

        points = point((0,0), size = 0)
    show(points + parametric_plot((sin(t + delta), sin(omega_2*t)), (t,0,tmax),adaptive_tolerance = tolerance, thickness = .2, figsize = 8, axes  = show_axes))

We can see that it appears to rotate in a higher-dimensional space as we vary $\delta$, despite our problem initially appearing 2-dimensional. These patterns show up on oscilloscopes (this rotation here comes from slightly irrational $\omega_1/\omega_2$, which end up being equivalent to varying $\delta$.) Regardless we can actually plot a third component of our $4$-dimensional system and see that we really do have a rotating figure in a higher dimensional space. Our 2d picture we started with doesn't tell the whole story.

The $2d$ picture we have plots the $x_1$ and $x_2$ components of the vectors of our solution. If we plotted the $p_1$ component we would get the 3d rotation but our curve will still intersect itself. So we choose a rather strange $3d$ slice of the 4d system where we plot the coordinates $$x_1, x_2, p_1 + \epsilon p_2$$ just to demonstrate that the solution curve actually doesn't intersect itself, and the rotation we see in the 2d picture actually has a physical interpretation in terms of a rotation in the coordinates $x_1$ and $p_1$.

@interact
def lissajous3d(omega_2 = input_box(default = 1.5, label = "$\\omega_2$"), 
             delta = slider(vmax = 6.28, vmin = 0, step_size = 0.01,default = pi/4, label = "$\\delta$"),
             tmax = slider(vmax = 500, vmin = 10, step_size = 5, default = 20, label = "max $t$"),
             show_axes = checkbox(default = False, label = "Show axes?"),
             t_values = input_box(default = [2,4,6,8,10], label = "$t$ values"),
             epsilon = slider(vmax = 0.2,vmin = 0, default = 0.05, step_size = 0.02, label = '$\\epsilon$'),
             plot_points = slider(vmax = 2000, vmin = 150, default = 200, step_size = 50, label = "plot points"),
             show_t = checkbox(default = False, label = "Plot $t$ values?"),
             ):
    t = var('t')
    if show_t:
        maxt = max(t_values)
        points = sum(point((sin(t0 + delta), sin(omega_2*t0), cos(t0 + delta) + epsilon * cos(omega_2*t0)), size = 8, color = Color(1,1.0-t0/maxt,0)) for t0 in t_values)
    else: 

        points = point((0,0), size = 0)
    show(points + parametric_plot3d((sin(t + delta), sin(omega_2*t), cos(t + delta)  + epsilon * cos(omega_2*t)), (t,0,tmax),plot_points = plot_points, thickness = .2, figsize = 8, axes  = show_axes))

Since it is hard for sage to update quickly enough to animate this and see rotation as we vary $\delta$, here is an animation of the 3d rotation for a few values of $\omega_2/\omega_1$ and the corresponding 2d Lissajous figures.

def animated_3d_lissajous(omega_2, maxt,num_t_values = 5, epsilon=0.01, num_frames = 30, frame =True):
    curve(t, delta) = (sin(t + delta), sin(omega_2*t), cos(t + delta) + epsilon *cos(omega_2*t ))
    curve2d(t, delta) = (sin(t + delta), sin(omega_2*t))
    t_values = list(srange(0,maxt - 0.01, maxt/num_t_values))
    plt1 = [sum(point(curve(t0,delta), size = 8, color = Color(1,1.0-t0/maxt,0)) for t0 in t_values) 
                + parametric_plot3d(curve(t,delta), (t, 0,maxt), frame = frame,figsize = 8)
                 for delta in srange(0, 2*pi, 2*pi/num_frames)]
    a = animate(plt1)
    plt2 = [sum(points(curve2d(t0, delta), size = 22, color = Color(1,1.0-t0/maxt,0)) for t0 in t_values) + 
                parametric_plot(curve2d(t,delta), (t,0,maxt), figsize= 8, thickness = 0.3, axes = False)
                for delta in srange(0,2*pi, 2*pi/num_frames)]
    b = animate(plt2)
    return((a,b))

def save_animated_3d_lissajous(omega_2, maxt,num_t_values = 5, epsilon=0.01, num_frames = 45, frame =True):
    curve(t, delta) = (sin(t + delta), sin(omega_2*t), cos(t + delta) + epsilon *cos(omega_2*t ))
    curve2d(t, delta) = (sin(t + delta), sin(omega_2*t))
    t_values = list(srange(0,maxt - 0.01, maxt/num_t_values))
    plt1 = [sum(point(curve(t0,delta), size = 8, color = Color(1,1.0-t0/maxt,0)) for t0 in t_values) 
                + parametric_plot3d(curve(t,delta), (t, 0,maxt), frame = frame,figsize = 8)
                 for delta in srange(0, 2*pi, 2*pi/num_frames)]
    a = animate(plt1)
    plt2 = [sum(points(curve2d(t0, delta), size = 22, color = Color(1,1.0-t0/maxt,0)) for t0 in t_values) + 
                parametric_plot(curve2d(t,delta), (t,0,maxt), figsize= 8, thickness = 0.3, axes = False)
                for delta in srange(0,2*pi, 2*pi/num_frames)]
    b = animate(plt2)
    a.apng(savefile = "./ratio" + str(RDF(omega_2))[:4] + "_3d_lissajous.apng")
    b.apng(savefile = "./ratio" + str(RDF(omega_2))[:4] + "_2d_lissajous.apng")
           
omega2s = [1,2, 3/2,4/3, 5/3, 5/4, 6/5]
tmaxs = [10,20,20,20,30,35, 35]
for (o2, tmx) in zip(omega2s, tmaxs):
    save_animated_3d_lissajous(o2, tmx)

Ratio $\omega_2 : \omega_1$	3d Figure	2d Figure
$1:1$
$2:1$
$3:2$
$4:3$
$5:3$
$5:4$
$6:5$