API

resonance/system.py

class resonance.system.System

This is the abstract base class for any single or multi degree of freedom system. It can be sub-classed to make a custom system or the necessary methods can be added dynamically.

add_measurement(name, func)

Creates a new measurement entry in the measurements attribute that uses the provided function to compute the measurement given a subset of the constants, coordinates, speeds, other measurements, and time.

Parameters:
  • name (string) – This must be a valid Python variable name and it should not clash with any names in the constants, coordinates, or speeds dictionary. This string can be different that the function name.
  • func (function) – This function must only have existing constant, coordinate, speed, or measurement names, and/or the special name 'time' in the function signature. These can be a subset of the available choices in constants, coordinates, speeds, measurements and any order in the signature is permitted. The function must be able to operate on both inputs that are a collection of floats or a collection of equal length 1D NumPy arrays and floats, i.e. the function must be vectorized. So be sure to use NumPy vectorized functions inside your function, i.e. numpy.sin() instead of math.sin(). The measurement function you create should return a item, either a scalar or array, that gives the values of the measurement.

Examples

>>> import numpy as np
>>> from resonance.linear_systems import SingleDoFLinearSystem
>>> sys = SingleDoFLinearSystem()
>>> sys.constants['m'] = 1.0  # kg
>>> sys.constants['c'] = 0.2  # kg*s
>>> sys.constants['k'] = 10.0  # N/m
>>> sys.coordinates['x'] = 1.0  # m
>>> sys.speeds['v'] = 0.25  # m/s
>>> def can_coeffs(m, c, k):
...     return m, c, k
...
>>> sys.canonical_coeffs_func = can_coeffs
>>> def force(x, v, c, k, time):
...    return -k * x - c * v + 5.0 * time
...
>>> # The measurement function you create must be vectorized, such
>>> # that it works with both floats and 1D arrays. For example with
>>> # floats:
>>> force(1.0, 0.5, 0.2, 10.0, 0.1)
-9.6
>>> # And with 1D arrays:
>>> force(np.array([1.0, 1.0]), np.array([0.25, 0.25]), 0.2, 10.0,
...       np.array([0.1, 0.2]))
array([-9.55, -9.05])
>>> sys.add_measurement('f', force)
>>> sys.measurements['f']  # time is 0.0 by default
-10.05
>>> sys.constants['k'] = 20.0  # N/m
>>> sys.measurements['f']
-20.05
>>> # Note that you should use NumPy functions to ensure your
>>> # measurement is vectorized.
>>> def force_mag(force):
...     return np.abs(force)
...
>>> force_mag(-10.05)
10.050000000000001
>>> force_mag(np.array([-10.05, -20.05]))
array([ 10.05,  20.05])
>>> sys.add_measurement('fmag', force_mag)
>>> sys.measurements['fmag']
20.05
animate_configuration(fps=30, **kwargs)

Returns a matplotlib animation function based on the configuration plot and the configuration plot update function.

Parameters:
  • fps (integer) – The frames per second that should be displayed in the animation. The latest trajectory will be resampled via linear interpolation to create the correct number of frames. Note that the frame rate will depend on the CPU speed of the computer. You’ll likely have to adjust this by trial and error to get something that matches well for your computer if you want the animation to run in real time.
  • **kwargs – Any extra keyword arguments will be passed to matplotlib.animation.FuncAnimation(). The interval keyword argument will be ignored.
config_plot_func

The configuration plot function arguments should be any of the system’s constants, coordinates, measurements, or 'time'. No other arguments are valid. The function has to return the matplotlib figure as the first item but can be followed by any number of mutable matplotlib objects that you may want to change during an animation. Refer to the matplotlib documentation for tips on creating figures.

Examples

>>> sys = SingleDoFLinearSystem()
>>> sys.constants['radius'] = 5.0
>>> sys.constants['center_y'] = 10.0
>>> sys.coordinates['center_x'] = 0.0
>>> def plot(radius, center_x, center_y, time):
...     fig, ax = plt.subplots(1, 1)
...     circle = Circle((center_x, center_y), radius=radius)
...     ax.add_patch(circle)
...     ax.set_title(time)
...     return fig, circle, ax
...
>>> sys.config_plot_function = plot
>>> sys.plot_configuration()
config_plot_update_func

The configuration plot update function arguments should be any of the system’s constants, coordinates, measurements, or ‘time’ in any order with the returned values from the config_plot_func as the last arguments in the exact order as in the configuration plot return statement. No other arguments are valid. Nothing need be returned from the function. See the matplotlib animation documentation for tips on creating these update functions.

Examples

>>> sys = SingleDoFLinearSystem()
>>> sys.constants['radius'] = 5.0
>>> sys.constants['center_y'] = 10.0
>>> sys.coordinates['center_x'] = 0.0
>>> def plot(radius, center_x, center_y, time):
...     fig, ax = plt.subplots(1, 1)
...     circle = Circle((center_x, center_y), radius=radius)
...     ax.add_patch(circle)
...     ax.set_title(time)
...     return fig, circle, ax
...
>>> sys.config_plot_function = plot
>>> def update(center_y, center_x, time, circle, ax):
...     # NOTE : that circle and ax have to be the last arguments and be
...     # in the same order as returned from plot()
...     circle.set_xy((center_x, center_y))
...     ax.set_title(time)
...
>>> sys.config_plot_update_func = update
>>> sys.animate_configuration()
constants

A dictionary containing the all of the system’s constants, i.e. parameters that do not vary with time.

Examples

>>> sys = System()
>>> sys.constants
{}
>>> sys.constants['mass'] = 5.0
>>> sys.constants
{'mass': 5.0}
>>> del sys.constants['mass']
>>> sys.constants
{}
>>> sys.constants['mass'] = 5.0
>>> sys.constants['length'] = 10.0
>>> sys.constants
{'mass': 5.0, 'length': 10.0}
coordinates

A dictionary containing the system’s generalized coordinates, i.e. coordinate parameters that vary with time. These values will be used as initial conditions in simulations.

Examples

>>> sys = System()
>>> sys.coordinates['angle'] = 0.0
>>> sys.coordinates
{'angle': 0.0}
free_response(final_time, initial_time=0.0, sample_rate=100, integrator='rungakutta4', **kwargs)

Returns a data frame with monotonic time values as the index and columns for each coordinate and measurement at the time value for that row. Note that this data frame is stored on the system as the variable .result until this method is called again, which will overwrite it.

Parameters:
  • final_time (float) – A value of time in seconds corresponding to the end of the simulation.
  • initial_time (float, optional) – A value of time in seconds corresponding to the start of the simulation.
  • sample_rate (integer, optional) – The sample rate of the simulation in Hertz (samples per second). The time values will be reported at the initial time and final time, i.e. inclusive, along with times space equally based on the sample rate.
  • integrator (string, optional) – Either rungakutta4 or lsoda. The rungakutta4 option is a very simple implementation and the sample_rate directly affects the accuracy and quality of the result. The lsoda makes use of SciPy’s odeint function which switches between two integrators for stiff and non-stiff portions of the simulation and is variable step so the sample rate does not affect the quality and accuracy of the result. This has no affect on single degree of freedom linear systems, as their solutions are computed analytically.
Returns:

df – A data frame indexed by time with all of the coordinates and measurements as columns.

Return type:

pandas.DataFrame

measurements

A dictionary containing the all of the system’s measurements, i.e. parameters that are functions of the constants, coordinates, speeds, and other measurements.

plot_configuration()

Returns a matplotlib figure generated by the function assigned to the config_plot_func attribute. You may need to call matplotlib.pyplot.show() to display the figure.

Returns:
  • fig (matplotlib.figure.Figure) – The first returned object is always a figure.
  • *args (matplotlib objects) – Any matplotlib objects can be returned after the figure.
speeds

A dictionary containing the system’s generalized speeds, i.e. speed parameters that vary with time. These values will be used as initial conditions in simulations.

Examples

>>> sys = System()
>>> sys.speeds['angle_vel'] = 0.0
>>> sys.speeds
{'angle_vel': 0.0}
states

An ordered dictionary containing the system’s state variables and values. The coordinates are always ordered before the speeds and the individual order of the values depends on the order they were added to coordinates and speeds.

Examples

>>> sys = System()
>>> sys.coordinates['angle'] = 0.2
>>> sys.speeds['angle_vel'] = 0.1
>>> sys.states
{'angle': 0.2, 'angle_vel': 0.1}
>>> list(sys.states.keys())
['angle', 'angle_vel']
>>> list(sys.states.values())
[0.2, 0.1]

resonance/linear_systems.py

class resonance.linear_systems.AutomobileLateralSystem

Bases: resonance.linear_systems.MultiDoFLinearSystem

class resonance.linear_systems.BallChannelPendulumSystem

Bases: resonance.linear_systems.MultiDoFLinearSystem

class resonance.linear_systems.BaseExcitationSystem

Bases: resonance.linear_systems.SingleDoFLinearSystem

This system represents a mass connected to a moving massless base via a spring and damper in parallel. The motion of the mass is subject to viscous damping. The system is described by:

constants
mass, m [kg]
The suspended mass.
damping, c [kg / s]
The viscous linear damping coefficient which represents any energy dissipation from things like air resistance, friction, etc.
stiffness, k [N / m]
The linear elastic stiffness of the spring.
coordinates
position, x [m]
The absolute position of the mass.
speeds
velocity, x_dot [m / s]
The absolute velocity of the mass.
periodic_base_displacing_response(twice_avg, cos_coeffs, sin_coeffs, frequency, final_time, initial_time=0.0, sample_rate=100, force_col_name='forcing_function', displace_col_name='displacing_function')

