# Creating and Exercising a Custom Single Degree of Freedom System¶

Note

You can download this example as a Python script: custom-sdof-system.py 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
msd_sys.constants
```

```
{'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
```

```
msd_sys.coordinates
```

```
{'x': 1.0}
```

```
msd_sys.speeds
```

```
{'v': 0.0}
```

```
msd_sys.states
```

```
{'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:

```
msd_sys.canonical_coefficients()
```

```
(1.0, 1.2566370614359172, 39.47841760435743)
```

The period of the natural frequency can be computed with:

```
msd_sys.period()
```

```
1.005037815259212
```

All information about the system can be displayed:

```
msd_sys
```

```
System name: SingleDoFLinearSystem
Canonical coefficients function defined: True
Configuration plot function defined: False
Configuration update function defined: False
Constants
=========
m = 1.00000
fn = 1.00000
zeta = 0.10000
Coordinates
===========
x = 1.00000
Speeds
======
v = d(x)/dt = 0.00000
Measurements
============
```

## 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
https://pandas.pydata.org/pandas-docs/stable/getting_started/dsintro.html for
an introduction to `DataFrame`

.

```
trajectories = msd_sys.free_response(5.0)
trajectories
```

x | x_acc | v | |
---|---|---|---|

time | |||

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
https://pandas.pydata.org/pandas-docs/stable/user_guide/visualization.html.

```
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)
trajectories
```

x | x_acc | v | KE | |
---|---|---|---|---|

time | ||||

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')
ax.add_patch(square)
# 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_aspect('equal')
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
https://matplotlib.org/api/_as_gen/matplotlib.animation.FuncAnimation.html 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
HTML(animation.to_jshtml(fps=30))
```