Chaos Pendulum

Note

You can download this example as a Python script: chaos-pendulum.py or Jupyter notebook: chaos-pendulum.ipynb.

This example gives a simple demostration of chaotic behavior in a simple two body system. The system is made up of a slender rod that is connected to the ceiling at one end with a revolute joint that rotates about the \(\hat{\mathbf{n}}_y\) unit vector. At the other end of the rod a flat plate is attached via a second revolute joint allowing the plate to rotate about the rod’s axis with aligns with the \(\hat{\mathbf{a}_z}\) unit vector.

../_images/chaos-pendulum.svg

Setup

import numpy as np
import matplotlib.pyplot as plt
import sympy as sm
import sympy.physics.mechanics as me
from pydy.system import System
from pydy.viz import Cylinder, Plane, VisualizationFrame, Scene

%matplotlib inline

me.init_vprinting(use_latex='mathjax')

Define Variables

First define the system constants:

  • \(m_A\): Mass of the slender rod.

  • \(m_B\): Mass of the plate.

  • \(l_B\): Distance from \(N_o\) to \(B_o\) along the slender rod’s axis.

  • \(w\): The width of the plate.

  • \(h\): The height of the plate.

  • \(g\): The acceleratoin due to gravity.

mA, mB, lB, w, h, g = sm.symbols('m_A, m_B, L_B, w, h, g')

There are two time varying generalized coordinates:

  • \(\theta(t)\): The angle of the slender rod with respect to the ceiling.

  • \(\phi(t)\): The angle of the plate with respect to the slender rod.

The two generalized speeds will then be defined as:

  • \(\omega(t)=\dot{\theta}\): The angular rate of the slender rod with respect to the ceiling.

  • \(\alpha(t)=\dot{\phi}\): The angluer rate of the plate with respect to the slender rod.

theta, phi = me.dynamicsymbols('theta, phi')
omega, alpha = me.dynamicsymbols('omega, alpha')

The kinematical differential equations are defined in this fashion for the KanesMethod class:

\[\begin{split}0 = \omega - \dot{\theta}\\ 0 = \alpha - \dot{\phi}\end{split}\]
kin_diff = (omega - theta.diff(), alpha - phi.diff())
kin_diff
\[\displaystyle \left( \omega - \dot{\theta}, \ \alpha - \dot{\phi}\right)\]

Define Orientations

There are three reference frames. These are defined as such:

N = me.ReferenceFrame('N')
A = me.ReferenceFrame('A')
B = me.ReferenceFrame('B')

The frames are oriented with respect to each other by simple revolute rotations. The following lines set the orientations:

A.orient(N, 'Axis', (theta, N.y))
B.orient(A, 'Axis', (phi, A.z))

Define Positions

Three points are necessary to define the problem:

  • \(N_o\): The fixed point which the slender rod rotates about.

  • \(A_o\): The center of mass of the slender rod.

  • \(B_o\): The center of mass of the plate.

No = me.Point('No')
Ao = me.Point('Ao')
Bo = me.Point('Bo')

The two centers of mass positions can be set relative to the fixed point, \(N_o\).

lA = (lB - h / 2) / 2
Ao.set_pos(No, lA * A.z)
Bo.set_pos(No, lB * A.z)

Specify the Velocities

The generalized speeds should be used in the definition of the linear and angular velocities when using Kane’s method. For simple rotations and the defined kinematical differential equations the angular rates are:

A.set_ang_vel(N, omega * N.y)
B.set_ang_vel(A, alpha * A.z)

Once the angular velocities are specified the linear velocities can be computed using the two point velocity thereom, starting with the origin point having a velocity of zero.

No.set_vel(N, 0)

Ao.v2pt_theory(No, N, A)
\[\displaystyle \left(\frac{L_{B}}{2} - \frac{h}{4}\right) \omega\mathbf{\hat{a}_x}\]
Bo.v2pt_theory(No, N, A)
\[\displaystyle L_{B} \omega\mathbf{\hat{a}_x}\]

Inertia

The central inertia of the symmetric slender rod with respect to its reference frame is a function of its length and its mass.

IAxx = sm.S(1) / 12 * mA * (2 * lA)**2
IAyy = IAxx
IAzz = 0

IA = (me.inertia(A, IAxx, IAyy, IAzz), Ao)

This gives the inertia tensor:

IA[0].to_matrix(A)
\[\begin{split}\displaystyle \left[\begin{matrix}\frac{m_{A} \left(L_{B} - \frac{h}{2}\right)^{2}}{12} & 0 & 0\\0 & \frac{m_{A} \left(L_{B} - \frac{h}{2}\right)^{2}}{12} & 0\\0 & 0 & 0\end{matrix}\right]\end{split}\]

The central inerita of the symmetric plate with respect to its reference frame is a function of its width and height.

IBxx = sm.S(1)/12 * mB * h**2
IByy = sm.S(1)/12 * mB * (w**2 + h**2)
IBzz = sm.S(1)/12 * mB * w**2

IB = (me.inertia(B, IBxx, IByy, IBzz), Bo)

IB[0].to_matrix(B)
\[\begin{split}\displaystyle \left[\begin{matrix}\frac{h^{2} m_{B}}{12} & 0 & 0\\0 & \frac{m_{B} \left(h^{2} + w^{2}\right)}{12} & 0\\0 & 0 & \frac{m_{B} w^{2}}{12}\end{matrix}\right]\end{split}\]

All of the information to define the two rigid bodies are now available. This information is used to create an object for the rod and the plate.

rod = me.RigidBody('rod', Ao, A, mA, IA)

plate = me.RigidBody('plate', Bo, B, mB, IB)

Loads

The only loads in this problem is the force due to gravity that acts on the center of mass of each body. These forces are specified with a tuple containing the point of application and the force vector.

rod_gravity = (Ao, mA * g * N.z)
plate_gravity = (Bo, mB * g * N.z)

Equations of motion

Now that the kinematics, kinetics, and inertia have all been defined the KanesMethod class can be used to generate the equations of motion of the system. In this case the independent generalized speeds, independent generalized speeds, the kinematical differential equations, and the inertial reference frame are used to initialize the class.

kane = me.KanesMethod(N, q_ind=(theta, phi), u_ind=(omega, alpha), kd_eqs=kin_diff)

The equations of motion are then generated by passing in all of the loads and bodies to the kanes_equations method. This produces \(f_r\) and \(f_r^*\).

bodies = (rod, plate)
loads = (rod_gravity, plate_gravity)

fr, frstar = kane.kanes_equations(bodies, loads)

sm.trigsimp(fr)
\[\begin{split}\displaystyle \left[\begin{matrix}g \left(- \frac{L_{B} m_{A}}{2} - L_{B} m_{B} + \frac{h m_{A}}{4}\right) \sin{\left(\theta \right)}\\0\end{matrix}\right]\end{split}\]
sm.trigsimp(frstar)
\[\begin{split}\displaystyle \left[\begin{matrix}\frac{m_{B} w^{2} \alpha \omega \sin{\left(2 \phi \right)}}{12} - \left(\frac{L_{B}^{2} m_{A}}{3} + L_{B}^{2} m_{B} - \frac{L_{B} h m_{A}}{3} + \frac{h^{2} m_{A}}{12} + \frac{h^{2} m_{B}}{12} + \frac{m_{B} w^{2} \cos^{2}{\left(\phi \right)}}{12}\right) \dot{\omega}\\- \frac{m_{B} w^{2} \left(\omega^{2} \sin{\left(2 \phi \right)} + 2 \dot{\alpha}\right)}{24}\end{matrix}\right]\end{split}\]

Simulation

The equations of motion can now be simulated numerically. Values for the constants, initial conditions, and time are provided to the System class along with the symbolic KanesMethod object.

sys = System(kane)

sys.constants = {lB: 0.2, # meters
                 h: 0.1, # meters
                 w: 0.2, # meters
                 mA: 0.01, # kilograms
                 mB: 0.1, # kilograms
                 g: 9.81} # meters per second squared

sys.initial_conditions = {theta: np.deg2rad(45),
                          phi: np.deg2rad(0.5),
                          omega: 0,
                          alpha: 0}

sys.times = np.linspace(0.0, 10.0, num=300)

The trajectories of the states are found with the integrate method.

x = sys.integrate()

The angles can be plotted to see how they change with respect to time given the initial conditions.

def plot():
    plt.figure()
    plt.plot(sys.times, np.rad2deg(x[:, :2]))
    plt.legend([sm.latex(s, mode='inline') for s in sys.coordinates])

plot()
../_images/chaos-pendulum_30_0.png

Chaotic Behavior

Now change the intial condition of the plat angle just slighty to see if the behvior of the system is similar.

sys.initial_conditions[phi] = np.deg2rad(1.0)
x = sys.integrate()
plot()
../_images/chaos-pendulum_31_0.png

Seems all good, very similar behavior. But now set the rod angle to \(90^\circ\) and try the same slight change in plate angle.

sys.initial_conditions[theta] = np.deg2rad(90)
sys.initial_conditions[phi] = np.deg2rad(0.5)
x = sys.integrate()
plot()
../_images/chaos-pendulum_32_0.png

First note that the plate behaves wildly. What happens when the initial plate angle is altered slightly.

sys.initial_conditions[phi] = np.deg2rad(1.0)
x = sys.integrate()
plot()
../_images/chaos-pendulum_33_0.png

The behavior does not look similar to the previous simulation. This is an example of chaotic behavior. The plate angle can not be reliably predicted because slight changes in the initial conditions cause the behavior of the system to vary widely.

Visualization

Finally, the system can be animated by attached a cylinder and a plane shape to the rigid bodies. To properly align the coordinate axes of the shapes with the bodies, simple rotations are used.

rod_shape = Cylinder(2 * lA, 0.005, color='red', name='rod')
plate_shape = Plane(w, h, color='blue', name='plate')

v1 = VisualizationFrame('rod',
                        A.orientnew('rod', 'Axis', (sm.pi / 2, A.x)),
                        Ao,
                        rod_shape)

v2 = VisualizationFrame('plate',
                        B.orientnew('plate', 'Body', (sm.pi / 2, sm.pi / 2, 0), 'XZX'),
                        Bo,
                        plate_shape)

scene = Scene(N, No, v1, v2, system=sys)

The following method opens up a simple gui that shows a 3D animatoin of the system.

scene.display_jupyter(axes_arrow_length=1.0)