.. _jspydymassspringdamperrst:

====
pydy
====


.. only:: html

    **Links:** :download:`notebook <js_pydy_mass_spring_damper.ipynb>`, :downloadlink:`html <js_pydy_mass_spring_damper2html.html>`, :download:`PDF <js_pydy_mass_spring_damper.pdf>`, :download:`python <js_pydy_mass_spring_damper.py>`, :downloadlink:`slides <js_pydy_mass_spring_damper.slides.html>`, :githublink:`GitHub|_doc/notebooks/2016/pydata/js_pydy_mass_spring_damper.ipynb|*`


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

.. code:: ipython3

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




.. image:: js_pydy_mass_spring_damper_2_0.svg



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

.. code:: ipython3

    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.

.. code:: ipython3

    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: -
:math:`x(t)`: distance of the particle from the ceiling - :math:`v(t)`:
speed of the particle - :math:`m`: mass of the particle - :math:`c`:
damping coefficient of the damper - :math:`k`: stiffness of the spring -
:math:`g`: acceleration due to gravity - :math:`t`: time

.. code:: ipython3

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

.. code:: ipython3

    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, :math:`C`.

.. code:: ipython3

    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, :math:`O`, and the
second one, :math:`P` which is aligned with the particle as it moves.

.. code:: ipython3

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

The velocity of point :math:`O` in the ceiling is zero.

.. code:: ipython3

    O.set_vel(ceiling, 0)

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

.. code:: ipython3

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




.. math::

    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.

.. code:: ipython3

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




.. math::

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



Now we can use Newton’s second law, :math:`0=F-ma`, to form the equation
of motion of the system.

.. code:: ipython3

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




.. math::

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



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

.. code:: ipython3

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




.. math::

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

.. code:: ipython3

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

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

.. code:: ipython3

    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 :math:`F_r` and
:math:`F_r^*`.

.. code:: ipython3

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




.. math::

    \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
:math:`M\frac{d}{dt}[q,u]^T=f(q, u)` and we can extract the mass matrix,
:math:`M`, and the forcing functions, :math:`f`.

.. code:: ipython3

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




.. math::

    \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
:math:`\frac{d}{dt}[q,u]^T=M^{-1}f(\dot{u}, u, q)`, which is the same as
previously found.

.. code:: ipython3

    M.inv() * f




.. math::

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

.. code:: ipython3

    from pydy.system import System

.. code:: ipython3

    sys = System(kane)

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

.. code:: ipython3

    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.

.. code:: ipython3

    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.

.. code:: ipython3

    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.

.. code:: ipython3

    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 :math:`P` and the
ceiling reference frame (the sphere does not rotate with respect to the
ceiling).

.. code:: ipython3

    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 :math:`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
:math:`\hat{c}_x` unit vector.

.. code:: ipython3

    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.

.. code:: ipython3

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

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

.. code:: ipython3

    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.

.. code:: ipython3

    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.

.. code:: ipython3

    scene.display_ipython()

.. code:: ipython3

    %load_ext version_information

.. code:: ipython3

    %version_information pydy, numpy, scipy






.. raw:: html

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