# Simulation of Ball drop and Spring mass damper system
"Simulation of dynamic systems for dummies".


This is a very simple description of how to do time simulations of a dynamic system using SciPy ODE (Ordinary Differnetial Equation) Solver.


```python
from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt
```

## Simulation of a static system to introduce ODEint

Define a method that takes a system state and describe how this state will change in time. The method does this by returning time derivatives for each state. The ODE solver will use these time derivatives to calculate new states, for the next time step.

Here is a method that takes a system to simulate a train that travels with constant speed:

(The system has only one state, the position of the train)


```python
V_start = 150*10**3/3600  # [m/s] Train velocity at start

def train(states,t):
    
    # states:
    # [x]
    
    x = states[0]  # Position of train
    dxdt = V_start # The position state will change by the speed of the train
    
    # Time derivative of the states:
    d_states_dt = np.array([dxdt])
    
    return d_states_dt
```


```python
x_start = 0  # [m] Train position at start

# The states at start of the simulation, the train is traveling with constant speed V at position x = 0.
states_0 = np.array([x_start])

# Create a time vector for the simulation:
t = np.linspace(0,10,100)

# Simulate with the "train" method and start states for the times in t:
states = odeint(func = train,y0 = states_0,t = t)

# The result is the time series of the states:
x = states[:,0]

```


```python
fig,ax = plt.subplots()
ax.plot(t,x,label = 'Train position')

ax.set_title('Train traveling at constant speed')
ax.set_xlabel('time [s]')
ax.set_ylabel('x [m]')
a = ax.legend()

```

The speed can hower be a state too:


```python
def train_2_states(states,t):
    
    # states:
    # [x,V]
    
    x = states[0]  # Position of train
    V = states[1]  # Speed of train
    dxdt = V  # The position state will change by the speed of the train
    dVdt = 0  # The velocity will not change (No acceleration)
    
    # Time derivative of the states:
    d_states_dt = np.array([dxdt,dVdt])
    
    return d_states_dt
```


```python
# The states at start of the simulation, the train is traveling with constant speed V at position x = 0.
states_0 = np.array([x_start,V_start])

# Create a time vector for the simulation:
t = np.linspace(0,10,100)

# Simulate with the "train" method and start states for the times in t:
states = odeint(func = train_2_states,y0 = states_0,t = t)

# The result is the time series of the states:
x = states[:,0]
dxdt = states[:,1]

```


```python
fig,axes = plt.subplots(ncols = 2)
fig.set_size_inches(11,5)
ax = axes[0]
ax.plot(t,x,label = 'Train position')

ax.set_title('Train traveling at constant speed')
ax.set_xlabel('time [s]')
ax.set_ylabel('x [m]')
a = ax.legend()

ax = axes[1]
ax.plot(t,dxdt,label = 'Train speed')

ax.set_title('Train traveling at constant speed')
ax.set_xlabel('time [s]')
ax.set_ylabel('dx/dt [m/s]')
a = ax.legend()

```

## Ball drop
Here is a system where the speed is not constant.
A simulation of a ball drop under the influence of gravity force.


```python
g = 9.81
m = 1

def ball_drop(states,t):
    
    # states:
    # [x,v]
    # F = g*m = m*dv/dt
    # --> dv/dt = (g*m) / m
    
    x = states[0]
    dxdt = states[1]
    
    dvdt = (g*m) / m
    
    d_states_dt = np.array([dxdt,dvdt])
    
    return d_states_dt
```


```python
states_0 = np.array([0,0])
t = np.linspace(0,10,100)
states = odeint(func = ball_drop,y0 = states_0,t = t)
x = states[:,0]
dxdt = states[:,1]

fig,axes = plt.subplots(ncols = 2)
fig.set_size_inches(11,5)
ax = axes[0]
ax.plot(t,x,label = 'Ball position')

ax.set_title('Ball drop')
ax.set_xlabel('time [s]')
ax.set_ylabel('x [m]')
a = ax.legend()

ax = axes[1]
ax.plot(t,dxdt,label = 'Ball speed')

ax.set_title('Ball drop')
ax.set_xlabel('time [s]')
ax.set_ylabel('dx/dt [m/s]')
a = ax.legend()
```

Simulating in air, where the ball has a resistance due aerodynamic drag.


```python
cd = 0.01

def ball_drop_air(states,t):
    
    # states:
    # [x,u]
    # F = g*m - cd*u = m*du/dt
    # --> du/dt = (g*m - cd*u**2) / m
    
    x = states[0]
    u = states[1]
    dxdt = u
    
    dudt = (g*m - cd*u**2) / m
    
    d_states_dt = np.array([dxdt,dudt])
    
    return d_states_dt
```


```python
states = odeint(func = ball_drop_air,y0 = states_0,t = t)
x_air = states[:,0]
dxdt_air = states[:,1]

fig,axes = plt.subplots(ncols = 2)
fig.set_size_inches(11,5)
ax = axes[0]
ax.plot(t,x,label = 'Vacuum')
ax.plot(t,x_air,label = 'Air')

ax.set_title('Ball drop in vacuum and air')
ax.set_xlabel('time [s]')
ax.set_ylabel('x [m]')
a = ax.legend()

ax = axes[1]
ax.plot(t,dxdt,label = 'Vacuum')
ax.plot(t,dxdt_air,label = 'Air')

ax.set_title('Ball drop in vacuum and air')
ax.set_xlabel('time [s]')
ax.set_ylabel('dx/dt [m/s]')
a = ax.legend()

```

The very classical dynamic system with a spring, a mass and a damper.



```python
k = 3  # The stiffnes of the spring (relates to position)
c = 0.1  # Damping term (relates to velocity)
m = 0.1  # The mass (relates to acceleration)

def spring_mass_damp(states,t):
    
    # states:
    # [x,v]
    # F = -k*x -c*v = m*dv/dt
    # --> dv/dt = (-kx -c*v) / m
    
    x = states[0]
    dxdt = states[1]
    
    dvdt = (-k*x -c*dxdt) / m
    
    d_states_dt = np.array([dxdt,dvdt])
    
    return d_states_dt
```


```python
y0 = np.array([1,0])
t = np.linspace(0,10,100)
states = odeint(func = spring_mass_damp,y0 = y0,t = t)
x = states[:,0]
dxdt = states[:,1]
```


```python
fig,ax = plt.subplots()
ax.plot(t,x)
ax.set_title('Spring mass damper simulation')
ax.set_xlabel('time [s]')
a = ax.set_ylabel('x [m]')

```

Also add a gravity force


```python
g = 9.81

def spring_mass_damp_g(states,t):
    
    # states:
    # [x,v]
    # F = g*m -k*x -c*v = m*dv/dt
    # --> dv/dt = (g*m -kx -c*v) / m
    
    x = states[0]
    dxdt = states[1]
    
    dvdt = (g*m -k*x -c*dxdt) / m
    
    d_states_dt = np.array([dxdt,dvdt])
    
    return d_states_dt
```


```python
states_g = odeint(func = spring_mass_damp_g,y0 = y0,t = t)
x_g = states_g[:,0]
dxdt_g = states_g[:,1]
```


```python
fig,ax = plt.subplots()
ax.plot(t,x,label = 'No gravity force')
ax.plot(t,x_g,label = 'Gravity force')

ax.set_title('Spring mass damper simulation with and without gravity')
ax.set_xlabel('time [s]')
ax.set_ylabel('x [m]')
a = ax.legend()
```

## SymPy solution


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


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


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


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


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


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


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


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




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




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


```python

```


```python

```

19536918fc8d7365498ac0e6afd51f2dd72b08e5
ipynb
Jupyter Notebook
notebooks/01_dynamics_experiments/spring_mass_damper.ipynb
martinlarsalbert/ForceMan
b07541f1c04d0ed103d7a5f66256812ebca2824c
MIT
notebooks/01_dynamics_experiments/spring_mass_damper.ipynb
martinlarsalbert/ForceMan
b07541f1c04d0ed103d7a5f66256812ebca2824c
MIT
2019-12-15T17:23:21.000Z
2019-12-15T17:23:21.000Z
notebooks/01_dynamics_experiments/spring_mass_damper.ipynb
martinlarsalbert/ForceMan
b07541f1c04d0ed103d7a5f66256812ebca2824c
MIT
Qwen/Qwen-72B
1. YES
2. YES