Modeling and Simulation in Python

Chapter 10 Example: Rigid Pendulum

Copyright 2017 Allen Downey

License: Creative Commons Attribution 4.0 International

In [1]:
# If you want the figures to appear in the notebook, 
# and you want to interact with them, use
# %matplotlib notebook

# If you want the figures to appear in the notebook, 
# and you don't want to interact with them, use
# %matplotlib inline

# If you want the figures to appear in separate windows, use
# %matplotlib qt5

# tempo switch from one to another, you have to select Kernel->Restart

%matplotlib notebook

from modsim import *


This notebook solves the Spider-Man problem from spiderman.ipynb, demonstrating a different development process for physical simulations.

In pendulum_sympy, we derive the equations of motion for a rigid pendulum without drag, yielding:

$ \ddot{x} = \frac{x}{x^{2} + y^{2}} \left(g y - vx^{2} - vy^{2}\right) $

$ \ddot{y} = - \frac{1}{x^{2} + y^{2}} \left(g x^{2} + y \left(vx^{2} + vy^{2}\right)\right) $

We'll use the same conditions we saw in spiderman.ipynb

In [2]:
condition = Condition(g = 9.8,
                      m = 75,
                      area = 1,
                      rho = 1.2,
                      v_term = 60,
                      duration = 30,
                      length0 = 100,
                      angle = (270 - 45),
                      k = 20)

Now here's a version of make_system that takes a Condition object as a parameter.

make_system uses the given value of v_term to compute the drag coefficient C_d.

In [3]:
def make_system(condition):
    """Makes a System object for the given conditions.
    condition: Condition with height, g, m, diameter, 
               rho, v_term, and duration
    returns: System with init, g, m, rho, C_d, area, and ts
    theta = np.deg2rad(angle)
    x, y = pol2cart(theta, length0)
    P = Vector(x, y)
    V = Vector(0, 0)
    init = State(x=P.x, y=P.y, vx=V.x, vy=V.y)
    C_d = 2 * m * g / (rho * area * v_term**2)
    ts = linspace(0, duration, 501)
    return System(init=init, g=g, m=m, rho=rho,
                  C_d=C_d, area=area, length0=length0,
                  k=k, ts=ts)

Let's make a System

In [4]:
system = make_system(condition)

init x -70.71067811865477 dimensionless y -...
g 9.8
m 75
rho 1.2
C_d 0.340278
area 1
length0 100
k 20
ts [0.0, 0.06, 0.12, 0.18, 0.24, 0.3, 0.36, 0.42,...

In [5]:

x -70.71067811865477 dimensionless
y -70.71067811865474 dimensionless
vx 0 dimensionless
vy 0 dimensionless

To write the slope function, we can get the expressions for ax and ay directly from SymPy and plug them in.

In [6]:
def slope_func(state, t, system):
    """Computes derivatives of the state variables.
    state: State (x, y, x velocity, y velocity)
    t: time
    system: System object with length0, m, k
    returns: sequence (vx, vy, ax, ay)
    x, y, vx, vy = state

    ax = x*(g*y - vx**2 - vy**2)/(x**2 + y**2)
    ay = -(g*x**2 + y*(vx**2 + vy**2))/(x**2 + y**2)

    return vx, vy, ax, ay

As always, let's test the slope function with the initial conditions.

In [7]:
slope_func(system.init, 0, system)

(<Quantity(0, 'dimensionless')>,
 <Quantity(0, 'dimensionless')>,
 <Quantity(4.9, 'dimensionless')>,
 <Quantity(-4.900000000000003, 'dimensionless')>)

And then run the simulation.

In [8]:
%time run_odeint(system, slope_func)

CPU times: user 132 ms, sys: 0 ns, total: 132 ms
Wall time: 130 ms

Visualizing the results

We can extract the x and y components as Series objects.

In [9]:
xs = system.results.x
ys = system.results.y

The simplest way to visualize the results is to plot x and y as functions of time.

In [10]:
plot(xs, label='x')
plot(ys, label='y')

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

We can plot the velocities the same way.

In [11]:
vxs = system.results.vx
vys = system.results.vy

In [12]:
plot(vxs, label='vx')
plot(vys, label='vy')

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

Another way to visualize the results is to plot y versus x. The result is the trajectory through the plane of motion.

In [13]:
plot(xs, ys, label='trajectory')

decorate(xlabel='x position (m)',
         ylabel='y position (m)')

We can also animate the trajectory. If there's an error in the simulation, we can sometimes spot it by looking at animations.

In [17]:
decorate(xlabel='x position (m)',
         ylabel='y position (m)',
         xlim=[-100, 100],
         ylim=[-200, -50],

for x, y in zip(xs, ys):
    plot(x, y, 'bo', update=True)

Here's a function that encapsulates that code and runs the animation in (approximately) real time.

In [15]:
def animate2d(xs, ys, speedup=1):
    """Animate the results of a projectile simulation.
    xs: x position as a function of time
    ys: y position as a function of time
    speedup: how much to divide `dt` by
    # get the time intervals between elements
    ts = xs.index
    dts = np.diff(ts)
    dts = np.append(dts, 0)

    # decorate the plot
    decorate(xlabel='x position (m)',
             ylabel='y position (m)',
             xlim=[xs.min(), xs.max()],
             ylim=[ys.min(), ys.max()],

    # loop through the values
    for x, y, dt in zip(xs, ys, dts):
        plot(x, y, 'bo', update=True)
        sleep(dt / speedup)

In [16]:
animate2d(system.results.x, system.results.y)

In [ ]: