Creating and Exercising a Custom Single Degree of Freedom System


You can download this example as a Python script: or Jupyter notebook: custom-sdof-system.ipynb.

Creating a new system

The first step is to import a “blank” SingleDoFLinearSystem and initialize it.

from resonance.linear_systems import SingleDoFLinearSystem

msd_sys = SingleDoFLinearSystem()

Now define the constant variables for the system. In this case, the single degree of freedom system will be described by its mass, natural frequency, and damping ratio.

msd_sys.constants['m'] = 1.0  # kg
msd_sys.constants['fn'] = 1.0  # Hz
msd_sys.constants['zeta'] = 0.1  # unitless

{'m': 1.0, 'fn': 1.0, 'zeta': 0.1}

Define the coordinate and speed. The software assumes that the speed is defined as the time derivative of the coordinate, i.e. \(v = \dot{x}\).

msd_sys.coordinates['x'] = 1.0  # m
msd_sys.speeds['v'] = 0.0  # m/s
{'x': 1.0}
{'v': 0.0}
{'x': 1.0, 'v': 0.0}

Now that the coordinate, speed, and constants are defined the equations of motion can be defined. For a single degree of freedom system a Python function must be defined that uses the system’s constants to compute the coefficients to the canonical second order equation, \(m \dot{v} + c v + k x = 0\).

The inputs to the function are the constants. The variable names should match those defined above on the system.

import numpy as np

def calculate_canonical_coefficients(m, fn, zeta):
    """Returns the system's mass, damping, and stiffness coefficients given
    the system's constants."""
    wn = 2*np.pi*fn
    k = m*wn**2
    c = zeta*2*wn*m
    return m, c, k

msd_sys.canonical_coeffs_func = calculate_canonical_coefficients

Once this function is defined and added to the system \(m,c,k\) can be computed using:

(1.0, 1.2566370614359172, 39.47841760435743)

The period of the natural frequency can be computed with:


All information about the system can be displayed:

System name: SingleDoFLinearSystem

Canonical coefficients function defined: True
Configuration plot function defined: False
Configuration update function defined: False

m = 1.00000
fn = 1.00000
zeta = 0.10000

x = 1.00000

v = d(x)/dt = 0.00000


Simulating the free response

The free_response() function simulates the now fully defined system given as an initial value problem. One or both of the coordinates and speeds must be set to provide a free response. The following shows the response to both \(x\) and \(v\) being set to some initial values.

msd_sys.coordinates['x'] = -5.0
msd_sys.speeds['v'] = 8.0

free_response() returns a Pandas DataFrame with the time values as the index and columns for the coordinate, speed, and additionally the time derivative of the speed (acceleration in this case). See for an introduction to DataFrame.

trajectories = msd_sys.free_response(5.0)
x x_acc v
0.00 -5.000000 187.338992 8.000000
0.01 -4.910728 181.496514 9.844723
0.02 -4.803311 175.015195 11.627801
0.03 -4.678398 167.928429 13.343009
0.04 -4.536697 160.271570 14.984469
... ... ... ...
4.96 -0.217071 8.839775 -0.214970
4.97 -0.218780 8.796369 -0.126761
4.98 -0.219609 8.719006 -0.039156
4.99 -0.219566 8.608415 0.047509
5.00 -0.218663 8.465444 0.132905

501 rows × 3 columns

There are a variety of plotting methods associated with the DataFrame that can be used to quickly plot the trajectories of the coordinate, speed, and acceleration. See more about plotting DataFrames at

axes = trajectories.plot(subplots=True)

Response to change in constants

This system is parameterized by its mass, natural frequency, and damping ratio. It can be useful to plot the trajectories of position for different values of \(\zeta\) for example.

Set the initial conditions back to simply stretching the spring 1 meter:

msd_sys.coordinates['x'] = 1.0
msd_sys.speeds['v'] = 0.0

Now change \(\zeta\) to different values and simulate the free response to see the different damping regimes:

Un-damped, \(\zeta=0\)

msd_sys.constants['zeta'] = 0.0  # Unitless
trajectories = msd_sys.free_response(5.0)
axes = trajectories['x'].plot()

Under-damped, \(0<\zeta<1\)

msd_sys.constants['zeta'] = 0.5  # Unitless
trajectories = msd_sys.free_response(5.0)
axes = trajectories['x'].plot()

Critically damped, \(\zeta=1\)

msd_sys.constants['zeta'] = 1.0  # Unitless
trajectories = msd_sys.free_response(5.0)
axes = trajectories['x'].plot()

Over-damped, \(\zeta>1\)

msd_sys.constants['zeta'] = 2.0  # Unitless
trajectories = msd_sys.free_response(5.0)
axes = trajectories['x'].plot()

Adding measurements

It is often useful to calculate the trajectories of other quantities. Systems in resonance allow “measurements” to be defined. These measurements are functions of the constants, coordinates, speeds, and/or time. To create a new measurement, create a function that returns the quantity of interest. Here a measurement function is defined that calculates the kinetic energy (\(\frac{1}{2}mv^2\)) of the system and then added to the system with variable name KE.

def calculate_kinetic_energy(m, v):
   return m*v**2/2

msd_sys.add_measurement('KE', calculate_kinetic_energy)

Once added, the measurement will be computed and added to the DataFrame containing the trajectories:

msd_sys.constants['zeta'] = 0.5  # Unitless
trajectories = msd_sys.free_response(5.0)
x x_acc v KE
0.00 1.000000e+00 -39.478418 8.881784e-16 3.944305e-31
0.01 9.980674e-01 -36.999522 -3.823857e-01 7.310941e-02
0.02 9.924348e-01 -34.530060 -7.400220e-01 2.738162e-01
0.03 9.833490e-01 -32.078904 -1.073048e+00 5.757160e-01
0.04 9.710551e-01 -29.654314 -1.381689e+00 9.545318e-01
... ... ... ... ...
4.96 4.648115e-08 0.000006 -1.189485e-06 7.074377e-13
4.97 3.487003e-08 0.000006 -1.132570e-06 6.413573e-13
4.98 2.383263e-08 0.000006 -1.074788e-06 5.775850e-13
4.99 1.337624e-08 0.000006 -1.016414e-06 5.165491e-13
5.00 3.505453e-09 0.000006 -9.577074e-07 4.586018e-13

501 rows × 4 columns

and can be plotted like any other column:

axes = trajectories['KE'].plot()

Plotting the configuration

resonance systems can plot and animate at the system’s configuration. To do so, a custom function that generates a configuration plot using matplotlib must be defined and associated with the system. Below a plot is created to show an orange block representing the mass and a spring attached to the block. The spring() function conveniently provides the x and y data needed to plot the spring.

import matplotlib.pyplot as plt
from resonance.functions import spring

# create a new constant to describe the block's dimension, l
msd_sys.constants['l'] = 0.2  # m

def create_configuration_figure(x, l):

    # create a figure with one or more axes
    fig, ax = plt.subplots()

    # the `spring()` function creates the x and y data for plotting a simple
    # spring
    spring_x_data, spring_y_data = spring(0.0, x, l/2, l/2, l/8, n=3)
    lines = ax.plot(spring_x_data, spring_y_data, color='purple')
    spring_line = lines[0]

    # add a square that represents the mass
    square = plt.Rectangle((x, 0.0), width=l, height=l, color='orange')

    # add a vertical line representing the spring's attachment point
    ax.axvline(0.0, linewidth=4.0, color='black')

    # set axis limits and aspect ratio such that the entire motion will appear
    ax.set_ylim((-l/2, 3*l/2))
    ax.set_xlim((-np.abs(x) - l, np.abs(x) + l))

    ax.set_xlabel('$x$ [m]')
    ax.set_ylabel('$y$ [m]')

    # this function must return the figure as the first item
    # but you also may return any number of objects that you'd like to have
    # access to modify, e.g. for an animation update

    return fig, ax, spring_line, square

# associate the function with the system
msd_sys.config_plot_func = create_configuration_figure

Now the configuration plot can be generated with plot_configuration(). This returns the same results as the function defined above.

fig, ax, spring_line, square = msd_sys.plot_configuration()

Animating the configuration

Reset to un-damped motion and simulate again

msd_sys.constants['zeta'] = 0.1
trajectories = msd_sys.free_response(5.0)

To animate the configuration, create a function that updates the various matplotlib objects using any constants, coordinates, speeds, and/or the special variable time. The last input arguments to this function must be all of the extra outputs of plot_configuration() (excluding the figure which is the first output). The order of these must match the order of the plot_configuration() outputs.

def update_configuration(x, l, time,  # any variables you need for updating
                         ax, spring_line, square):  # returned items from plot_configuration() in same order

    ax.set_title('{:1.2f} [s]'.format(time))

    xs, ys = spring(0.0, x, l/2, l/2, l/8, n=3)
    spring_line.set_data(xs, ys)

    square.set_xy((x, 0.0))

msd_sys.config_plot_update_func = update_configuration

Now that the update function is associated, animate_configuration() will create the animation. Here the frames-per-second are set to an explicit value.

animation = msd_sys.animate_configuration(fps=30)

If using the notebook interactively with %matplotlib widget set, the animation above will play. But animate_configuration() returns a matplotlib FuncAnimation object which has other options that allow the generation of different formats, see for options. One option is to create a Javascript/HTML versions that displays nicely in the notebook with different play options:

from IPython.display import HTML