# Modeling and Simulation in Python

Chapter 25

``````

In [1]:

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

``````

### Teapots and Turntables

Tables in Chinese restaurants often have a rotating tray or turntable that makes it easy for customers to share dishes. These turntables are supported by low-friction bearings that allow them to turn easily and glide. However, they can be heavy, especially when they are loaded with food, so they have a high moment of inertia.

Suppose I am sitting at a table with a pot of tea on the turntable directly in front of me, and the person sitting directly opposite asks me to pass the tea. I push on the edge of the turntable with 1 Newton of force until it has turned 0.5 radians, then let go. The turntable glides until it comes to a stop 1.5 radians from the starting position. How much force should I apply for a second push so the teapot glides to a stop directly opposite me?

The following figure shows the scenario, where `F` is the force I apply to the turntable at the perimeter, perpendicular to the moment arm, `r`, and `tau` is the resulting torque. The blue circle near the bottom is the teapot.

We'll answer this question in these steps:

1. We'll use the results from the first push to estimate the coefficient of friction for the turntable.

2. We'll use that coefficient of friction to estimate the force needed to rotate the turntable through the remaining angle.

Our simulation will use the following parameters:

1. The radius of the turntable is 0.5 meters, and its weight is 7 kg.

2. The teapot weights 0.3 kg, and it sits 0.4 meters from the center of the turntable.

As usual, I'll get units from Pint.

``````

In [2]:

m = UNITS.meter
s = UNITS.second
kg = UNITS.kilogram
N = UNITS.newton

``````
``````

Out[2]:

newton

``````

And store the parameters in a `Params` object.

``````

In [3]:

mass_disk=7*kg,
mass_pot=0.3*kg,
force=1*N,
torque_friction=0.2*N*m,
t_end=20*s)

``````
``````

Out[3]:

.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}

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

text-align: right;
}

values

0.5 meter

mass_disk
7 kilogram

0.4 meter

mass_pot
0.3 kilogram

force
1 newton

torque_friction
0.2 meter * newton

theta_end

t_end
20 second

``````

`make_system` creates the initial state, `init`, and computes the total moment of inertia for the turntable and the teapot.

``````

In [4]:

def make_system(params):
"""Make a system object.

params: Params object

returns: System object
"""
mass_disk, mass_pot = params.mass_disk, params.mass_pot

I_disk = mass_disk * radius_disk**2 / 2

return System(params, init=init, I=I_disk+I_pot)

``````

Here's the `System` object we'll use for the first phase of the simulation, while I am pushing the turntable.

``````

In [5]:

system1 = make_system(params)

``````
``````

Out[5]:

.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}

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

text-align: right;
}

values

0.5 meter

mass_disk
7 kilogram

0.4 meter

mass_pot
0.3 kilogram

force
1 newton

torque_friction
0.2 meter * newton

theta_end

t_end
20 second

init

I
0.923 kilogram * meter ** 2

``````

### Simulation

When I stop pushing on the turntable, the angular acceleration changes abruptly. We could implement the slope function with an `if` statement that checks the value of `theta` and sets `force` accordingly. And for a coarse model like this one, that might be fine. But we will get more accurate results if we simulate the system in two phases:

1. During the first phase, force is constant, and we run until `theta` is 0.5 radians.

2. During the second phase, force is 0, and we run until `omega` is 0.

Then we can combine the results of the two phases into a single `TimeFrame`.

Here's the slope function we'll use:

``````

In [6]:

def slope_func(state, t, system):
"""Computes the derivatives of the state variables.

state: State object
t: time
system: System object

returns: sequence of derivatives
"""
theta, omega = state
torque_friction, I = system.torque_friction, system.I

torque = radius_disk * force - torque_friction
alpha = torque / I

return omega, alpha

``````

As always, we'll test the slope function before running the simulation.

``````

In [7]:

slope_func(system1.init, 0, system1)

``````
``````

Out[7]:

0.32502708559046584 <Unit('newton / kilogram / meter')>)

``````

Here's an event function that stops the simulation when `theta` reaches `theta_end`.

``````

In [8]:

def event_func1(state, t, system):
"""Stops when theta reaches theta_end.

state: State object
t: time
system: System object

returns: difference from target
"""
theta, omega = state
return theta - system.theta_end

``````
``````

