import numpy
from matplotlib import pyplot

%matplotlib inline

Going further in understanding coupled oscillators¶

Last time, we investigated a system consisting of two masses $m$ coupled to each other by a spring of spring constant $k$. These masses were also coupled to walls with springs having the same spring constant $k$, so the whole system had the following configuration:

WALL] /\/\/\/\/ [BLOCK 1] /\/\/\/\/ [BLOCK 2] /\/\/\/\/ [WALL

Assume, for simplicity, that the springs are at their natural lengths when the system is sitting at rest. You were asked to use Euler's method to determine what the motion of this system would look like given particular initial conditions, and then you were asked to investigate some properties of the motion under general initial conditions. In this session, we will extend our analysis of coupled oscillators in the following ways:

We will

Learn a new alogirithm for computing the numerical solutions to differential equations that is more stable than Euler's Method
Test this new algorithm on the system consisting of a single harmonic oscillator to see how it behaves differently
Use this new algorithm to reproduce the behavior that we saw with a pair of coupled oscillators
Gain a deeper understanding of what is happening with the pair of coupled oscillators by analyzing that system using linear algebra
Extend our linear algebraic analysis of the pair of coupled oscillators to predict the behavior of a system in which there are three or more coupled oscillators
Check our predictions by numerically solving the differential equations for three or more oscillators

A new numerical ODE solving algorithm¶

Recall that Euler's method works as follows:

Start with the equation of motion of your system which will be of the form

$$ \ddot x = \frac{F}{m}. $$

It's a second-order ODE for the unknown function $x(t)$. If the system consists of an object moving in more than one dimension, then we would write the position and force as vectors. If there is more than one object in the system, then there will be one such equation of motion for each.

The next step is to introduce the velocity $v(t)$ as a new unknown function and convert your original equation of motion into a pair first-order differential equations of the form:

\begin{align} \dot x(t) &= v(t) \\ \dot v(t) &= \frac{F}{m} \end{align}

Then, choose some initial conditions $x_0$ and $v_0$, choose a time step $\Delta t$, and choose a number of steps. Generate a discrete approximation to $x(t)$ and $v(t)$ in the form of a sequence of $x$ values and $v$ values by using the following iteration scheme:

\begin{align} x_{n+1} = x_n + v_n\, \Delta t \\ v_{n+1} = v_n + \frac{F_n}{m}\, \Delta t \\ \end{align}

Amazingly, by altering this scheme only slightly, we will able to vastly improve its performance on the sort of oscillating systems we've been working with. The new scheme is called the Euler-Cromer algorithm, and it works like this:

\begin{align} v_{n+1} &= v_n + \frac{F_n}{m}\, \Delta \\ x_{n+1} &= x_n + v_{n+1}\,\Delta t \end{align}

Exercise.¶

What are the differences between the Euler-Cromer scheme and the Euler scheme?

Answer.¶

Euler-Cromer uses the result of $v_{n+1}$ to calculate the next $x_{n+1}$ also velocity does not depend on position.

Exercise.¶

In a code cell below, make two copies of your code from two sessions ago in which you solved for and plotted the position versus time of a single harmonic oscillator using the Euler scheme.
Alter one copy of the code so that it runs the Euler-Cromer scheme instead.
Run both versions of the code to generate plots of the position as a function of time for a harmonic oscillator, and plot these positions versus time on the same plot. Make sure to use enough time steps that you see many periods.
What is the difference in the plots between the Euler scheme and the Euler-Cromer scheme?

### Answer.

m = 2.0
k = 1.0
x_0 = 1.0
v_0 = 0.0
delta_t = 0.5
num_steps = 50

t_array = numpy.zeros(num_steps + 1)
x_array = numpy.zeros(num_steps + 1)
v_array = numpy.zeros(num_steps + 1)
v_ecarray = numpy.zeros(num_steps + 1)
x_ecarray = numpy.zeros(num_steps + 1)

t_array[0] = 0.0
x_array[0] = 1.0
v_array[0] = 0.0
v_ecarray[0] = 0.0
x_ecarray[0] = 1.0
for n in range(num_steps):
    
    t_array[n + 1] = t_array[n] + delta_t
    x_array[n + 1] = x_array[n] + v_array[n] * delta_t
    v_array[n + 1] = v_array[n] - (k / m) * x_array[n] * delta_t
    
#EC

    v_ecarray[n + 1] = v_ecarray[n] + (-1 * k * x_ecarray[n] / m)  * delta_t
    x_ecarray[n + 1] = x_ecarray[n] + v_ecarray[n+1] * delta_t
    
pyplot.plot(t_array, x_array, 'green') #plot of Euler's method
pyplot.plot(t_array, x_ecarray, 'red') #plot of Euler-Cromer's method

# Euler-Cromer's graph is more concave. Euler's method produces more error

[<matplotlib.lines.Line2D at 0x7f752c443668>]

Understanding coupled oscillators using linear algebra¶

As we'll see in a moment, we can use linear algebraic concepts to better predict and understand what the system of coupled oscillators will do. Doing linear algebra with python is pretty easy using numpy. For example, to specify a 2-by-2 matrix $A$, we would type

A = numpy.array([[1, 0], [0, 2]])

This defines the 2-by-2 matrix

$$ A = \begin{pmatrix} 1 & 0 \\ 0 & 2 \end{pmatrix} $$

What the code is doing is taking a list of lists (aka a nested list) [[1, 0], [0, 1]] and converting it into an object called a numpy array which can be manipulated like a matrix. For example, we can find its eigenvalues and corresponding eigenvectors as follows:

from numpy import linalg #imports a linear algebra package in numpy (once you import once, you don't have to do it again)

linalg.eig(A)

(array([ 1.,  2.]), array([[ 1.,  0.],
        [ 0.,  1.]]))

As you can see, this command outputs one array the gives the eigenvalue $1$ and $2$, and another array that gives the corresponding eigenvectors $\binom{1}{0}$ and $\binom{0}{1}$.

Analyzing coupled oscillators using linear algebra¶

We will now explore how eigenvalues and eigenvectors can give us insight into the dynamics of a system of couple oscillators through a sequence of exploratory exercises:

Write the equation of motion for the system of two coupled oscillators in matrix form as follows:

\begin{align} \begin{pmatrix} \ddot x_1 \\ \ddot x_2 \end{pmatrix} = -A\begin{pmatrix} x_1 \\ x_2 \end{pmatrix} \end{align}

where $A$ is a 2-by-2 matrix that depends on the $\omega = \sqrt{k/m}$.

Choose $\omega = 1$, and use python to numerically determine the eigenvalues and corresponding eigenvectors of the matrix $A$.
The first eigenvalue $\lambda$ will have a corresponding eigenvector $\binom{a}{b}$. Run the Euler-cromer scheme for the pair of coupled oscillators with initial data $x_{1,0} = a$, $x_{2,0} = b$, $v_{1,0} = 0$, and $v_{2,0} = 0$. Plot the time versus position of the masses compare the frequency of the oscillations of the two masses with the square root of the eigenvalue $\lambda$ -- see anything interesting?
Repeat for the second eigenvalue and its corresponding eigenvector.
The mode of oscillation corresponding to each eigenvector is called a normal mode. It turns out that every motion of the system of coupled oscillators can be built out of these normal modes of oscillation, although we won't see why this is the case quite yet.

Consider now a system of three coupled oscillators:

  WALL] /\/\/\/\/ [BLOCK 1] /\/\/\/\/ [BLOCK 2] /\/\/\/\/ [BLOCK 3] /\/\/\/\/ [WALL

Write the equation of motion for a system of three coupled oscillators in matrix form as follows:

\begin{align} \begin{pmatrix} \ddot x_1 \\ \ddot x_2 \\ \ddot x_3 \end{pmatrix} = -A\begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} \end{align}