Returns the trajectory of the system’s coordinates, speeds, accelerations, and measurements if a periodic function defined by a Fourier series is applied as displacement of the base in the same direction as the system’s coordinate. The displacing function is defined as:

                 N
y(t)  = a0 / 2 + ∑ (an * cos(n*ω*t) + bn * sin(n*ω*t))
                n=1

Where a0, a1…an, and b1…bn are the Fourier coefficients. If N=∞ then the Fourier series can describe any periodic function with a period (2*π)/ω.

Parameters:
  • twice_avg (float) – Twice the average value over one cycle, a0.
  • cos_coeffs (float or sequence of floats) – The N cosine Fourier coefficients: a1, …, aN.
  • sin_coeffs (float or sequence of floats) – The N sine Fourier coefficients: b1, …, bN.
  • frequency (float) – The frequency, ω, in radians per second corresponding to one full cycle of the function.
  • final_time (float) – A value of time in seconds corresponding to the end of the simulation.
  • initial_time (float, optional) – A value of time in seconds corresponding to the start of the simulation.
  • sample_rate (integer, optional) – The sample rate of the simulation in Hertz (samples per second). The time values will be reported at the initial time and final time, i.e. inclusive, along with times space equally based on the sample rate.
  • force_col_name (string, optional) – A valid Python identifier that will be used as the column name for the forcing function trajectory in the returned data frame.
  • displace_col_name (string, optional) – A valid Python identifier that will be used as the column name for the forcing function trajectory in the returned data frame.
Returns:

A data frame indexed by time with all of the coordinates, speeds, measurements, and forcing/displacing functions as columns.

Return type:

pandas.DataFrame

sinusoidal_base_displacing_response(amplitude, frequency, final_time, initial_time=0.0, sample_rate=100, force_col_name='forcing_function', displace_col_name='displacing_function')

Returns the trajectory of the system’s coordinates, speeds, accelerations, and measurements if a sinusoidal displacement function described by:

y(t) = Y * sin(ω*t)

is specified for the movement of the base in the direction of the system’s coordinate.

Parameters:
  • amplitude (float) – The amplitude of the displacement function, Y, in meters.
  • frequency (float) – The frequency, ω, in radians per second of the sinusoidal displacement.
  • final_time (float) – A value of time in seconds corresponding to the end of the simulation.
  • initial_time (float, optional) – A value of time in seconds corresponding to the start of the simulation.
  • sample_rate (integer, optional) – The sample rate of the simulation in Hertz (samples per second). The time values will be reported at the initial time and final time, i.e. inclusive, along with times space equally based on the sample rate.
  • force_col_name (string, optional) – A valid Python identifier that will be used as the column name for the forcing function trajectory in the returned data frame.
  • displace_col_name (string, optional) – A valid Python identifier that will be used as the column name for the forcing function trajectory in the returned data frame.
Returns:

A data frame indexed by time with all of the coordinates and measurements as columns.

Return type:

pandas.DataFrame

class resonance.linear_systems.BicycleSystem

Bases: resonance.linear_systems.MultiDoFLinearSystem

class resonance.linear_systems.BookOnCupSystem

Bases: resonance.linear_systems.SingleDoFLinearSystem

This system represents dynamics of a typical engineering textbook set atop a cylinder (a coffee cup) such that the book can vibrate without slip on the curvature of the cup. It is described by:

constants
thickness, t [meters]
the thickness of the book
length, l [meters]
the length of the edge of the book which is tangent to the cup’s surface
mass, m [kilograms]
the mass of the book
radius, r [meters]
the outer radius of the cup
coordinates
book_angle, theta [radians]
the angle of the book with respect to the gravity vector
speeds
book_angle_vel, theta [radians]
the angular rate of the book with respect to the gravity vector
class resonance.linear_systems.ClockPendulumSystem

Bases: resonance.linear_systems.SingleDoFLinearSystem

This system represents dynamics of a simple compound pendulum in which a rigid body is attached via a revolute joint to a fixed point. Gravity acts on the pendulum to bring it to an equilibrium state and there is no friction in the joint. It is described by:

constants
pendulum_mass, m [kg]
The mass of the compound pendulum.
inertia_about_joint, i [kg m**2]
The moment of inertia of the compound pendulum about the revolute joint.
joint_to_mass_center, l [m]
The distance from the revolute joint to the mass center of the compound pendulum.
acc_due_to_gravity, g [m/s**2]
The acceleration due to gravity.
coordinates
angle, theta [rad]
The angle of the pendulum relative to the direction of gravity. When theta is zero the pendulum is hanging down in it’s equilibrium state.
speeds
angle_vel, theta_dot [rad / s]
The angular velocity of the pendulum about the revolute joint axis.
class resonance.linear_systems.CompoundPendulumSystem