In [9]:

event_func1(system1.init, 0, system1)

``````
``````

Out[9]:

``````

Now we can run the first phase.

``````

In [10]:

results1, details1 = run_ode_solver(system1, slope_func, events=event_func1)
details1

``````
``````

Out[10]:

.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}

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

text-align: right;
}

values

success
True

message
A termination event occurred.

``````

And look at the results.

``````

In [11]:

results1.tail()

``````
``````

Out[11]:

.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}

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

text-align: right;
}

theta
omega

1.000000

1.200000

1.400000

1.600000

1.751961

``````

### Phase 2

Before we run the second phase, we have to extract the final time and state of the first phase.

``````

In [12]:

t_0 = results1.last_label() * s

``````
``````

Out[12]:

1.7519607843137255 second

``````

And make an initial `State` object for Phase 2.

``````

In [13]:

init2 = results1.last_row()

``````
``````

Out[13]:

.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}

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

text-align: right;
}

values

theta

omega

``````

And a new `System` object with zero force.

``````

In [14]:

system2 = System(system1, t_0=t_0, init=init2, force=0*N)

``````
``````

Out[14]:

.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}

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

text-align: right;
}

values

0.5 meter

mass_disk
7 kilogram

0.4 meter

mass_pot
0.3 kilogram

force
0 newton

torque_friction
0.2 meter * newton

theta_end

t_end
20 second

init
ome...

I
0.923 kilogram * meter ** 2

t_0
1.7519607843137255 second

``````

Here's an event function that stops when angular velocity is 0.

``````

In [15]:

def event_func2(state, t, system):
"""Stops when omega is 0.

state: State object
t: time
system: System object

returns: omega
"""
theta, omega = state
return omega

``````
``````

In [16]:

event_func2(system2.init, 0, system2)

``````
``````

Out[16]:

``````

Now we can run the second phase.

``````

In [17]:

slope_func(system2.init, system2.t_0, system2)

``````
``````

Out[17]:

-0.21668472372697725 <Unit('newton / kilogram / meter')>)

``````
``````

In [18]:

results2, details2 = run_ode_solver(system2, slope_func, events=event_func2)
details2

``````
``````

Out[18]:

.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}

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

text-align: right;
}

values

success
True

message
A termination event occurred.

``````

And check the results.

``````

In [19]:

results2.tail()

``````
``````

Out[19]:

.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}

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

text-align: right;
}

theta
omega

3.751961

3.951961

4.151961

4.351961

4.379902

``````

Pandas provides `combine_first`, which combines `results1` and `results2`.

``````

In [20]:

results = results1.combine_first(results2)
results.tail()

``````
``````

Out[20]:

.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}

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

text-align: right;
}

theta
omega

3.751961

3.951961

4.151961

4.351961

4.379902

``````

Now we can plot `theta` for both phases.

``````

In [21]:

def plot_theta(results):
plot(results.theta, label='theta')
decorate(xlabel='Time (s)',

plot_theta(results)

``````
``````

``````

And `omega`.

``````

In [22]:

def plot_omega(results):
plot(results.omega, label='omega', color='C1')
decorate(xlabel='Time (s)',

plot_omega(results)

``````
``````

``````
``````

In [23]:

subplot(2, 1, 1)
plot_theta(results)
subplot(2, 1, 2)
plot_omega(results)
savefig('figs/chap25-fig01.pdf')

``````
``````

Saving figure to file figs/chap25-fig01.pdf

``````

### Estimating friction

Let's take the code from the previous section and wrap it in a function.

``````

In [24]:

def run_two_phases(force, torque_friction, params):
"""Run both phases.

force: force applied to the turntable
torque_friction: friction due to torque
params: Params object

returns: TimeFrame of simulation results
"""
# put the specified parameters into the Params object
params = Params(params, force=force, torque_friction=torque_friction)

# run phase 1
system1 = make_system(params)
results1, details1 = run_ode_solver(system1, slope_func,
events=event_func1)

# get the final state from phase 1
t_0 = results1.last_label() * s
init2 = results1.last_row()

# run phase 2
system2 = System(system1, t_0=t_0, init=init2, force=0*N)
results2, details2 = run_ode_solver(system2, slope_func,
events=event_func2)

# combine and return the results
results = results1.combine_first(results2)
return TimeFrame(results)

``````

Let's test it with the same parameters.

``````

In [25]:

force = 1*N
torque_friction = 0.2*N*m
results = run_two_phases(force, torque_friction, params)
results.tail()

``````
``````

Out[25]:

.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}

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

text-align: right;
}

theta
omega

3.751961

3.951961

4.151961

4.351961

4.379902

``````

And check the results.

``````

In [26]:

theta_final = results.last_row().theta

``````
``````

Out[26]:

``````

Here's the error function we'll use with `root_bisect`.

It takes a hypothetical value for `torque_friction` and returns the difference between `theta_final` and the observed duration of the first push, 1.5 radian.

``````

In [27]:

def error_func1(torque_friction, params):
"""Error function for root_scalar.

torque_friction: hypothetical value
params: Params object

returns: offset from target value
"""
force = 1 * N
results = run_two_phases(force, torque_friction, params)
theta_final = results.last_row().theta
print(torque_friction, theta_final)
return theta_final - 1.5 * radian

``````

Testing the error function.

``````

In [28]:

guess1 = 0.1*N*m
error_func1(guess1, params)

``````
``````

0.1 meter * newton 2.491081016010591 radian

Out[28]:

``````
``````

In [29]:

guess2 = 0.3*N*m
error_func1(guess2, params)

``````
``````

0.3 meter * newton 0.8319164068661051 radian

Out[29]:

``````

And running `root_scalar`.

``````

In [30]:

res = root_bisect(error_func1, [guess1, guess2], params)

``````
``````

0.1 meter * newton 2.491081016010591 radian
0.3 meter * newton 0.8319164068661051 radian
0.2 meter * newton 1.247699599245727 radian
0.15000000000000002 meter * newton 1.664679140408952 radian
0.17500000000000002 meter * newton 1.4245119372074873 radian
0.16250000000000003 meter * newton 1.5351133644066652 radian
0.16875 meter * newton 1.4775125382367449 radian
0.16562500000000002 meter * newton 1.5057134910749645 radian
0.16718750000000002 meter * newton 1.4913395920086545 radian
0.16640625000000003 meter * newton 1.4984572231242101 radian
0.16601562500000003 meter * newton 1.502067905168568 radian
0.16621093750000004 meter * newton 1.5002582165443015 radian
0.16630859375000004 meter * newton 1.4993566348484146 radian
0.16625976562500006 meter * newton 1.4998071542109142 radian
0.16623535156250005 meter * newton 1.500032617476343 radian
0.16624755859375007 meter * newton 1.499919868872052 radian
0.16624145507812504 meter * newton 1.4999762389308366 radian
0.16623840332031253 meter * newton 1.5000044271426904 radian
0.1662399291992188 meter * newton 1.4999903327715458 radian
0.16623916625976565 meter * newton 1.499997379890813 radian
0.1662387847900391 meter * newton 1.5000009035001758 radian
0.16623897552490235 meter * newton 1.4999991416913503 radian
0.16623888015747074 meter * newton 1.5000000225947272 radian
0.16623892784118655 meter * newton 1.4999995821427796 radian
0.16623890399932864 meter * newton 1.4999998023686882 radian
0.1662388920783997 meter * newton 1.4999999124816914 radian

Out[30]:

.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}

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

text-align: right;
}

values

converged
True

root
0.16623888611793522 meter * newton

``````

The result is the coefficient of friction that yields a total rotation of 1.5 radian.

``````

In [31]:

torque_friction = res.root

``````
``````

Out[31]:

0.16623888611793522 meter newton

``````

Here's a test run with the estimated value.

``````

In [32]:

force = 1 * N
results = run_two_phases(force, torque_friction, params)
theta_final = get_last_value(results.theta)

``````
``````

Out[32]:

``````

Looks good.

### Animation

Here's a draw function we can use to animate the results.

``````

In [33]:

from matplotlib.patches import Circle
from matplotlib.patches import Arrow

def draw_func(state, t):
theta, omega = state

# draw a circle for the table

# draw a circle for the teapot
circle2 = Circle(center, 0.05, color='C1')

# make the aspect ratio 1
plt.axis('equal')

``````
``````

In [34]:

state = results.first_row()
draw_func(state, 0)

``````
``````

``````
``````

In [35]:

