Tutorial 7 -- Animation of a Dynamic System

In this tutorial, we want to produce an animation of a simple pendulum governed by the following ODE: $$ \ddot{\theta} + \gamma \dot{\theta} + \omega_0^2 \sin\theta = 0 \qquad\qquad \text{(1)} $$ where $\omega_0 = \sqrt{g/l}$ is the natural frequency of small amplitude oscillations, and $\gamma$ is the damping coefficient. Given initial values of the angle and angular velocity $(\theta_0, \dot{\theta}_0)$ at time $t=0$, the solution of the ODE gives the functions $(\theta(t), \dot{\theta}(t))$ versus time $t$. We want to find a numerical solution of (1) and use this solution to produce an animation in physical space.

To accomplish this, we will use elements of tutorials 5 and 6.

Step 1: we import the necessary modules

import numpy as np
from matplotlib import pyplot as plt
from matplotlib import animation
import scipy.integrate as integrate

Step 2: we define the parameters.

g = 9.81 # [m/s^2]
l = 0.5 # [m]
gamma = 0. # [1/s] damping coefficient
omega0 = np.sqrt(g/l)

Step 3: we define the initial conditions, that is, initial values for angle theta and angular velocity omega.

theta= 170 # [degree]  # initial angle
theta = theta*np.pi/180
omega = 0  # initial angular velocity

Step 4: we want to get an idea of a time scale as a function of the initial angle. To do this, we compute the period of oscillations Tper by a quadrature. We then use this value to set the maximum time Tmax of the simulation as a multiple of Tper.

def f(x):
    return 1/np.sqrt(np.cos(x)-np.cos(theta))

Tper = integrate.quad(f ,0,theta )
T0 = 2.0*np.sqrt(2.0)*Tper[0]/omega0
Tmax = 8*T0
print('Tmax=',Tmax)

Step 5: we define graphical and simulation parameters.

xmax =  1.05*l
ymax =  1.05*l
Nstep= 3600   # number of time steps / number of frames
dt= Tmax/Nstep  # time increment
d=0.05  # size of fixed support

Step 6: we set up the figure, the axis, and the line object we want to animate. This line object is empty for now. We also draw a rectangular object of size 2d x 2d centered at the origin.

fig = plt.figure()
ax = plt.axes(xlim=(-xmax, xmax), ylim=(-ymax, ymax))

line, = ax.plot([], [], 'o', ls='-', ms=8, markevery=[0,-1])
xR,yR = -d,-d # lower left corner of rectangle
dx,dy= 2*d, 2*d  # dimensions of rectangle
rect = plt.Rectangle((xR,yR), dx,dy, facecolor= 'black', edgecolor= 'black')
ax.add_patch(rect)
plt.axis('off')
ax.set_aspect('equal')

Step 7: we define a function which models the equations of motion

def F(X,t):
    [theta, omega] = X
    dthetadt = omega
    domegadt = -g/l* np.sin(theta) - gamma *omega
    return dthetadt, domegadt

Step 8: we define a function to set up the initialization of the plot of each frame.

def init():
    line.set_data([], [])
    return line,

Step 9: we define the animation function. At each value of integer i, we compute the pendulum position by integration of the ODE from time timeIn = i*dt to time timeOut= (i+1)*dt. To do this, we call scipy ODE integrator odeint. Once we have new values of variables theta and omega, we can compute the new Cartesian coordinates of the pendulum and update the line object.

def animate(i):
    
    global theta, omega
    timeIn= i*dt
    timeOut= (i+1)*dt
    Xinit = [theta, omega]
    odeout = integrate.odeint(F, Xinit, [timeIn, timeOut])[1]
    theta = odeout[0] 
    omega = odeout[1]
    theta = (theta + np.pi) % (2*np.pi) - np.pi # bring theta back to [-pi, pi]
    x,y = l*np.sin(theta), -l*np.cos(theta)
    line.set_data([0, x], [0,y]) # modify the pendulum's position    
    return line,

Step 10: we can now call the animator with the blit variable set to True. This computes the anim object. We then save the video in a mp4 file. We can view this video in our browser.

anim = animation.FuncAnimation(fig, animate, init_func=init,frames=Nstep, interval=20, blit=True)
anim.save('pendulum_animation_odeint.mp4', fps=30, extra_args=['-vcodec', 'libx264'])
plt.show()

The complete script can be found here. The animation with l =0.5 [m], theta = 270 [deg], omega= 0, gamma = 0, Tmax= 2*Tper and 800 frames gives this video:

If we now add damping with gamma =0.05 and Tmax= 8*Tper, 3600 frames gives this video

It is important to properly integrate the ODE. If we attempt to do the same simulation with a simplistic Euler’s method, we get an inaccurate motion, as shown by this video This was obtained with the same initial conditions as chosen in the first video of this tutorial. The script is here.

[Previous Tutorial]