Bases: resonance.linear_systems.SingleDoFLinearSystem

This system represents dynamics of a simple compound pendulum in which a rigid body is attached via a revolute joint to a fixed point. Gravity acts on the pendulum to bring it to an equilibrium state and there is no friction in the joint. It is described by:

constants
pendulum_mass, m [kg]
The mass of the compound pendulum.
inertia_about_joint, i [kg m**2]
The moment of inertia of the compound pendulum about the revolute joint.
joint_to_mass_center, l [m]
The distance from the revolute joint to the mass center of the compound pendulum.
acc_due_to_gravity, g [m/s**2]
The acceleration due to gravity.
coordinates
angle, theta [rad]
The angle of the pendulum relative to the direction of gravity. When theta is zero the pendulum is hanging down in it’s equilibrium state.
speeds
angle_vel, theta_dot [rad / s]
The angular velocity of the pendulum about the revolute joint axis.
class resonance.linear_systems.FourStoryBuildingSystem

Bases: resonance.linear_systems.MultiDoFLinearSystem

class resonance.linear_systems.MassSpringDamperSystem

Bases: resonance.linear_systems.SingleDoFLinearSystem

This system represents dynamics of a mass connected to a spring and damper (dashpot). The mass moves horizontally without friction and is acted on horizontally by the spring and damper in parallel. The system is described by:

constants
mass, M [kg]
The system mass.
damping, C [kg / s]
The viscous linear damping coefficient which represents any energy dissipation from things like air resistance, slip, etc.
stiffness, K [N / m]
The linear elastic stiffness of the spring.
coordinates

position, x [m]

speeds

velocity, x_dot [m / s]

class resonance.linear_systems.MultiDoFLinearSystem

Bases: resonance.nonlinear_systems.MultiDoFNonLinearSystem

This is the abstract base class for any multi degree of freedom linear system. It can be sub-classed to make a custom system or the necessary methods can be added dynamically.

canonical_coefficients()

Returns the mass, damping, and stiffness matrices in that order.

canonical_coeffs_func

A function that returns the three linear coefficient matrices of the left hand side of a set of canonical second order ordinary differential equations. This equation looks like the following:

Mv’ + Cv + Kx = F(t)

where:

  • M: mass matrix
  • C: damping matrix
  • K: stiffness matrix
  • x: the generalized coordinate vector
  • v: the generalized speed vector

The coefficients M, C, and K must be defined in terms of the system’s constants.

Example

This is an example of a simple double pendulum linearized about its equilibrium.

>>> sys = MulitDoFLinearSystem()
>>> sys.constants['g'] = 9.8  # m/s**2
>>> sys.constants['l1'] = 1.0  # m
>>> sys.constants['l2'] = 1.0  # m
>>> sys.constnats['m1'] = 0.5  # kg
>>> sys.constnats['m2'] = 0.5  # kg
>>> sys.coordinates['theta1'] = 0.3  # rad
>>> sys.coordinates['theta2'] = 0.0  # rad
>>> sys.speeds['omega1'] = 0.0  # rad/s
>>> sys.speeds['omega2'] = 0.0  # rad/s
>>> def coeffs(m1, m2, l1, l2, g):
...     # Represents a linear model of a simple double pendulum
...     M = np.array([[l1 * (m1 + m2), m2 * l2],
...                   [m2 * l2, m2 * l1]])
...     C = 0.0
...     K = np.array([[-g * (m1 + m2), 0],
...                   [0, -m2 * g]])
...     return M, C, K
>>> sys.canonical_coeffs_func = coeffs
forced_response(final_time, initial_time=0.0, sample_rate=100, integrator='rungakutta4', **kwargs)

Returns a data frame with monotonic time values as the index and columns for each coordinate and measurement at the time value for that row. Note that this data frame is stored on the system as the variable result until this method is called again, which will overwrite it.

Parameters:
  • final_time (float) – A value of time in seconds corresponding to the end of the simulation.
  • initial_time (float, optional) – A value of time in seconds corresponding to the start of the simulation.
  • sample_rate (integer, optional) – The sample rate of the simulation in Hertz (samples per second). The time values will be reported at the initial time and final time, i.e. inclusive, along with times space equally based on the sample rate.
  • integrator (string, optional) – Either rungakutta4 or lsoda. The rungakutta4 option is a very simple implementation and the sample_rate directly affects the accuracy and quality of the result. The lsoda makes use of SciPy’s odeint function which switches between two integrators for stiff and non-stiff portions of the simulation and is variable step so the sample rate does not affect the quality and accuracy of the result. This has no affect on single degree of freedom linear systems, as their solutions are computed analytically.
Returns:

df – A data frame indexed by time with all of the coordinates and measurements as columns.

Return type:

pandas.DataFrame

Notes

You must have defined a forcing_func for this to execute. If there is no forcing function this will return the free response.

forcing_func

A function that returns the right hand side forcing vector of the canonical second order linear ordinary differential equations. This equation looks like the following:

Mv’ + Cv + Kx = F(t)

where:

  • M: mass matrix
  • C: damping matrix
  • K: stiffness matrix
  • x: the generalized coordinate vector
  • v: the generalized speed vector

The coefficients M, C, and K must be defined in terms of the system’s constants.

Example

This is an example of a simple double pendulum linearized about its equilibrium. The angles, theta1 and theta2, are defined relative to the vertical and when both are zero the pendulum is in its hanging equilibrium. The forcing function applies sinusoidal torquing with respect to theta1 and theta2.

>>> sys = MulitDoFLinearSystem()
>>> sys.constants['g'] = 9.8  # m/s**2
>>> sys.constants['l1'] = 1.0  # m
>>> sys.constants['l2'] = 1.0  # m
>>> sys.constants['m1'] = 0.5  # kg
>>> sys.constants['m2'] = 0.5  # kg
>>> sys.coordinates['theta1'] = 0.3  # rad
>>> sys.coordinates['theta2'] = 0.0  # rad
>>> sys.speeds['omega1'] = 0.0  # rad/s
>>> sys.speeds['omega2'] = 0.0  # rad/s
>>> def coeffs(m1, m2, l1, l2, g):
...     # Represents a linear model of a simple double pendulum
...     M = np.array([[l1 * (m1 + m2), m2 * l2],
...                   [m2 * l2, m2 * l1]])
...     C = np.zeros_like(M)
...     K = np.array([[-g * (m1 + m2), 0],
...                   [0, -m2 * g]])
...     return M, C, K
>>> sys.canonical_coeffs_func = coeffs
>>> sys.constants['To'] = 1.0  # Nm
>>> sys.constants['beta'] = 0.01  # rad/s
>>> def forcing(To, beta, time):
...     T1 = To * np.cos(beta * time)
...     T2 = To * np.sin(beta * time)
...     return T1, T2
...
>>> sys.forcing_func = forcing
class resonance.linear_systems.SimplePendulumSystem

Bases: resonance.linear_systems.SingleDoFLinearSystem

This system represents dynamics of a simple pendulum in which a point mass is fixed on a massless pendulum arm of some length to a revolute joint. Gravity acts on the pendulum to bring it to an equilibrium state and there is no friction in the joint. It is described by:

constants
pendulum_mass, m [kg]
The mass of the compound pendulum.
pendulum_length, l [m]
The distance from the revolute joint to the point mass location.
acc_due_to_gravity, g [m/s**2]
The acceleration due to gravity.
coordinates
angle, theta [rad]
The angle of the pendulum relative to the direction of gravity. When theta is zero the pendulum is hanging down in it’s equilibrium state.
speeds
angle_vel, theta_dot [rad / s]
The angular velocity of the pendulum about the revolute joint axis.
class resonance.linear_systems.SimpleQuarterCarSystem

Bases: resonance.linear_systems.BaseExcitationSystem

This system represents a mass connected to a moving massless base via a spring and damper in parallel. The motion of the mass is subject to viscous damping. The system is described by:

constants
mass, m [kg]
The suspended mass.
damping, c [kg / s]
The viscous linear damping coefficient which represents any energy dissipation from things like air resistance, friction, etc.
stiffness, k [N / m]
The linear elastic stiffness of the spring.
coordinates
position, x [m]
The absolute position of the mass.
speeds
velocity, x_dot [m / s]
The absolute velocity of the mass.
class resonance.linear_systems.SingleDoFLinearSystem

Bases: resonance.linear_systems._LinearSystem

This is the abstract base class for any single degree of freedom linear system. It can be sub-classed to make a custom system or the necessary methods can be added dynamically.

frequency_response(frequencies, amplitude)

Returns the amplitude and phase shift for simple sinusoidal forcing of the system. The first holds the plot of the coordinate’s amplitude as a function of forcing frequency and the second holds a plot of the coordinate’s phase shift with respect to the forcing function.

Parameters:
  • frequencies (array_like, shape(n,)) –
  • amplitude (float) – The value of the forcing amplitude.
Returns:

  • amp_curve (ndarray, shape(n,)) – The amplitude values of the coordinate at different frequencies.
  • phase_curve (ndarray, shape(n,)) – The phase shift values in radians of the coordinate relative to the forcing.

frequency_response_plot(amplitude, log=False, axes=None)

Returns an array of two matplotlib axes. The first holds the plot of the coordinate’s amplitude as a function of forcing frequency and the second holds a plot of the coordinate’s phase shift with respect to the forcing function.

Parameters:
  • amplitude (float) – The value of the forcing amplitude.
  • log (boolean, optional) – If True, the amplitude will be plotted on a semi-log Y plot.
period()

Returns the (damped) period of oscillation of the coordinate in seconds.

periodic_forcing_response(twice_avg, cos_coeffs, sin_coeffs, frequency, final_time, initial_time=0.0, sample_rate=100, col_name='forcing_function')

Returns the trajectory of the system’s coordinates, speeds, accelerations, and measurements if a periodic forcing function defined by a Fourier series is applied as a force or torque in the same direction as the system’s coordinate. The forcing function is defined as:

                        N
F(t) or T(t) = a0 / 2 + ∑ (an * cos(n*ω*t) + bn * sin(n*ω*t))
                       n=1

Where a0, a1…an, and b1…bn are the Fourier coefficients. If N=∞ then the Fourier series can describe any periodic function with a period (2*π)/ω.

Parameters:
  • twice_avg (float) – Twice the average value over one cycle, a0.
  • cos_coeffs (float or sequence of floats) – The N cosine Fourier coefficients: a1, …, aN.
  • sin_coeffs (float or sequence of floats) – The N sine Fourier coefficients: b1, …, bN.
  • frequency (float) – The frequency, ω, in radians per second corresponding to one full cycle of the function.
  • final_time (float) – A value of time in seconds corresponding to the end of the simulation.
  • initial_time (float, optional) – A value of time in seconds corresponding to the start of the simulation.
  • sample_rate (integer, optional) – The sample rate of the simulation in Hertz (samples per second). The time values will be reported at the initial time and final time, i.e. inclusive, along with times space equally based on the sample rate.
  • col_name (string, optional) – A valid Python identifier that will be used as the column name for the forcing function trajectory in the returned data frame.
Returns:

A data frame indexed by time with all of the coordinates and measurements as columns.

Return type:

pandas.DataFrame

sinusoidal_forcing_response(amplitude, frequency, final_time, initial_time=0.0, sample_rate=100, col_name='forcing_function')

Returns the trajectory of the system’s coordinates, speeds, accelerations, and measurements if a sinusoidal forcing (or torquing) function defined by:

F(t) = Fo * cos(ω * t)

or

T(t) = To * cos(ω * t)

is applied to the moving body in the direction of the system’s coordinate.

Parameters:
  • amplitude (float) – The amplitude of the forcing/torquing function, Fo or To, in Newtons or Newton-Meters.
  • frequency (float) – The frequency, ω, in radians per second of the sinusoidal forcing.
  • final_time (float) – A value of time in seconds corresponding to the end of the simulation.
  • initial_time (float, optional) – A value of time in seconds corresponding to the start of the simulation.
  • sample_rate (integer, optional) – The sample rate of the simulation in Hertz (samples per second). The time values will be reported at the initial time and final time, i.e. inclusive, along with times space equally based on the sample rate.
  • col_name (string, optional) – A valid Python identifier that will be used as the column name for the forcing function trajectory in the returned data frame.
Returns:

A data frame indexed by time with all of the coordinates and measurements as columns.

Return type:

pandas.DataFrame

class resonance.linear_systems.TorsionalPendulumSystem

Bases: resonance.linear_systems.SingleDoFLinearSystem

This system represents dynamics of a simple torsional pendulum in which the torsionally elastic member’s axis is aligned with gravity and the axis of the torsion member passes through the mass center of an object attached to it’s lower end. The top of the torsion rod is rigidly attached to the “ceiling”. It is described by:

constants
rotational_inertia, I [kg m**2]
The moment of inertia of the object attached to the pendulum.
torsional_damping, C [N s / m]
The viscous linear damping coefficient which represents any energy dissipation from things like air resistance, slip, etc.
torsional_stiffness, K [N / m]
The linear elastic stiffness coefficient of the torsion member, typically a round slender rod.
coordinates

torsional_angle, theta [rad]

speeds

torsional_angle_vel, theta_dot [rad / s]

resonance/nonlinear_systems.py

class resonance.nonlinear_systems.BallChannelPendulumSystem

Bases: resonance.nonlinear_systems.MultiDoFNonLinearSystem

class resonance.nonlinear_systems.ClockPendulumSystem

Bases: resonance.nonlinear_systems.SingleDoFNonLinearSystem

This system represents dynamics of a compound pendulum representing a clock pendulum. It is made up of a thin long cylindrical rod with a thin disc bob on the end. Gravity acts on the pendulum to bring it to an equilibrium state and there is option Coulomb friction in the joint. It is described by:

constants
bob_mass, m_b [kg]
The mass of the bob (a thin disc) on the end of the pendulum.
bob_radius, r [m]
The radius of the bob (a thin disc) on the end of the pendulum.
rod_mass, m_r [kg]
The mass of the then cylindrical rod.
rod_length, l [m]
The length of the rod which connects the pivot joint to the center of the bob.
coeff_of_friction, mu [unitless]
The Coulomb coefficient of friction between the materials of the pivot joint.
joint_friction_radius, R [m]
The radius of the contact disc at the pivot joint. The joint is assumed to be two flat discs pressed together.
joint_clamp_force, F_N [N]
The clamping force pressing the two flat discs together at the pivot joint.
acc_due_to_gravity, g [m/s**2]
The acceleration due to gravity.
coordinates
angle, theta [rad]
The angle of the pendulum relative to the direction of gravity. When theta is zero the pendulum is hanging down in it’s equilibrium state.
speeds
angle_vel, theta_dot [rad / s]
The angular velocity of the pendulum about the revolute joint axis.
class resonance.nonlinear_systems.MultiDoFNonLinearSystem

Bases: resonance.system.System

This is the abstract base class for any single degree of freedom nonlinear system. It can be sub-classed to make a custom system or the necessary methods can be added dynamically.

diff_eq_func

A function that returns the time derivatives of the coordinates and speeds, i.e. computes the right hand side of the explicit first order differential equations. This equation looks like the following for linear motion:

dx
-- = f(t, q1, ..., qn, u1, ..., un, p1, p2, ..., pO)
dt

where:

  • x: [q1, …, qn, u1, …, un], the “state vector”
  • t: a time value
  • q: the coordinates
  • u: the speeds
  • p: any number of constants, O is the number of constants

Your function should be able to operate on 1d arrays as inputs, i.e. use numpy math functions in your function, e.g. numpy.sin instead of math.sin. Besides the constants, coordinates, and speeds, there is a special variable time that you can use to give the current value of time inside your function.

Note

The function has to return the derivatives of the states in the order of the state attribute.

Warning

Do not use measurements as a function argument. This may cause causality issues and is not yet supported. You are unlikely to get a correct answer if you use a measurement in this function.

Example

>>> sys = SingleDoFNonLinearSystem()
>>> sys.constants['gravity'] = 9.8  # m/s**2
>>> sys.constants['length'] = 1.0  # m
>>> sys.constants['mass'] = 0.5  # kg
>>> sys.constants['omega_b'] = 0.1  # rad/s
>>> sys.coordinates['theta'] = 0.3  # rad
>>> sys.speeds['omega'] = 0.0  # rad/s
>>> sys.states
{'theta': 0.3, 'omega': 0.0}  # note the order!
>>> def rhs(theta, omega, gravity, length, mass, omega_b, time):
...     # Represents a linear model of a simple pendulum under
...     # sinusoidal torquing.
...     #  m * l**2 ω' + m * g * l * sin(θ) = sin(ω_b * t)
...     thetad = omega
...     omegad = (np.sin(omega_b * time) -
...               m*g*l*np.sin(theta)) / m / l**2
...     return thetad, omegad  # in order of sys.states
>>> sys.diff_eq_func = rhs
class resonance.nonlinear_systems.SingleDoFNonLinearSystem

Bases: resonance.nonlinear_systems.MultiDoFNonLinearSystem

resonance/functions.py

class resonance.functions.Phasor(init, frequency=0, growth_rate=0)

Phasor that can be advanced in time with rotation and growth rates.

Parameters:
  • init (complex) – Initial phasor in rectangular form (Re + jIm)
  • frequency (float, optional) – Rotation rate in rad/s.
  • growth_rate (float, optional) – Exponential growth rate (decay if < 0).
t

Current time.

Type:float
re

Current real component of the phasor.

Type:float
im

Current imaginary component of the phasor.

Type:float
radius

Current radius of the phasor.

Type:float
angle

Current angle of the phasor.

Type:float
trace_t

History of time values (since most recent clear()).

Type:list
trace_re

History of real component values (since most recent clear()).

Type:list
trace_im

History of imaginary component values (since most recent clear()).

Type:list
advance(dt)

Advance the phasor by a time step dt.

clear()

Clear trajectories.

classmethod from_eig(eigvec_component, eigval)

Creates a phasor from an eigenvalue/eigenvector component pair.

Parameters:
  • eigvec_component (complex) – A single eigenvector component representing the phasor’s initial real/imaginary parts.
  • eigval (complex) – The eigenvector, which specifies the phasor’s growth rate (real part) and rotational frequency (imaginary part).
class resonance.functions.PhasorAnimation(fig, t, phasors, re_range=(-1, 1), im_range=(-1, 1), repeat=True, repeat_delay=0, time_stretch=1, blit=True)

Animation for demonstrating rotating phasors.

