# Modeling and Simulation in Python

Case study

Copyright 2017 Allen Downey

License: [Creative Commons Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0)



```python
# Configure Jupyter so figures appear in the notebook
%matplotlib inline

# Configure Jupyter to display the assigned value after an assignment
%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'

# import functions from the modsim.py module
from modsim import *
```

## Yo-yo

Suppose you are holding a yo-yo with a length of string wound around its axle, and you drop it while holding the end of the string stationary.  As gravity accelerates the yo-yo downward, tension in the string exerts a force upward.  Since this force acts on a point offset from the center of mass, it exerts a torque that causes the yo-yo to spin.



This figure shows the forces on the yo-yo and the resulting torque.  The outer shaded area shows the body of the yo-yo.  The inner shaded area shows the rolled up string, the radius of which changes as the yo-yo unrolls.

In this model, we can't figure out the linear and angular acceleration independently; we have to solve a system of equations:

$\sum F = m a $

$\sum \tau = I \alpha$

where the summations indicate that we are adding up forces and torques.

As in the previous examples, linear and angular velocity are related because of the way the string unrolls:

$\frac{dy}{dt} = -r \frac{d \theta}{dt} $

In this example, the linear and angular accelerations have opposite sign.  As the yo-yo rotates counter-clockwise, $\theta$ increases and $y$, which is the length of the rolled part of the string, decreases.

Taking the derivative of both sides yields a similar relationship between linear and angular acceleration:

$\frac{d^2 y}{dt^2} = -r \frac{d^2 \theta}{dt^2} $

Which we can write more concisely:

$ a = -r \alpha $

This relationship is not a general law of nature; it is specific to scenarios like this where there is rolling without stretching or slipping.

Because of the way we've set up the problem, $y$ actually has two meanings: it represents the length of the rolled string and the height of the yo-yo, which decreases as the yo-yo falls.  Similarly, $a$ represents acceleration in the length of the rolled string and the height of the yo-yo.

We can compute the acceleration of the yo-yo by adding up the linear forces:

$\sum F = T - mg = ma $

Where $T$ is positive because the tension force points up, and $mg$ is negative because gravity points down.

Because gravity acts on the center of mass, it creates no torque, so the only torque is due to tension:

$\sum \tau = T r = I \alpha $

Positive (upward) tension yields positive (counter-clockwise) angular acceleration.

Now we have three equations in three unknowns, $T$, $a$, and $\alpha$, with $I$, $m$, $g$, and $r$ as known parameters.  It is simple enough to solve these equations by hand, but we can also get SymPy to do it for us.




```python
from sympy import init_printing, symbols, Eq, solve

init_printing()
```


```python
T, a, alpha, I, m, g, r = symbols('T a alpha I m g r')
```


```python
eq1 = Eq(a, -r * alpha)
```


```python
eq2 = Eq(T - m * g, m * a)
```


```python
eq3 = Eq(T * r, I * alpha)
```


```python
soln = solve([eq1, eq2, eq3], [T, a, alpha])
```


```python
soln[T]
```


```python
soln[a]
```


```python
soln[alpha]
```


The results are

$T      = m g I / I^*   $

$a      = -m g r^2 / I^* $

$\alpha = m g r / I^*    $

where $I^*$ is the augmented moment of inertia, $I + m r^2$.

You can also see [the derivation of these equations in this video](https://www.youtube.com/watch?v=chC7xVDKl4Q).

To simulate the system, we don't really need $T$; we can plug $a$ and $\alpha$ directly into the slope function.


```python
radian = UNITS.radian
m = UNITS.meter
s = UNITS.second
kg = UNITS.kilogram
N = UNITS.newton
```




newton



**Exercise:**  Simulate the descent of a yo-yo.  How long does it take to reach the end of the string?

I provide a `Params` object with the system parameters:

* `Rmin` is the radius of the axle.  `Rmax` is the radius of the axle plus rolled string.

* `Rout` is the radius of the yo-yo body.  `mass` is the total mass of the yo-yo, ignoring the string.  

* `L` is the length of the string.

* `g` is the acceleration of gravity.


```python
params = Params(Rmin = 8e-3 * m,
                Rmax = 16e-3 * m,
                Rout = 35e-3 * m,
                mass = 50e-3 * kg,
                L = 1 * m,
                g = 9.8 * m / s**2,
                t_end = 1 * s)
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>values</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>Rmin</th>
      <td>0.008 meter</td>
    </tr>
    <tr>
      <th>Rmax</th>
      <td>0.016 meter</td>
    </tr>
    <tr>
      <th>Rout</th>
      <td>0.035 meter</td>
    </tr>
    <tr>
      <th>mass</th>
      <td>0.05 kilogram</td>
    </tr>
    <tr>
      <th>L</th>
      <td>1 meter</td>
    </tr>
    <tr>
      <th>g</th>
      <td>9.8 meter / second ** 2</td>
    </tr>
    <tr>
      <th>t_end</th>
      <td>1 second</td>
    </tr>
  </tbody>
</table>
</div>



Here's a `make_system` function that computes `I` and `k` based on the system parameters.

I estimated `I` by modeling the yo-yo as a solid cylinder with uniform density ([see here](https://en.wikipedia.org/wiki/List_of_moments_of_inertia)).

In reality, the distribution of weight in a yo-yo is often designed to achieve desired effects.  But we'll keep it simple.


```python
def make_system(params):
    """Make a system object.
    
    params: Params with Rmin, Rmax, Rout, 
                              mass, L, g, t_end
    
    returns: System with init, k, Rmin, Rmax, mass,
                         I, g, ts
    """
    L, mass = params.L, params.mass
    Rout, Rmax, Rmin = params.Rout, params.Rmax, params.Rmin 
    
    init = State(theta = 0 * radian,
                 omega = 0 * radian/s,
                 y = L,
                 v = 0 * m / s)
    
    I = mass * Rout**2 / 2
    k = (Rmax**2 - Rmin**2) / 2 / L / radian    
    
    return System(params, init=init, I=I, k=k)
```

Testing `make_system`


```python
system = make_system(params)
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>values</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>Rmin</th>
      <td>0.008 meter</td>
    </tr>
    <tr>
      <th>Rmax</th>
      <td>0.016 meter</td>
    </tr>
    <tr>
      <th>Rout</th>
      <td>0.035 meter</td>
    </tr>
    <tr>
      <th>mass</th>
      <td>0.05 kilogram</td>
    </tr>
    <tr>
      <th>L</th>
      <td>1 meter</td>
    </tr>
    <tr>
      <th>g</th>
      <td>9.8 meter / second ** 2</td>
    </tr>
    <tr>
      <th>t_end</th>
      <td>1 second</td>
    </tr>
    <tr>
      <th>init</th>
      <td>theta               0 radian
omega    0.0 radi...</td>
    </tr>
    <tr>
      <th>I</th>
      <td>3.0625000000000006e-05 kilogram * meter ** 2</td>
    </tr>
    <tr>
      <th>k</th>
      <td>9.6e-05 meter / radian</td>
    </tr>
  </tbody>
</table>
</div>




```python
system.init
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>values</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>theta</th>
      <td>0 radian</td>
    </tr>
    <tr>
      <th>omega</th>
      <td>0.0 radian / second</td>
    </tr>
    <tr>
      <th>y</th>
      <td>1 meter</td>
    </tr>
    <tr>
      <th>v</th>
      <td>0.0 meter / second</td>
    </tr>
  </tbody>
</table>
</div>



Write a slope function for this system, using these results from the book:

$ r = \sqrt{2 k y + R_{min}^2} $ 

$ T      = m g I / I^*  $

$ a      = -m g r^2 / I^* $

$ \alpha  = m g r / I^*  $

where $I^*$ is the augmented moment of inertia, $I + m r^2$.



```python
# Solution

def slope_func(state, t, system):
    """Computes the derivatives of the state variables.
    
    state: State object with theta, omega, y, v
    t: time
    system: System object with Rmin, k, I, mass
    
    returns: sequence of derivatives
    """
    theta, omega, y, v = state
    g, k, Rmin = system.g, system.k, system.Rmin
    I, mass = system.I, system.mass
    
    r = sqrt(2*k*y + Rmin**2)
    alpha = mass * g * r / (I + mass * r**2)
    a = -r * alpha
        
    return omega, alpha, v, a        
```

Test your slope function with the initial paramss.


```python
# Solution

slope_func(system.init, 0*s, system)
```




    (0.0 <Unit('radian / second')>,
     180.54116292458264 <Unit('1 / radian ** 0.5 / second ** 2')>,
     0.0 <Unit('meter / second')>,
     -2.888658606793322 <Unit('meter / radian / second ** 2')>)



Write an event function that will stop the simulation when `y` is 0.


```python
# Solution

def event_func(state, t, system):
    """Stops when y is 0.
    
    state: State object with theta, omega, y, v
    t: time
    system: System object with Rmin, k, I, mass
    
    returns: y
    """
    theta, omega, y, v = state
    return y
```

Test your event function:


```python
# Solution

event_func(system.init, 0*s, system)
```




1 meter



Then run the simulation.


```python
# Solution

results, details = run_ode_solver(system, slope_func, events=event_func, max_step=0.05)
details
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>values</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>sol</th>
      <td>None</td>
    </tr>
    <tr>
      <th>t_events</th>
      <td>[[0.879217870162702]]</td>
    </tr>
    <tr>
      <th>nfev</th>
      <td>134</td>
    </tr>
    <tr>
      <th>njev</th>
      <td>0</td>
    </tr>
    <tr>
      <th>nlu</th>
      <td>0</td>
    </tr>
    <tr>
      <th>status</th>
      <td>1</td>
    </tr>
    <tr>
      <th>message</th>
      <td>A termination event occurred.</td>
    </tr>
    <tr>
      <th>success</th>
      <td>True</td>
    </tr>
  </tbody>
</table>
</div>



Check the final state.  If things have gone according to plan, the final value of `y` should be close to 0.


```python
# Solution

results.tail()
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>theta</th>
      <th>omega</th>
      <th>y</th>
      <th>v</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0.706147</th>
      <td>44.0182</td>
      <td>121.32</td>
      <td>0.327291</td>
      <td>-1.76573</td>
    </tr>
    <tr>
      <th>0.756147</th>
      <td>50.268</td>
      <td>128.609</td>
      <td>0.236982</td>
      <td>-1.84499</td>
    </tr>
    <tr>
      <th>0.806147</th>
      <td>56.8724</td>
      <td>135.496</td>
      <td>0.142962</td>
      <td>-1.91405</td>
    </tr>
    <tr>
      <th>0.856147</th>
      <td>63.8093</td>
      <td>141.885</td>
      <td>0.0457628</td>
      <td>-1.97198</td>
    </tr>
    <tr>
      <th>0.879218</th>
      <td>67.1147</td>
      <td>144.631</td>
      <td>9.02056e-17</td>
      <td>-1.99469</td>
    </tr>
  </tbody>
</table>
</div>



Plot the results.

`theta` should increase and accelerate.


```python
def plot_theta(results):
    plot(results.theta, color='C0', label='theta')
    decorate(xlabel='Time (s)',
             ylabel='Angle (rad)')
plot_theta(results)
```

`y` should decrease and accelerate down.


```python
def plot_y(results):
    plot(results.y, color='C1', label='y')

    decorate(xlabel='Time (s)',
             ylabel='Length (m)')
    
plot_y(results)
```

Plot velocity as a function of time; is the yo-yo accelerating?


```python
# Solution

v = results.v # m / s
plot(v)
decorate(xlabel='Time (s)',
         ylabel='Velocity (m/s)')
```

Use `gradient` to estimate the derivative of `v`.  How does the acceleration of the yo-yo compare to `g`?


```python
# Solution

a = gradient(v)
plot(a)
decorate(xlabel='Time (s)',
         ylabel='Acceleration (m/$s^2$)')
```


```python

```

2c0f1a3033d6fb6b19f8b8d47b3dbd58fad8aa8a
ipynb
Jupyter Notebook
soln/yoyo_soln.ipynb
kanhaiyap/ModSimPy
af16c079ec398ff9b3822d3dcda75873ce900ced
MIT
2019-04-27T22:43:12.000Z
2019-11-11T15:12:23.000Z
soln/yoyo_soln.ipynb
ffriass/ModSimPy
c36a476a20042acb33773e47d12aea5b0c413e60
MIT
2019-10-09T18:50:22.000Z
2022-03-21T01:39:48.000Z
soln/yoyo_soln.ipynb
ffriass/ModSimPy
c36a476a20042acb33773e47d12aea5b0c413e60
MIT
Qwen/Qwen-72B
1. YES
2. YES