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.