This chapter is available as a Jupyter notebook where you can read the text, run the code, and work on the exercises. Click here to access the notebooks:

In this chapter and the next we’ll model systems that involve rotating objects. In general, rotation is complicated. In three dimensions, objects can rotate around three axes and many objects are easier to spin around some axes than others. If the configuration of an object changes over time, it might become easier or harder to spin, which explains the surprising dynamics of gymnasts, divers, ice skaters, etc. And when you apply a twisting force to a rotating object, the effect is often contrary to intuition. For an example, see this video on gyroscopic precession:

We will not take on the physics of rotation in all its glory; rather, we will focus on simple scenarios where all rotation and all twisting forces are around a single axis. In that case, we can treat some vector quantities as if they were scalars; that is, simple numbers.

The fundamental ideas in these examples are angular velocity, angular acceleration, torque, and moment of inertia. If you are not already familiar with these concepts, don’t worry; I will define them as we go along, and I will point to additional reading.

The Physics of Toilet Paper#

As an example of a system with rotation, we’ll simulate the manufacture of a roll of toilet paper, as shown in this video: Starting with a cardboard tube at the center, we will roll up 47 m of paper, a typical length for a roll of toilet paper in the U.S. (see

The following figure shows a diagram of the system: \(r\) represents the radius of the roll at a point in time. Initially, \(r\) is the radius of the cardboard core, \(R_{min}\). When the roll is complete, \(r\) is \(R_{max}\).

Diagram of a roll of toilet paper, showing change in paper length as a result of a small rotation, .

I’ll use \(\theta\) to represent the total rotation of the roll in radians. In the diagram, \(d\theta\) represents a small increase in \(\theta\), which corresponds to a distance along the circumference of \(r~d\theta\).

I’ll use \(y\) to represent the total length of paper that’s been rolled. Initially, \(\theta=0\) and \(y=0\). For each small increase in \(\theta\), there is a corresponding increase in \(y\):

\[dy = r~d\theta\]

If we divide both sides by a small increase in time, \(dt\), we get a differential equation for \(y\) as a function of time.

\[\frac{dy}{dt} = r \frac{d\theta}{dt}\]

As we roll up the paper, \(r\) increases. Assuming it increases by a fixed amount per revolution, we can write

\[dr = k~d\theta\]

where \(k\) is an unknown constant we’ll have to figure out. Again, we can divide both sides by \(dt\) to get a differential equation in time:

\[\frac{dr}{dt} = k \frac{d\theta}{dt}\]

Finally, let’s assume that \(\theta\) increases at a constant rate of \(\omega = 300\) rad/s (about 2900 revolutions per minute):

\[\frac{d\theta}{dt} = \omega\]

This rate of change is called an angular velocity. Now we have a system of differential equations we can use to simulate the system.

Setting Parameters#

Here are the parameters of the system:

Rmin = 0.02      # m
Rmax = 0.055     # m
L = 47           # m
omega = 300      # rad / s

Rmin and Rmax are the initial and final values for the radius, r. L is the total length of the paper. omega is the angular velocity in radians per second.

Figuring out k is not easy, but we can estimate it by pretending that r is constant and equal to the average of Rmin and Rmax:

Ravg = (Rmax + Rmin) / 2

In that case, the circumference of the roll is also constant:

Cavg = 2 * np.pi * Ravg

And we can compute the number of revolutions to roll up length L, like this.

revs = L / Cavg

Converting rotations to radians, we can estimate the final value of theta.

theta = 2 * np.pi * revs

Finally, k is the total change in r divided by the total change in theta.

k_est = (Rmax - Rmin) / theta

At the end of the chapter, we’ll derive k analytically, but this estimate is enough to get started.

Simulating the System#

The state variables we’ll use are theta, y, and r. Here are the initial conditions:

init = State(theta=0, y=0, r=Rmin)

And here’s a System object with init and t_end:

system = System(init=init, t_end=10)

Now we can use the differential equations from the previous section to write a slope function:

def slope_func(t, state, system):
    theta, y, r = state
    dydt = r * omega
    drdt = k_est * omega
    return omega, dydt, drdt

As usual, the slope function takes a time stamp, a State object, and a System object.

The job of the slope function is to compute the time derivatives of the state variables. The derivative of theta is angular velocity, omega. The derivatives of y and r are given by the differential equations we derived.

And as usual, we’ll test the slope function with the initial conditions.

slope_func(0, system.init, system)
(300, 6.0, 0.008377659574468085)

We’d like to stop the simulation when the length of paper on the roll is L. We can do that with an event function that passes through 0 when y equals L:

def event_func(t, state, system):
    theta, y, r = state
    return L - y

We can test it with the initial conditions:

event_func(0, system.init, system)

Now let’s run the simulation:

results, details = run_solve_ivp(system, slope_func,
'A termination event occurred.'

Here are the last few time steps.

theta y r
4.010667 1203.200000 44.277760 0.05360
4.052444 1215.733333 44.951740 0.05395
4.094222 1228.266667 45.630107 0.05430
4.136000 1240.800000 46.312860 0.05465
4.177778 1253.333333 47.000000 0.05500

The time it takes to complete one roll is about 4.2 seconds, which is consistent with what we see in the video.


The final value of y is 47 meters, as expected.

final_state = results.iloc[-1] 

The final value of r is 0.55 m, which is Rmax.


The total number of rotations is close to 200, which seems plausible.

radians = final_state.theta
rotations = radians / 2 / np.pi

As an exercise, we’ll see how fast the paper is moving. But first, let’s take a closer look at the results.

Plotting the Results#

Here’s what theta looks like over time.

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

theta grows linearly, as we should expect with constant angular velocity.

Here’s what r looks like over time.

def plot_r(results):
    results.r.plot(color='C2', label='r')

    decorate(xlabel='Time (s)',
             ylabel='Radius (m)')

r also increases linearly.

Here’s what y looks like over time.

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

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

Since the derivative of y depends on r, and r is increasing, y grows with increasing slope.

In the next section, we’ll see that we could have solved these differential equations analytically. However, it is often useful to start with simulation as a way of exploring and checking assumptions.

Analytic Solution#

Since angular velocity is constant:

\[\frac{d\theta}{dt} = \omega \quad\quad (1)\]

We can find \(\theta\) as a function of time by integrating both sides:

\[\theta(t) = \omega t\]

Similarly, we can solve this equation

\[\frac{dr}{dt} = k \omega\]

to find

\[r(t) = k \omega t + R_{min}\]

Then we can plug the solution for \(r\) into the equation for \(y\):

\[\begin{split}\begin{aligned} \frac{dy}{dt} & = r \omega \quad\quad (2) \\ & = \left[ k \omega t + R_{min} \right] \omega \nonumber\end{aligned}\end{split}\]

Integrating both sides yields:

\[y(t) = \left[ k \omega t^2 / 2 + R_{min} t \right] \omega\]

So \(y\) is a parabola, as you might have guessed.

We can also use these equations to find the relationship between \(y\) and \(r\), independent of time, which we can use to compute \(k\). Dividing Equations 1 and 2 yields

\[\frac{dr}{dy} = \frac{k}{r}\]

Separating variables yields

\[r~dr = k~dy\]

Integrating both sides yields

\[r^2 / 2 = k y + C\]

Solving for \(y\), we have

\[y = \frac{1}{2k} (r^2 - C) \label{eqn3}\]

When \(y=0\), \(r=R_{min}\), so

\[R_{min}^2 / 2 = C\]

When \(y=L\), \(r=R_{max}\), so

\[L = \frac{1}{2k} (R_{max}^2 - R_{min}^2)\]

Solving for \(k\) yields

\[k = \frac{1}{2L} (R_{max}^2 - R_{min}^2) \label{eqn4}\]

Plugging in the values of the parameters yields 2.8e-5 m/rad, the same as the “estimate” we computed in Section xxx.

k = (Rmax**2 - Rmin**2) / (2 * L)

In this case the estimate turns out to be exact.


This chapter introduces rotation, starting with an example where angular velocity is constant. We simulated the manufacture of a roll of toilet paper, then we solved the same problem analytically.

In the next chapter, we’ll see a more interesting example where angular velocity is not constant. And we’ll introduce three new concepts: torque, angular acceleration, and moment of inertia.

But first, you might want to work on the following exercise.


This chapter is available as a Jupyter notebook where you can read the text, run the code, and work on the exercises. You can access the notebooks at

Exercise 1#

Since we keep omega constant, the linear velocity of the paper increases with radius. We can use gradient to estimate the derivative of results.y.

dydt = gradient(results.y)

Here’s what the result looks like.

decorate(xlabel='Time (s)',
         ylabel='Linear velocity (m/s)')

With constant angular velocity, linear velocity is increasing, reaching its maximum at the end.

max_linear_velocity = dydt.iloc[-1]

Now suppose this peak velocity is the limiting factor; that is, we can’t move the paper any faster than that. In that case, we might be able to speed up the process by keeping the linear velocity at the maximum all the time.

Write a slope function that keeps the linear velocity, dydt, constant, and computes the angular velocity, omega, accordingly. Then, run the simulation and see how much faster we could finish rolling the paper.

Hide code cell content
# Solution

def slope_func(t, state, system):
    theta, y, r = state
    dydt = max_linear_velocity
    omega = dydt / r
    drdt = k * omega
    return omega, dydt, drdt
Hide code cell content
# Solution

slope_func(0, system.init, system)
(822.375000000005, 16.4475000000001, 0.022965259308510778)
Hide code cell content
# Solution

results, details = run_solve_ivp(system, slope_func, 
'A termination event occurred.'
Hide code cell content
# Solution

t_final = results.index[-1]
Hide code cell content
# Solution

Hide code cell content
# Solution

Hide code cell content
# Solution