Coupled Harmonic Oscillator

Teng-Jui Lin

Content adapted from Mr. Ben Trey’s coupled harmonic oscillator internet. The application of scipy.integrate.solve_ivp() and phase portrait is related to UW AMATH 301, Beginning Scientific Computing.

• Phase portraits

  • Coupled harmonic oscillator

  • Coupled harmonic oscillator internet

Coupled harmonic oscillator

Problem Statement. Consider a coupled harmonic oscillator with two blocks $m_1$ and $m_2$ attached with each other and two fixed walls with springs with spring constants $k_1,k_2,k_3$. It is supported by a smooth horizontal surface. The equilibrium length of the springs are $l_1,l_2,l_3$, and the distance between the two walls is $L$. The damping coefficient of each block is $c_1$ and $c_2$. The position of the blocks is $x_1$ and $x_2$ with origin at the equilibrium position of $m_1$ and $m_2$, respectively; the velocity of the blocks is ${\dot{x}}_1$ and ${\dot{x}}_2$; the acceleration of the blocks is ${\ddot{x}}_1$ and ${\ddot{x}}_2$. External forces on the blocks are $F_{\text{ext},1}$ and $F_{\text{ext},2}$. Define positive axis to the right.

Suppose $m_1$ is initially pulled away from its equilibrium position and released without external forces. Solve the position of the blocks $x_1$ and $x_2$ over time with the initial conditions $x_1(0)=-0.2\mathrm{m},x_2(0)=0.2\mathrm{m},v_1(0)=0,v_2(0)=0$, provided the following system specifications:

• Masses: $m_1=m_2=1\mathrm{kg}$

• Spring constants: $k_1=k_2=k_3=100\mathrm{N}/\mathrm{m}$

• Equilibrium length of springs: $l_1=l_2=l_3=1\mathrm{m}$

• Damping constants: $c_1=c_2=3\mathrm{kg}/\mathrm{s}$

Deliverables:

• Generate a plot of $x_1(t)$ and $x_2(t)$ over time

• Generate an animation of the system

Solution. From the above specifications, we can get some additional information. Since we have no external forces, so we have $F_{\text{ext},1}=F_{\text{ext},2}=0$. We assume the blocks are point masses, so the distance between the walls is $L=l_1+l_2+l_3=3\mathrm{m}$.

Now, we analyze the forces on each block. The net vertical force is zero because gravitational force and the normal force cancels. Friction can be neglected because of smooth surface. We then look at the horizontal forces for $m_1$ and $m_2$.

For $m_1$, we have Newton’s second law and the horizontal forces

$$\begin{array}{r}\begin{array}{rlrl}\text{Newton's second law:}& & \sum F_1& =m_1{\ddot{x}}_1\\ \text{Spring 1:}& & F_{\text{spring 1}}& =-k_1{x}_1\\ \text{Spring 2:}& & F_{\text{spring 2}}& =-k_2(x_1-x_2)\\ \text{Damping:}& & F_{\text{damping}}& =-c_1{\dot{x}}_1\\ \text{External:}& & F_{\text{ext,1}}& =0\end{array}\end{array}$$

which form the ODE

$$\begin{array}{r}\begin{array}{rl}\sum F_1& =F_{\text{spring 1}}+F_{\text{spring 2}}+F_{\text{damping}}+F_{\text{ext,1}}\\ m_1{\ddot{x}}_1& =-k_1{x}_1-k_2(x_1-x_2)-c_1{\dot{x}}_1\end{array}\end{array}$$

For $m_2$, we have Newton’s second law and the horizontal forces

$$\begin{array}{r}\begin{array}{rlrl}\text{Newton's second law:}& & \sum F_2& =m_2{\ddot{x}}_2\\ \text{Spring 2:}& & F_{\text{spring 2}}& =-k_2(x_2-x_1)\\ \text{Spring 3:}& & F_{\text{spring 3}}& =-k_3{x}_2\\ \text{Damping:}& & F_{\text{damping}}& =-c_2{\dot{x}}_2\\ \text{External:}& & F_{\text{ext,2}}& =0\end{array}\end{array}$$

which form the ODE

$$\begin{array}{r}\begin{array}{rl}\sum F_2& =F_{\text{spring 2}}+F_{\text{spring 3}}+F_{\text{damping}}+F_{\text{ext,2}}\\ m_2{\ddot{x}}_2& =-k_2(x_2-x_1)-k_3{x}_2-c_2{\dot{x}}_2\end{array}\end{array}$$

The above two ODEs are coupled system of linear second-order ODEs. They can be written into a system of

$$\begin{array}{r}\{\begin{array}{l}{\dot{x}}_1=v_1\\ {\dot{v}}_1={\displaystyle \frac{-k_1{x}_1-k_2(x_1-x_2)-c_1{\dot{x}}_1}{m_1}}\\ {\dot{x}}_2=v_2\\ {\dot{v}}_2={\displaystyle \frac{-k_2(x_2-x_1)-k_3{x}_2-c_2{\dot{x}}_2}{m_2}}\end{array}\end{array}$$

that can be solved by scipy.integrate.solve_ivp().

Solve the system

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import scipy
from scipy import integrate

# model params
m1 = 1  # kg
m2 = 1  # kg
k1 = 100  # N/m
k2 = 100  # N/m
k3 = 100  # N/m
c1 = 3  # kg/s
c2 = 3  # kg/s
l1 = 1  # m
l2 = 1  # m
l3 = 1  # m

# time array
t_initial = 0
t_final = 2
t = np.linspace(t_initial, t_final, 1000)
t_len = len(t)

# initial values (x is relative to eqm length)
x1_init = -0.2
v1_init = 0
x2_init = 0.2
v2_init = 0
initial_values = [x1_init, v1_init, x2_init, v2_init]

# external forces
F1 = lambda t : 0
F2 = lambda t : 0

# def F1(t):
#     '''External force on m1.'''
#     omega = 100
#     A = np.sqrt(k1 + k2)+1000
#     phi = 1
#     return A*np.sin(omega*t + phi)

# def F1(t):
#     if(int(t/3)%2 == 0):
#         omega = 1e4
#     else:
#         omega = np.sqrt(k1 + k2)
#     A = 20
#     phi = 1
#     return A*np.sin(omega*t + phi)

## equivalent, simpler form for coupled ode systems
## note: this form not convenient for phase portrait
def coupled_ode_syst(t, Y):
    # t -> independent variable
    # Y -> functions evaluated at independent variable
    x1, v1, x2, v2 = Y
    derivatives = np.array([
        v1,
        (-k1*x1 - k2*(x1 - x2) - c1*v1 + F1(t))/m1,
        v2,
        (-k2*(x2 - x1) - k3*x2 - c2*v2 + F2(t))/m2,
    ])
    return derivatives

ode_soln = scipy.integrate.solve_ivp(coupled_ode_syst, [t_initial, t_final], initial_values, t_eval=t).y
x1, v1, x2, v2 = ode_soln

Position plot

def custom_plot_settings():
    %config InlineBackend.figure_format = 'retina'
    %matplotlib inline

    plt.rcParams.update({
        'font.family': 'Arial',  # Times New Roman, Calibri
        'font.weight': 'normal',
        'mathtext.fontset': 'cm',
        'font.size': 18,

        'lines.linewidth': 2,

        'axes.linewidth': 2,
        'axes.spines.top': False,
        'axes.spines.right': False,
        'axes.titleweight': 'bold',
        'axes.titlesize': 18,
        'axes.labelweight': 'bold',

        'xtick.major.size': 8,
        'xtick.major.width': 2,
        'ytick.major.size': 8,
        'ytick.major.width': 2,

        'figure.dpi': 80,
        'savefig.dpi': 300,

        'legend.framealpha': 1, 
        'legend.edgecolor': 'black',
        'legend.fancybox': False,
        'legend.fontsize': 14,

        'animation.html': 'html5',
    })

custom_plot_settings()

fig, axs = plt.subplots(2, 1, figsize=(5, 6), sharex=True)
axs[0].plot(t, x1, label='$x_1(t)$')
axs[0].plot(t, x2, label='$x_2(t)$')
axs[0].set_ylabel('$x(t)$')
# axs[0].set_ylim(0.5, 2.5)

axs[1].plot(t, v1, label='$v_1(t)$')
axs[1].plot(t, v2, label='$v_2(t)$')
axs[1].set_xlabel('$t$')
axs[1].set_ylabel('$v(t)$')
# axs[1].set_ylim(-3.5, 3.5)

for i in range(2):
    axs[i].set_xlim(t_initial, t_final)
    axs[i].plot([t_initial, t_final], [0, 0], '--', color='grey', lw=0.5, zorder=0)  # zero ref

Animation of the system

Note on animation: for local reproducible results, download ffmpeg and add to path variable. For reproducible results online, use Google Colab and run the command below.

# Run the command in Google Colab for reproducible results online
# !apt install ffmpeg

def spring_xy(left_bound, right_bound, vertical_shift=0.25, harmonic=15, amplitude=0.1):
    '''
    Returns array of x and y coordinates of spring animation object
    as a standing wave with wave function 
        y = A sin(2pi(x-x0)/lambda) + b 
    where 
        lambda = 2L/n
    '''
    interval_len = right_bound - left_bound
    wave_len = 2*(interval_len)/harmonic
    spring_x = np.linspace(left_bound, right_bound, 100)
    spring_y = amplitude*np.sin(2*np.pi*(spring_x - left_bound)/wave_len) + vertical_shift
    return np.array([spring_x, spring_y])

# animation object params
block_width = 0
block_height = 0.5
block_center_offset = block_width / 2

wall1_x = 0
block1_eqm_x = wall1_x + l1
block2_eqm_x = block1_eqm_x + l2
wall2_x = block2_eqm_x + l3

# ## interactive plot 
# ## uncomment for testing
# ## used for checking time series before making animation
# ## can be used to test all animation below

# # plot settings
# custom_plot_settings()
# %matplotlib qt

# # plot static portion
# fig, ax = plt.subplots(figsize=(5, 2))
# ax.set_xlim(wall1_x, wall2_x)
# ax.set_ylim(0, 1)
# ax.spines['right'].set_visible(True)
# plt.tight_layout()

# # plot empty framework
# block1 = ax.add_patch(patches.Rectangle((block1_eqm_x+block_center_offset, 0), block_width, block_height, fill=False, color='black', zorder=10, lw=2)) 
# block2 = ax.add_patch(patches.Rectangle((block2_eqm_x+block_center_offset, 0), block_width, block_height, fill=False, color='black', zorder=10, lw=2)) 
# spring1, = ax.plot([], [], color='black')
# spring2, = ax.plot([], [], color='purple')
# spring3, = ax.plot([], [], color='black')
# title = ax.set_title('')

# # animation parameters
# t_range = t_len  # manually set t range, default t_len
# anim_time = 10  # s
# fps = 60
# frame_num = int(fps * anim_time)

# # update changes each frame
# for n in range(frame_num):
#     time_points = round(t_range/frame_num)
#     frame_final_time = min(time_points*n+time_points, t_range-1)  # avoid index out of range
#     title.set_text(f't = {t[frame_final_time] :.2f}')
#     # spring 1
#     left_bound = wall1_x
#     right_bound = block1_eqm_x + x1[frame_final_time] - block_center_offset
#     spring1.set_data(*spring_xy(left_bound, right_bound))
#     # spring 2
#     left_bound = block1_eqm_x + x1[frame_final_time] + block_center_offset
#     right_bound = block2_eqm_x + x2[frame_final_time] - block_center_offset
#     spring2.set_data(*spring_xy(left_bound, right_bound))
#     # spring 3
#     left_bound = block2_eqm_x + x2[frame_final_time] + block_center_offset
#     right_bound = wall2_x
#     spring3.set_data(*spring_xy(left_bound, right_bound))
#     # blocks
#     block1.set_xy([block1_eqm_x + x1[frame_final_time] - block_center_offset, 0])
#     block2.set_xy([block2_eqm_x + x2[frame_final_time] - block_center_offset, 0])
#     plt.pause(0.0001)

def make_animation(t_range=t_len, anim_time=4, fps=60):
    '''
    This function is notebook-specific and not meant to generalize to other settings.
    Makes animation of coupled harmonic oscillator with damping.
    Warning: Many parameters are taken from the global namespace. They need to be defined before use.
    '''
    # back to static plot and animations
    custom_plot_settings()

    # plot static portion
    fig, ax = plt.subplots(figsize=(5, 2))
    ax.set_xlim(wall1_x, wall2_x)
    ax.set_ylim(0, 1)
    ax.spines['right'].set_visible(True)
    plt.tight_layout()

    # plot empty framework
    block1 = ax.add_patch(patches.Rectangle((block1_eqm_x+block_center_offset, 0), block_width, block_height, fill=False, color='black', zorder=10, lw=2)) 
    block2 = ax.add_patch(patches.Rectangle((block2_eqm_x+block_center_offset, 0), block_width, block_height, fill=False, color='black', zorder=10, lw=2)) 
    spring1, = ax.plot([], [], color='black')
    spring2, = ax.plot([], [], color='purple')
    spring3, = ax.plot([], [], color='black')
    title = ax.set_title('')

    def draw_frame(n):
        '''
        Commands to update parameters.
        Here, the spring and block objects and the title.
        '''
        time_points = round(t_range/frame_num)
        frame_final_time = min(time_points*n+time_points, t_range-1)  # avoid index out of range
        title.set_text(f't = {t[frame_final_time] :.2f}')
        # spring 1
        left_bound = wall1_x
        right_bound = block1_eqm_x + x1[frame_final_time] - block_center_offset
        spring1.set_data(*spring_xy(left_bound, right_bound))
        # spring 2
        left_bound = block1_eqm_x + x1[frame_final_time] + block_center_offset
        right_bound = block2_eqm_x + x2[frame_final_time] - block_center_offset
        spring2.set_data(*spring_xy(left_bound, right_bound))
        # spring 3
        left_bound = block2_eqm_x + x2[frame_final_time] + block_center_offset
        right_bound = wall2_x
        spring3.set_data(*spring_xy(left_bound, right_bound))
        # blocks
        block1.set_xy([block1_eqm_x + x1[frame_final_time] - block_center_offset, 0])
        block2.set_xy([block2_eqm_x + x2[frame_final_time] - block_center_offset, 0])
        return fig,

    # create animation of given time length
    # note here we fit all the data points into the given animation time
    from matplotlib import animation
    frame_num = int(fps * anim_time)
    anim = animation.FuncAnimation(fig, draw_frame, frames=frame_num, interval=1000/fps, blit=True)

    plt.close()  # disable showing initial frame
    return anim

# convert animation to video (time-limiting step)
from IPython.display import HTML
anim = make_animation()  # uses custom function above
HTML(anim.to_html5_video() + '<style>video{width: 400px !important; height: auto;}</style>')

Internet

message='Mr. Trey'
messageBinary=[0]*len(message)

#turning string into ascii
for n in range(len(message)):
    messageBinary[n]=bin(ord(message[n]))

#removing ob
for n in range(len(message)):
    messageBinary[n]=messageBinary[n][2:]
    if len(messageBinary[n])<7:
        messageBinary[n]='0'+messageBinary[n]
    
messageToSend=''.join(messageBinary)
print('Message in binary sent '+messageToSend+'\n')

Message in binary sent 10011011110010010111001000001010100111001011001011111001

# model params
m1 = 1  # kg
m2 = 1  # kg
k1 = 100  # N/m
k2 = 100  # N/m
k3 = 100  # N/m
c1 = 3  # kg/s
c2 = 3  # kg/s
l1 = 1  # m
l2 = 1  # m
l3 = 1  # m

# time array
t_initial = 0
t_final = 3*len(messageToSend) - 1
t = np.linspace(t_initial, t_final, 1000)
t_len = len(t)

# initial values (x is relative to eqm length)
x1_init = 0
v1_init = 0
x2_init = 0
v2_init = 0
initial_values = [x1_init, v1_init, x2_init, v2_init]

def F1(t):
    if(messageToSend[int(t/3)] == '0'):
        omega = 1e4
    else:
        omega = np.sqrt(k1 + k2)
    A = 20
    phi = 1
    return A*np.sin(omega*t + phi)

def coupled_ode_syst(t, Y):
    # youtube
    # t -> independent variable
    # Y -> functions evaluated at independent variable
    x1, v1, x2, v2 = Y
    derivatives = np.array([
        v1,
        (-k1*x1 - k2*(x1 - x2) - c1*v1 + F1(t))/m1,
        v2,
        (-k2*(x2 - x1) - k3*x2 - c2*v2)/m2,
    ])
    return derivatives

ode_soln = scipy.integrate.solve_ivp(coupled_ode_syst, [t_initial, t_final], initial_values, t_eval=t).y
x1, v1, x2, v2 = ode_soln

fig, axs = plt.subplots(2, 1, figsize=(10, 4), sharex=True)
axs[0].plot(t, x1, label='$x_1(t)$', alpha=1)
axs[0].plot(t, x2, label='$x_2(t)$', alpha=0.5)
axs[0].set_ylabel('$x(t)$')
# axs[0].set_ylim(0.5, 2.5)

axs[1].plot(t, v1, label='$v_1(t)$', alpha=1)
axs[1].plot(t, v2, label='$v_2(t)$', alpha=0.5)
axs[1].set_xlabel('$t$')
axs[1].set_ylabel('$v(t)$')
# axs[1].set_ylim(-3.5, 3.5)

for i in range(2):
    axs[i].set_xlim(t_initial, t_final)
    axs[i].plot([t_initial, t_final], [0, 0], '--', color='grey', lw=0.5, zorder=0)  # zero ref

# convert animation to video (time-limiting step)
from IPython.display import HTML
anim = make_animation()  # uses custom function above
HTML(anim.to_html5_video() + '<style>video{width: 400px !important; height: auto;}</style>')