Choose $\omega = 1$, and us python to numerically determine the eigenvalues and corresponding eigenvectors of the matrix $A$.
The eigenvalue $\lambda$ will have a corresponding eigenvector $(a, b, c)$. Run the Euler-cromer scheme for the triple of coupled oscillators with initial data $x_{1,0} = a$, $x_{2,0} = b$, $x_{3,0} = c$, $v_{1,0} = 0$, $v_{2,0} = 0$, $v_{3,0} = 0$. Plot the time versus position of the masses and compare the frequency of the oscillations of the masses with the square root of the eigenvalue $\lambda$ -- see anything interesting?
Repeat for the second eigenvalue and its corresponding eigenvector.
Repeat for the third eigenvalue and its corresponding eigenvector.
Each of these eigenvalues and its corresponding eigenvector is a normal mode of this system of three coupled oscillators, and just as in the case of two oscillators, an arbitrary motion of the system as a whole can be built out of these normal modes. This analysis can actually be extended to any number of oscillators one desires. What would the matrix $A$ look like in general?

\begin{align} \begin{pmatrix} \ddot x_1 \\ \ddot x_2 \end{pmatrix} = -\begin{pmatrix} 2k/m & -k/m\\ -k/m & 2k/m \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} \end{align}

k = 1.0
m = 1.0
w = (k / m)**(1/2)


A = numpy.array([[2*w**2, -1*w**2], [-1*w**2, 2*w**2]])

linalg.eig(A)

(array([ 3.,  1.]), array([[ 0.70710678,  0.70710678],
        [-0.70710678,  0.70710678]]))

print(A[1,1])

2.0

m = 1.0
k = 0.5
x_1 = [0.707]
v_1 = [0.0]
x_2 = [0.707]
v_2 = [0.0]
delta_t = 0.5
num_steps = 50
j = 1

t_array = numpy.zeros(num_steps + 1)
x_1array = numpy.zeros(num_steps + 1)
x_2array = numpy.zeros(num_steps + 1)
v_1array = numpy.zeros(num_steps + 1)
v_2array = numpy.zeros(num_steps + 1)

for j in range(1):
    t_array[0] = 0.0
    x_1array[0] = x_1[j]
    v_1array[0] = v_1[j]
    x_2array[0] = x_2[j]
    v_2array[0] = v_2[j]
    n = 0
    for n in range(num_steps):
        t_array[n + 1] = t_array[n] + delta_t
        v_1array[n + 1] = v_1array[n] + (( 2*k*x_1array[n] - k*x_2array[n] ) / m)  * delta_t
        v_2array[n + 1] = v_2array[n] + ( ( 2*k*x_2array[n] - k*x_1array[n] ) / m)  * delta_t
        x_1array[n + 1] = x_1array[n] + v_1array[n+1] * delta_t
        x_2array[n + 1] = x_2array[n] + v_2array[n+1] * delta_t
        
pyplot.plot(x_1array, t_array, 'red')

[<matplotlib.lines.Line2D at 0x7f752c1f9e48>]