Coupled Oscillators

This example demonstrates a system of two coupled damped harmonic oscillators using two ODE blocks.

The system consists of two masses connected by springs with damping. Each oscillator is coupled to the other through a coupling spring constant.

The equations of motion for the coupled system are:

\[ \begin{align}\begin{aligned}m_1 \ddot{x}_1 = -k_1 x_1 - c_1 \dot{x}_1 - F_e\\m_2 \ddot{x}_2 = -k_2 x_2 - c_2 \dot{x}_2 + F_e\end{aligned}\end{align} \]

With the force-coupling given as:

\[F_e = k_{12}(x_1 - x_2)\]

where \(m_1, m_2\) are the masses, \(k_1, k_2\) are the spring constants, \(c_1, c_2\) are the damping coefficients, and \(k_{12}\) is the coupling spring constant between the two oscillators, leading to external forces \(F_e\).

As a block diagram it would look like this:

Now let’s implement this system in PathSim:

First let’s import the Simulation and Connection classes and the required blocks from the block library:

import numpy as np
import matplotlib.pyplot as plt

# Apply PathSim docs matplotlib style for consistent, theme-friendly figures
plt.style.use('../pathsim_docs.mplstyle')

from pathsim import Simulation, Connection
from pathsim.blocks import ODE, Function, Scope

System Parameters

Next, let’s define the system parameters:

# Mass parameters
m1 = 1.0
m2 = 1.5

# Spring constants
k1 = 2.0
k2 = 3.0
k12 = 0.5  # coupling spring constant

# Damping coefficients
c1 = 0.05
c2 = 0.1

# Initial conditions [position, velocity]
x1_0 = np.array([2.0, 0.0])  # oscillator 1 starts displaced
x2_0 = np.array([0.0, 0.0])  # oscillator 2 starts at rest

Block Construction

Now we define the differential equations for each oscillator and create the ODE blocks:

Each ODE block takes a function that defines the right-hand side of the differential equation, an initial condition, and optionally a Jacobian for improved convergence with implicit solvers.

# Define the differential equation for oscillator 1
def osc1_func(x1, u, t):
    """
    x1 = [position, velocity] of oscillator 1
    u = [force] external force due to coupling
    """
    f_e = u[0] # external force

    dx1_dt = x1[1]  # velocity
    dv1_dt = (-k1*x1[0] - c1*x1[1] - f_e) / m1  # acceleration

    return np.array([dx1_dt, dv1_dt])

# Define the differential equation for oscillator 2
def osc2_func(x2, u, t):
    """
    x2 = [position, velocity] of oscillator 2
    u = [force] external force due to coupling
    """
    f_e = u[0] # external force

    dx2_dt = x2[1]  # velocity
    dv2_dt = (-k2*x2[0] - c2*x2[1] - f_e) / m2  # acceleration

    return np.array([dx2_dt, dv2_dt])

# Define function for coupling of the oscillators
def coupling_func(x1, x2):
    f = k12 * (x1 - x2)
    return f, -f

# Create the ODE blocks
osc1 = ODE(osc1_func, x1_0)
osc2 = ODE(osc2_func, x2_0)

# Create Function block for coupling
fn = Function(coupling_func)

# Create a scope to visualize both oscillators
sc1 = Scope(labels=[r"$x_1(t)$ - Oscillator 1", r"$x_2(t)$ - Oscillator 2"])
sc2 = Scope(labels=[r"$f_e(t)$ - Coupling"])

blocks = [osc1, osc2, fn, sc1, sc2]

Connections

Now we connect the blocks. The key aspect of this system is the coupling between the two oscillators:

The Connection class defines the signal flow between blocks. Each oscillator sends its state to the other oscillator as input, creating the coupling. We also connect both oscillators to the Scope for visualization.

connections = [
    Connection(osc1[0], fn[0], sc1[0]),
    Connection(osc2[0], fn[1], sc1[1]),
    Connection(fn[0], osc1[0], sc2[0]),
    Connection(fn[1], osc2[0]),
]

Simulation Setup

Finally, we create the simulation:

We instantiate the Simulation with the blocks and connections. For this non-stiff system, the default SSPRK22 solver (a 2nd order explicit Runge-Kutta method) works well.

sim = Simulation(blocks, connections, dt=0.01, log=True)

10:57:45 - INFO - LOGGING (log: True)
10:57:45 - INFO - BLOCKS (total: 5, dynamic: 2, static: 3, eventful: 0)
10:57:45 - INFO - GRAPH (nodes: 5, edges: 7, alg. depth: 2, loop depth: 0, runtime: 0.068ms)

Running the Simulation

Now let’s run the simulation and visualize the results:

# Run the simulation
sim.run(duration=75, reset=True)

# Plot the results
fig, ax = sc1.plot()
fig, ax = sc2.plot()
plt.show()

10:57:45 - INFO - RESET (time: 0.0)
10:57:45 - INFO - STARTING -> TRANSIENT (Duration: 75.00s)
10:57:45 - INFO - --------------------   1% | 0.0s<0.9s | 8567.0 it/s
10:57:46 - INFO - ####----------------  20% | 0.2s<0.7s | 8857.1 it/s
10:57:46 - INFO - ########------------  40% | 0.4s<0.5s | 8606.0 it/s
10:57:46 - INFO - ############--------  60% | 0.5s<0.4s | 8493.4 it/s
10:57:46 - INFO - ################----  80% | 0.7s<0.2s | 8166.1 it/s
10:57:46 - INFO - #################### 100% | 0.9s<--:-- | 8389.4 it/s
10:57:46 - INFO - FINISHED -> TRANSIENT (total steps: 7500, successful: 7500, runtime: 908.44 ms)

Analysis

The plot shows the position of both oscillators over time. Notice how:

1. Energy Transfer: The initially displaced oscillator 1 transfers energy to oscillator 2 through the coupling spring.

2. Damped Motion: Both oscillators gradually lose energy due to damping, eventually settling to rest.

3. Coupled Dynamics: The motion of each oscillator is influenced by the other, creating a complex interplay.

This example demonstrates how PathSim can elegantly handle coupled differential equations using multiple ODE blocks that exchange information through connections, making it easy to model systems with interacting components.