animate(results, draw_func)

``````
``````

``````

### Exercises

Now finish off the example by estimating the force that delivers the teapot to the desired position.

Write an error function that takes `force` and `params` and returns the offset from the desired angle.

``````

In [36]:

# Solution

def error_func2(force, params):
"""Error function for root_scalar.

force: hypothetical value
params: Params object

returns: offset from target value
"""
# notice that this function uses the global value of torque_friction
results = run_two_phases(force, torque_friction, params)
theta_final = get_last_value(results.theta)
print(force, theta_final)
remaining_angle = np.pi - 1.5
return theta_final - remaining_angle * radian

``````

Test the error function with `force=1`

``````

In [37]:

# Solution

guess1 = 0.5 * N
error_func2(guess1, params)

``````
``````

Out[37]:

``````
``````

In [38]:

# Solution

guess2 = 2 * N
error_func2(guess2, params)

``````
``````

Out[38]:

``````

And run `root_bisect` to find the desired force.

``````

In [39]:

# Solution

res = root_bisect(error_func2, [guess1, guess2], params)

``````
``````

Out[39]:

.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}

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

text-align: right;
}

values

converged
True

root
1.0941684395074844 newton

``````
``````

In [40]:

force = res.root
results = run_two_phases(force, torque_friction, params)
theta_final = get_last_value(results.theta)

``````
``````

Out[40]:

``````
``````

In [41]:

remaining_angle = np.pi - 1.5

``````
``````

Out[41]:

1.6415926535897931

``````

Exercise: Now suppose my friend pours 0.1 kg of tea and puts the teapot back on the turntable at distance 0.3 meters from the center. If I ask for the tea back, how much force should they apply, over an arc of 0.5 radians, to make the teapot glide to a stop back in front of me? You can assume that torque due to friction is proportional to the total mass of the teapot and the turntable.

``````

In [42]:

# Solution

mass_before = params.mass_pot + params.mass_disk

``````
``````

Out[42]:

7.3 kilogram

``````
``````

In [43]:

# Solution

``````
``````

Out[43]:

.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}

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

text-align: right;
}

values

0.5 meter

mass_disk
7 kilogram

0.3 meter

mass_pot
0.2 kilogram

force
1 newton

torque_friction
0.2 meter * newton

theta_end

t_end
20 second

``````
``````

In [44]:

# Solution

mass_after = params2.mass_pot + params2.mass_disk

``````
``````

Out[44]:

7.2 kilogram

``````
``````

In [45]:

# Solution

torque_friction

``````
``````

Out[45]:

0.16623888611793522 meter newton

``````
``````

In [46]:

# Solution

torque_friction2 = torque_friction * mass_after / mass_before

``````
``````

Out[46]:

0.16396164110262107 meter newton

``````
``````

In [47]:

# Solution

guess = 2 * N
results = run_two_phases(guess, torque_friction2, params2)
results.tail()

``````
``````

Out[47]:

.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}

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

text-align: right;
}

theta
omega

5.630969

5.830969

6.030969

6.230969

6.287870

``````
``````

In [48]:

# Solution

subplot(2, 1, 1)
plot_theta(results)
subplot(2, 1, 2)
plot_omega(results)

``````
``````

``````
``````

In [49]:

# Solution

# Note: this is so similar to error_func2, it would be better
# to generalize it, but for expediency, I will make a modified
# verison.

def error_func3(force, params):
"""Error function for root_scalar.

force: hypothetical value
params: Params object

returns: offset from target value
"""
# notice that this function uses the global value of torque_friction2
results = run_two_phases(force, torque_friction2, params)
theta_final = get_last_value(results.theta)
print(force, theta_final)
return theta_final - remaining_angle

``````
``````

In [50]:

# Solution

guess1 = 1 * N
error_func3(guess1, params)

``````
``````

Out[50]:

``````
``````

In [51]:

# Solution

guess2 = 3 * N
error_func3(guess2, params)

``````
``````

Out[51]:

``````
``````

In [52]:

# Solution

res = root_bisect(error_func3, [guess1, guess2], params2)

``````
``````

Out[52]:

.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}

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

text-align: right;
}

values

converged
True

root
2.0648153424263 newton

``````
``````

In [53]:

# Solution

error_func3(res.root, params)

``````
``````

Out[53]: