```
In [1]:
```from IPython.display import SVG
SVG(filename='mass_spring_damper.svg')

```
Out[1]:
```

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

```
In [2]:
```import sympy as sym
import sympy.physics.mechanics as me

```
In [3]:
```from sympy.physics.vector import init_vprinting
init_vprinting()

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

```
In [4]:
```x, v = me.dynamicsymbols('x v')

```
In [5]:
```m, c, k, g, t = sym.symbols('m c k g t')

```
In [6]:
```ceiling = me.ReferenceFrame('C')

```
In [7]:
```O = me.Point('O')
P = me.Point('P')

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

```
In [8]:
```O.set_vel(ceiling, 0)

```
In [9]:
```P.set_pos(O, x * ceiling.x)
P.set_vel(ceiling, v * ceiling.x)
P.vel(ceiling)

```
Out[9]:
```

```
In [10]:
```damping = -c * P.vel(ceiling)
stiffness = -k * P.pos_from(O)
gravity = m * g * ceiling.x
forces = damping + stiffness + gravity
forces

```
Out[10]:
```

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

```
In [11]:
```zero = me.dot(forces - m * P.acc(ceiling), ceiling.x)
zero

```
Out[11]:
```

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

```
Out[12]:
```

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

```
In [13]:
```mass = me.Particle('mass', P, m)

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

```
In [14]:
```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^*$.

```
In [15]:
```fr, frstar = kane.kanes_equations([(P, forces)], [mass])
fr, frstar

```
Out[15]:
```

```
In [16]:
```M = kane.mass_matrix_full
f = kane.forcing_full
M, f

```
Out[16]:
```

```
In [17]:
```M.inv() * f

```
Out[17]:
```

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.

```
In [18]:
```from pydy.system import System

```
In [19]:
```sys = System(kane)

```
```

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

`linspace`

is often useful for this.

```
In [21]:
```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.

```
In [22]:
```x_trajectory = sys.integrate()

PyDy has a native module `pydy.viz`

which is used to visualize a System in an interactive 3D GUI.

```
In [23]:
```from pydy.viz import *

```
```

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

```
In [24]:
```bob = Sphere(2.0, color="red", material="metal")
bob_vframe = VisualizationFrame(ceiling, P, bob)

```
In [25]:
```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)

`System`

onto a canvas.
It takes a ReferenceFrame and Point as arguments.

```
In [26]:
```scene = Scene(ceiling, O, system=sys)

```
```

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

```
In [27]:
```scene.visualization_frames = [bob_vframe, ceiling_vframe]

```
In [28]:
```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.

```
In [29]:
```scene.display_ipython()