# pydy

[pydy](https://github.com/pydy/pydy) simulates physical systems.

Example from [pydy/../mass_spring_damper](https://github.com/pydy/pydy/tree/master/examples/mass_spring_damper).

[documentation](http://pydy.readthedocs.io/en/stable/index.html) [examples](https://github.com/pydy/pydy/tree/master/examples) [sources](https://github.com/pydy/pydy)

Here we will derive the equations of motion for the classic mass-spring-damper system under the influence of gravity. The following figure gives a pictorial description of the problem.


```python
from IPython.display import SVG
SVG(filename='pydy.svg')
```




    

    



Start by loading in the core functionality of both SymPy and Mechanics.


```python
import sympy as sym
import sympy.physics.mechanics as me
```

We can make use of the pretty printing of our results by loading SymPy's printing extension, in particular we will use the vector printing which is nice for mechanics objects.


```python
from sympy.physics.vector import init_vprinting
init_vprinting(use_latex='mathjax')
```

We'll start by defining the variables we will need for this problem:
- $x(t)$: distance of the particle from the ceiling
- $v(t)$: speed of the particle
- $m$: mass of the particle
- $c$: damping coefficient of the damper
- $k$: stiffness of the spring
- $g$: acceleration due to gravity
- $t$: time


```python
x, v = me.dynamicsymbols('x v')
```


```python
m, c, k, g, t = sym.symbols('m c k g t')
```

Now, we define a Newtonian reference frame that represents the ceiling which the particle is attached to, $C$.


```python
ceiling = me.ReferenceFrame('C')
```

We will need two points, one to represent the original position of the particle which stays fixed in the ceiling frame, $O$, and the second one, $P$ which is aligned with the particle as it moves.


```python
O = me.Point('O')
P = me.Point('P')
```

The velocity of point $O$ in the ceiling is zero.


```python
O.set_vel(ceiling, 0)
```

Point $P$ can move downward in the $y$ direction and its velocity is specified as $v$ in the downward direction.


```python
P.set_pos(O, x * ceiling.x)
P.set_vel(ceiling, v * ceiling.x)
P.vel(ceiling)
```




$$v\mathbf{\hat{c}_x}$$



There are three forces acting on the particle. Those due to the acceleration of gravity, the damper, and the spring.


```python
damping = -c * P.vel(ceiling)
stiffness = -k * P.pos_from(O)
gravity = m * g * ceiling.x
forces = damping + stiffness + gravity
forces
```




$$(- c v + g m - k x)\mathbf{\hat{c}_x}$$



Now we can use Newton's second law, $0=F-ma$, to form the equation of motion of the system.


```python
zero = me.dot(forces - m * P.acc(ceiling), ceiling.x)
zero
```




$$- c v + g m - k x - m \dot{v}$$



We can then form the first order equations of motion by solving for $\frac{dv}{dt}$ and introducing the kinematical differential equation, $v=\frac{dx}{dt}$.


```python
dv_by_dt = sym.solve(zero, v.diff(t))[0]
dx_by_dt = v
dv_by_dt, dx_by_dt
```




$$\left ( \frac{1}{m} \left(- c v + g m - k x\right), \quad v\right )$$



Forming the equations of motion can also be done with the automated methods available in the Mechanics package: `LagrangesMethod` and `KanesMethod`. Here we will make use of Kane's method to find the same equations of motion that we found manually above. First, define a particle that represents the mass attached to the damper and spring.


```python
mass = me.Particle('mass', P, m)
```

Now we can construct a `KanesMethod` object by passing in the generalized coordinate, $x$, the generalized speed, $v$, and the kinematical differential equation which relates the two, $0=v-\frac{dx}{dt}$.


```python
kane = me.KanesMethod(ceiling, q_ind=[x], u_ind=[v], kd_eqs=[v - x.diff(t)])
```

Now Kane's equations can be computed, and we can obtain $F_r$ and $F_r^*$.


```python
fr, frstar = kane.kanes_equations([(P, forces)], [mass])
fr, frstar
```




$$\left ( \left[\begin{matrix}- c v + g m - k x\end{matrix}\right], \quad \left[\begin{matrix}- m \dot{v}\end{matrix}\right]\right )$$



The equations are also available in the form $M\frac{d}{dt}[q,u]^T=f(q, u)$ and we can extract the mass matrix, $M$, and the forcing functions, $f$.


```python
M = kane.mass_matrix_full
f = kane.forcing_full
M, f
```




$$\left ( \left[\begin{matrix}1 & 0\\0 & m\end{matrix}\right], \quad \left[\begin{matrix}v\\- c v + g m - k x\end{matrix}\right]\right )$$



Finally, we can form the first order differential equations of motion $\frac{d}{dt}[q,u]^T=M^{-1}f(\dot{u}, u, q)$, which is the same as previously found.


```python
M.inv() * f
```




$$\left[\begin{matrix}v\\\frac{1}{m} \left(- c v + g m - k x\right)\end{matrix}\right]$$



# Simulating the system

Now that we have defined the mass-spring-damper system, we are going to simulate it.

PyDy's `System` is a wrapper that holds the Kanes object to integrate the equations of motion using numerical values of constants.


```python
from pydy.system import System
```


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

Now, we specify the numerical values of the constants and the initial values of states in the form of a dict.


```python
sys.constants = {m:10.0, g:9.8, c:5.0, k:10.0}
sys.initial_conditions = {x:0.0, v:0.0}
```

We must generate a time vector over which the integration will be carried out. NumPy's `linspace` is often useful for this.


```python
from numpy import linspace
sys.times = linspace(0.0, 10.0, 100)
```

The trajectory of the states over time can be found by calling the `.integrate()` method.


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

# Visualizing the System

PyDy has a native module `pydy.viz` which is used to visualize a System in an interactive 3D GUI.


```python
from pydy.viz import *
```

For visualizing the system, we need to create shapes for the objects we wish to visualize, and map each of them 
to a `VisualizationFrame`, which holds the position and orientation of the object. First create a sphere to represent the bob and attach it to the point $P$ and the ceiling reference frame (the sphere does not rotate with respect to the ceiling).


```python
bob = Sphere(2.0, color="red", material="metal")
bob_vframe = VisualizationFrame(ceiling, P, bob)
```

Now create a circular disc that represents the ceiling and fix it to the ceiling reference frame. The circle's default axis is aligned with its local $y$ axis, so we need to attach it to a rotated ceiling reference frame if we want the circle's axis to align with the $\hat{c}_x$ unit vector.


```python
ceiling_circle = Circle(radius=10, color="white", material="metal")
from numpy import pi
rotated = ceiling.orientnew("C_R", 'Axis', [pi / 2, ceiling.z])
ceiling_vframe = VisualizationFrame(rotated, O, ceiling_circle)
```

Now we initialize a Scene. A Scene contains all the information required to visualize a `System` onto a canvas.
It takes a ReferenceFrame and Point as arguments.


```python
scene = Scene(ceiling, O, system=sys)
```

We provide the VisualizationFrames, which we want to visualize as a list to scene.


```python
scene.visualization_frames = [bob_vframe, ceiling_vframe]
```

The default camera of Scene has the z axis of the base frame pointing out of the screen, and the y axis pointing up. We want the x axis to point downwards, so we supply a new camera that will achieve this.


```python
camera_frame = ceiling.orientnew('Camera Frame','Axis', [pi / 2, ceiling.z])
camera_point = O.locatenew('Camera Location', 100 * camera_frame.z)
primary_camera = PerspectiveCamera(camera_frame, camera_point)
scene.cameras = [primary_camera]
```

Now, we call the display method.


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


```python
%load_ext version_information
```


```python
%version_information pydy, numpy, scipy
```




<table><tr><th>Software</th><th>Version</th></tr><tr><td>Python</td><td>3.5.0 64bit [MSC v.1900 64 bit (AMD64)]</td></tr><tr><td>IPython</td><td>4.2.0</td></tr><tr><td>OS</td><td>Windows 10.0.10586</td></tr><tr><td>pydy</td><td>0.3.1</td></tr><tr><td>numpy</td><td>1.11.1rc1</td></tr><tr><td>scipy</td><td>0.17.1</td></tr><tr><td colspan='2'>Tue Jun 14 17:29:07 2016 Paris, Madrid (heure d?été)</td></tr></table>




```python

```

791a4936d0bddd39d53932a56ddb01610234bd88
ipynb
Jupyter Notebook
_doc/notebooks/2016/pydata/js_pydy_mass_spring_damper.ipynb
sdpython/jupytalk
34abdf128de24becb21a9f08f243c3a74dadbfd9
MIT
_doc/notebooks/2016/pydata/js_pydy_mass_spring_damper.ipynb
sdpython/jupytalk
34abdf128de24becb21a9f08f243c3a74dadbfd9
MIT
2016-11-13T19:52:35.000Z
2021-12-29T10:59:41.000Z
_doc/notebooks/2016/pydata/js_pydy_mass_spring_damper.ipynb
sdpython/jupytalk
34abdf128de24becb21a9f08f243c3a74dadbfd9
MIT
2016-09-10T10:44:50.000Z
2021-09-22T16:28:56.000Z
Qwen/Qwen-72B
1. YES
2. YES