Two axes are set up. On top, there is an s-plane to show the real and imaginary components of the phasors. The current phasor “vector” is shown with a thick line, the current endpoint of the vector is shown with a circle, thin lines show the projection of the real part of the phasor down to the bottom of the plane, and the time history of the endpoint of the vectors are shown.

On bottom, the phasors’ real components are plotted in time. The plot is rotated so that time is positive downward, and the x axes of the s-plane and the time plots are lined up. The current value is shown with a circle, thin lines show the projection from the top of the plot to the current value, and the time history is plotted.

Parameters:
  • fig (Figure) – matplotlib Figure object on which to animate.
  • t (array) – Array of time values at which to plot. Even time spacing is assumed.
  • phasors (list) – List of Phasor objects to advance and plot.
  • re_range (tuple, optional) – Limits of the real axis.
  • im_range (tuple, optional) – Limits of the imaginary axis.
  • repeat (bool, optional) – Specifies whether or not to repeat the animation once it finishes.
  • repeat_delay (float, optional) – Amount of time to wait before repeating the animation in milliseconds.
  • time_stretch (float, optional) – Multiplicative factor of the plotting interval. Increasing time_stretch effectively makes the animation slower without affecting the time units.
  • blit (bool, optional) – Specifies whether or not to use blitting.
new_frame_seq()

Return a new sequence of frame information.

resonance.functions.benchmark_par_to_canonical(p)

Returns the canonical matrices of the Whipple bicycle model linearized about the upright constant velocity configuration. It uses the parameter definitions from [Meijaard2007].

Parameters:p (dictionary) – A dictionary of the benchmark bicycle parameters. Make sure your units are correct, best to ue the benchmark paper’s units!
Returns:
  • M (ndarray, shape(2,2)) – The mass matrix.
  • C1 (ndarray, shape(2,2)) – The damping like matrix that is proportional to the speed, v.
  • K0 (ndarray, shape(2,2)) – The stiffness matrix proportional to gravity, g.
  • K2 (ndarray, shape(2,2)) – The stiffness matrix proportional to the speed squared, v**2.

References

[Meijaard2007]J. P. Meijaard, J. M. Papadopoulos, A. Ruina, and A. L. Schwab, “Linearized dynamics equations for the balance and steer of a bicycle: A benchmark and review,” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 463, no. 2084, pp. 1955–1982, Aug. 2007.
resonance.functions.centered_rectangle(xy, width, height, angle=0.0)

Returns the arguments for Rectangle given the x and y coordinates of the center of the rectangle.

Parameters:
  • xy (tuple of floats) – The x and y coordinates of the center of the rectangle.
  • width (float) – Width of the rectangle. When angle=0.0 this is along the x axis.
  • height (float) – Height of the rectangle. When angle=0.0 this is along the y axis.
  • angle (float) – Angle of rotation about the z axis in degrees.
Returns:

  • xy_ll (tuple of floats) – The x and y coordinates of the lower left hand corner of the rectangle.
  • width (float) – Width of the rectangle. When angle=0.0 this is along the x axis.
  • height (float) – Height of the rectangle. When angle=0.0 this is along the y axis.
  • angle (float) – Angle of rotation about the z axis in degrees.

resonance.functions.estimate_period(time, signal)

Computes the period of oscillation based on the given periodic signal.

Parameters:
  • time (array_like, shape(n,)) – An array of monotonically increasing time values.
  • signal (array_like, shape(n,)) – An array of values for the periodic signal at each time in t.
Returns:

period – An estimate of the period of oscillation.

Return type:

float

resonance.functions.spring(xA, xB, yA, yB, w, n=1, x=None, y=None)

Returns the x and y coordinates of the points that define a spring diagram between points (xA, yB) and (yA, yB).

Parameters:
  • xA (float) – x coordinate of the beginning of the spring.
  • xB (float) – x coordinate of the end of the spring.
  • yA (float) – y coordinate of the beginning of the spring.
  • yB (float) – y coordinate of the end of the spring.
  • w (float) – The width of the spring.
  • n (integer, optional) – Number of coils.
  • x (ndarray, shape(2*n + 2), optional) – Preallocated array for the results.
  • y (ndarray, shape(2*n + 2), optional) – Preallocated array for the results.
Returns:

  • x (ndarray, shape(2*n + 2)) – x coordinates of the points that define the ends of each line in the spring.
  • y (ndarray, shape(2*n + 2)) – y coordinates of the points that define the ends of each line in the spring.

Examples

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from resonance.functions import spring
>>> plt.axes().set_aspect('equal')
>>> for angle in np.arange(0, 2*np.pi, np.pi/4):
...     plt.plot(*spring(0.0, np.cos(angle), 0.0, np.sin(angle), 0.1, n=4))
...
>>> plt.show()