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.
In [1]:
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
In [2]:
%matplotlib nbagg
In [3]:
me.init_vprinting(use_latex='mathjax')
In [4]:
mA, mB, lB, w, h, g = sm.symbols('m_A, m_B, L_B, w, h, g')
There are two time varying generalized coordinates:
The two generalized speeds will then be defined as:
In [5]:
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:
In [6]:
kin_diff = (omega - theta.diff(), alpha - phi.diff())
kin_diff
Out[6]:
In [7]:
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:
In [8]:
A.orient(N, 'Axis', (theta, N.y))
B.orient(A, 'Axis', (phi, A.z))
In [9]:
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$.
In [10]:
lA = (lB - h / 2) / 2
Ao.set_pos(No, lA * A.z)
Bo.set_pos(No, lB * A.z)
In [11]:
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.
In [12]:
No.set_vel(N, 0)
In [13]:
Ao.v2pt_theory(No, N, A)
Out[13]:
In [14]:
Bo.v2pt_theory(No, N, A)
Out[14]:
In [15]:
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:
In [16]:
IA[0].to_matrix(A)
Out[16]:
The central inerita of the symmetric plate with respect to its reference frame is a function of its width and height.
In [17]:
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)
In [18]:
IB[0].to_matrix(B)
Out[18]:
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.
In [19]:
rod = me.RigidBody('rod', Ao, A, mA, IA)
In [20]:
plate = me.RigidBody('plate', Bo, B, mB, IB)
In [21]:
rod_gravity = (Ao, mA * g * N.z)
plate_gravity = (Bo, mB * g * N.z)
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.
In [22]:
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^*$.
In [23]:
bodies = (rod, plate)
loads = (rod_gravity, plate_gravity)
fr, frstar = kane.kanes_equations(loads, bodies)
In [24]:
sm.trigsimp(fr)
Out[24]:
In [25]:
sm.trigsimp(frstar)
Out[25]:
In [26]:
sys = System(kane)
In [27]:
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
In [28]:
sys.initial_conditions = {theta: np.deg2rad(45),
phi: np.deg2rad(0.5),
omega: 0,
alpha: 0}
In [29]:
sys.times = np.linspace(0, 10, 500)
The trajectories of the states are found with the integrate
method.
In [30]:
x = sys.integrate()
The angles can be plotted to see how they change with respect to time given the initial conditions.
In [31]:
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()
In [32]:
sys.initial_conditions[phi] = np.deg2rad(1.0)
x = sys.integrate()
plot()
Seems all good, very similar behavior. But now set the rod angle to $90^\circ$ and try the same slight change in plate angle.
In [33]:
sys.initial_conditions[theta] = np.deg2rad(90)
sys.initial_conditions[phi] = np.deg2rad(0.5)
x = sys.integrate()
plot()
First note that the plate behaves wildly. What happens when the initial plate angle is altered slightly.
In [34]:
sys.initial_conditions[phi] = np.deg2rad(1.0)
x = sys.integrate()
plot()
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.
In [35]:
rod_shape = Cylinder(2 * lA, 0.005, color='red')
plate_shape = Plane(h, w, color='blue')
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.
In [36]:
scene.display_ipython()
In [37]:
%load_ext version_information
%version_information numpy, sympy, scipy, matplotlib, pydy
Out[37]: