# Introduction

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.



# Setup


```python
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
```


```python
%matplotlib nbagg
```


```python
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.


```python
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.


```python
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:

$$0 = \omega - \dot{\theta}\\
0 = \alpha - \dot{\phi}$$


```python
kin_diff = (omega - theta.diff(), alpha - phi.diff())
kin_diff
```




$$\left ( \omega - \dot{\theta}, \quad \alpha - \dot{\phi}\right )$$



# Define Orientations

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


```python
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:


```python
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.


```python
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$.


```python
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:


```python
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.


```python
No.set_vel(N, 0)
```


```python
Ao.v2pt_theory(No, N, A)
```




$$\left(\frac{L_{B}}{2} - \frac{h}{4}\right) \omega\mathbf{\hat{a}_x}$$




```python
Bo.v2pt_theory(No, N, A)
```




$$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.


```python
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:


```python
IA[0].to_matrix(A)
```




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



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


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


```python
IB[0].to_matrix(B)
```




$$\left[\begin{matrix}\frac{h^{2} m_{B}}{12} & 0 & 0\\0 & \frac{m_{B}}{12} \left(h^{2} + w^{2}\right) & 0\\0 & 0 & \frac{m_{B} w^{2}}{12}\end{matrix}\right]$$



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.


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


```python
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.


```python
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.


```python
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^*$.


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

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


```python
sm.trigsimp(fr)
```




$$\left[\begin{matrix}g \left(- \frac{L_{B} m_{A}}{2} - L_{B} m_{B} + \frac{h m_{A}}{4}\right) \operatorname{sin}\left(\theta\right)\\0\end{matrix}\right]$$




```python
sm.trigsimp(frstar)
```




$$\left[\begin{matrix}\frac{m_{B} w^{2}}{12} \alpha \omega \operatorname{sin}\left(2 \phi\right) - \left(\frac{L_{B}^{2} m_{A}}{3} + L_{B}^{2} m_{B} - \frac{L_{B} h}{3} m_{A} + \frac{h^{2} m_{A}}{12} + \frac{h^{2} m_{B}}{12} + \frac{m_{B} w^{2}}{12} \operatorname{cos}^{2}\left(\phi\right)\right) \dot{\omega}\\- \frac{m_{B} w^{2}}{24} \left(\omega^{2} \operatorname{sin}\left(2 \phi\right) + 2 \dot{\alpha}\right)\end{matrix}\right]$$



# 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.


```python
sys = System(kane)
```


```python
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
```


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


```python
sys.times = np.linspace(0, 10, 500)
```

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


```python
x = sys.integrate()
```

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


```python
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()
```


    <IPython.core.display.Javascript object>






# Chaotic Behavior

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


```python
sys.initial_conditions[phi] = np.deg2rad(1.0)
x = sys.integrate()
plot()  
```


    <IPython.core.display.Javascript object>






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


```python
sys.initial_conditions[theta] = np.deg2rad(90)
sys.initial_conditions[phi] = np.deg2rad(0.5)
x = sys.integrate()
plot()
```


    <IPython.core.display.Javascript object>






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


```python
sys.initial_conditions[phi] = np.deg2rad(1.0)
x = sys.integrate()
plot()
```


    <IPython.core.display.Javascript object>






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.


```python
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.


```python
scene.display_ipython()
```


```python
%load_ext version_information
%version_information numpy, sympy, scipy, matplotlib, pydy
```




<table><tr><th>Software</th><th>Version</th></tr><tr><td>Python</td><td>3.5.3 64bit [GCC 4.4.7 20120313 (Red Hat 4.4.7-1)]</td></tr><tr><td>IPython</td><td>4.2.0</td></tr><tr><td>OS</td><td>Linux 4.10.0 21 generic x86_64 with debian stretch sid</td></tr><tr><td>numpy</td><td>1.12.1</td></tr><tr><td>sympy</td><td>1.0</td></tr><tr><td>scipy</td><td>0.19.0</td></tr><tr><td>matplotlib</td><td>2.0.2</td></tr><tr><td>pydy</td><td>0.4.0dev</td></tr><tr><td colspan='2'>Mon May 29 13:08:51 2017 PDT</td></tr></table>



3ce27df1496ba02a2e7a4b713c81e51022f81033
ipynb
Jupyter Notebook
examples/chaos_pendulum/chaos_pendulum.ipynb
JanMV/pydy
22c6a3965853bb3641c63e493717976d775034e2
BSD-3-Clause
2015-01-31T11:43:22.000Z
2022-03-15T02:18:21.000Z
examples/chaos_pendulum/chaos_pendulum.ipynb
JanMV/pydy
22c6a3965853bb3641c63e493717976d775034e2
BSD-3-Clause
2015-01-17T16:56:42.000Z
2022-02-08T05:27:08.000Z
examples/chaos_pendulum/chaos_pendulum.ipynb
JanMV/pydy
22c6a3965853bb3641c63e493717976d775034e2
BSD-3-Clause
2015-02-03T13:02:45.000Z
2021-12-21T12:57:21.000Z
Qwen/Qwen-72B
1. YES
2. YES