Using Python and a couple of popular modules, SciPy and Matplotlib, you can simulate physical systems like the double pendulum shown above.

Each system will be represented as a class that consists of two parts:

• The differential equations of the system
• An initial state

We can plug these two items into a numerical integrator supplied by SciPy, odeint, and render the output using matplotlib.

The Single Pendulum

We'll start with the basic case, the single uncoupled pendulum. With Z representing the state of the system, and theta being the angular position of the pendulum arm,

$\tag{1} \vec{Z} = \begin{bmatrix} \theta \\ \dot{\theta} \end{bmatrix}$

And the first derivative,

$\tag{2} \dot{\vec{Z}} = \begin{bmatrix} \dot{\theta} \\ -\frac{g}{l}\sin{\theta} \end{bmatrix}$

For the initial state, all that is required is a starting angle.

The Python class we are creating has three methods:

1. init : Defines the initial state as well as the physical parameters (strength of gravity, length of arm)
2. equation : Returns the system of equations outlined above. We have this because the odeint method requires a callable function that returns the system of equations.
3. solve_ODE : Uses the odeint function to numerically solve the differential equations and then converts the data to x,y coordinates.

class Pendulum():

 """init_state has the form [theta, thetaDot]"""
   def __init__(self,
              g = 9.8,
              L = 2.0):

     self.init_state = [np.radians(35.0), 0]
     self.time = np.arange(0, 50.0, 0.025) 
     self.g = g
     self.L = L

    def equation(self, z0,t):

       theta, thetaDot = z0
       dzdt = [thetaDot, -(self.g/self.L)*sin(theta)]
       return dzdt

    def solve_ODE(self):
      
       self.state = odeint(self.equation, self.init_state, self.time)
       x1 = sin(self.state[:, 0])*self.L
       y1 = -1*cos(self.state[:, 0])*self.L

       return x1, y1

Here is are a few animations using of the output. Each animation consists of a red and blue pendulum both starting from the same initial angle.

• The red pendulum represents the model we get from using the odeint and the model outlined above.
• The blue pendulum represents the data that we would get if we integrated the equations using the less accurate "small-angle approximation" that is often used to solve the equations without a computer integrator.

The Coupled Pendulum

The coupled pendulum consists of two pendulums linked by an elastic spring. Similarly, we start by defining the state equations,

Since there are now two pendulums, we must define a second variable to represent the displacement of the second arm, phi.

$\tag{3} \vec{Z} = \begin{bmatrix} \theta \\ \dot{\theta} \\ \phi \\ \dot{\phi} \end{bmatrix}$

The first derivative,

$\tag{4} \dot{\vec{Z}} = \begin{bmatrix} \dot{\theta} \\ \ddot{\theta} \\ \dot{\phi} \\ \ddot{\phi} \end{bmatrix}$

I will skip some math and provide the expressions:

$\tag{5} \ddot{\theta} = \frac{kl_2\sin{\phi}\sin{\theta}(m_1(l_1\dot{\theta}^{2}-g) - kl_1)}{m_1l_1\cos{\theta}}$ $\tag{6} \ddot{\phi} = \frac{kl_1\sin{\theta}\sin{\phi}(m_2(l_2\dot{\phi}^{2}-g) - kl_2)}{m_2l_2\cos{\phi}}$

Let's create a new class for the coupled pendulum with this new system.

class coupledPendulum():

   def __init__(self,
    g = 9.8,
    L1 = 1.5,
    L2 = 1.5,
    M1 = 1.0,
    M2 = 1.0,
    k = 0.5):

 self.init_state = [np.radians(25.0), 0, np.radians(0.0), 0]
 self.params = (L1, L2, M1, M2, g, k)
 self.time = np.arrange(0, 50.0, 0.025)

 def equation(self, z0,t):
    (L1, L2, M1, M2, g, k) = self.params
    theta, thetaDot, phi, phiDot = z0

    dzdt = [
     thetaDot,
    (sin(theta)*(M1*(L1*thetaDot*thetaDot-g)-k*L1)+k*L2*sin(phi))/(M1*L1*cos(theta)),
    phiDot,
    (sin(phi)*(M2*(L2*phiDot*phiDot-g)-k*L2)+k*L1*sin(theta))/(M2*L2*cos(phi))
    ]

    return dzdt 

 def solve_ODE(self):
    self.state = odeint(self.equation, self.init_state, self.time)

    """convert data into (x,y) coordinates"""
    x1 = sin(self.state[:, 0])*self.params[0]
    y1 = -1*cos(self.state[:, 0])*self.params[0]
    x2 = sin(self.state[:, 2])*self.params[1]
    y2 = -1*cos(self.state[:, 2])*self.params[1]
    return x1, y1, x2, y2    

The Double Pendulum

Lastly, there is the double pendulum. As with the coupled pendulum, we have two coordinates representing the angular displacements of the two different arms. In this case, the two arms are joined end-to-end.

$\tag{7} \vec{Z} = \begin{bmatrix} \theta \\ \dot{\theta} \\ \phi \\ \dot{\phi} \end{bmatrix}$

and again,

$\tag{8} \dot{\vec{Z}} = \begin{bmatrix} \dot{\theta} \\ \ddot{\theta} \\ \dot{\phi} \\ \ddot{\phi} \end{bmatrix}$

The expressions for ${\ddot{\theta}}$ and ${\ddot{\phi}}$ are pretty verbose and are given in the code below. If you are interested in how to get to them yourself, read about the Euler-Lagrange method here.

The double pendulum class:

class doublePendulum():

    def __init__(self,
                 g = 9.8,
                 L1 = 2.0,
                 L2 = 1.0,
                 M1 = 1.0,
                 M2 = 5.0):

        self.init_state = [np.radians(120.0), np.radians(0), np.radians(89.0), np.radians(0)]
        self.params = (L1, L2, M1, M2, g)
        self.time = np.arrange(0, 50.0, 0.025)

    def equation(self, y0, t):

        (L1, L2, M1, M2, g) = self.params
        theta, thetaDot, phi, phiDot = y0

        delsin = sin(theta - phi)
        delcos = cos(theta - phi)

        dydx = [thetaDot,
                (-M2 * L1 * thetaDot**2 * delsin * delcos + g * M2 * sin(phi) * delcos - M2 * L2 * phiDot**2 * delsin - (M1 + M2) * g * sin(theta))/(L1 * (M1 + M2) - M2 * L1 * delcos**2),
                phiDot,
                (M2 * L2 * phiDot**2 * delsin * delcos + g * sin(theta) * delcos * (M1 + M2) + L1 * thetaDot**2 * delsin * (M1 + M2) - g * sin(phi) * (M1 + M2))/(L2 * (M1 + M2) - M2 * L2 * delcos**2)
                ]

        return dydx

    def solve_ODE(self):

        self.state = odeint(self.equation, self.init_state, self.time)
        x1 = sin(self.state[:, 0])*self.params[0]
        y1 = -1*cos(self.state[:, 0])*self.params[0]
        x2 = x1 + sin(self.state[:, 2])*self.params[1]
        y2 = y1 + -1*cos(self.state[:, 2])*self.params[1]

        return x1, y1, x2, y2

Below are two animations of a double pendulum with two different initial states, an 89 degree release and a 90 degree release. You can see how quickly the two deviate from each other even with such a minor difference.