This tutorial aims at modelling and solving the yet classical but not so simple problem of the pendulum. A representiation is given bellow (source: Wikipedia).
The simple pendulumOn a mechanical point of view, the mass m is supposed to be concentrated at the lower end of the rigid arm. The length of the arm is noted l. The angle between the arm and the vertical direction is noted θ. A simple modelling using dynamics leads to:
Γ = P⃗ + T⃗
Where:
- Γ⃗ is the acceleration of the mass,
- P⃗ if the weight of the mass,
- T⃗ if the reaction force of the arm.
A projection of this equation along the direction perpendicular to the arm gives a more simple equation:
This equation is a second order, non linear ODE. The closed form solution is only known when the equation is linearized by assuming that θ is small enough to write that sin θ ≈ θ. In this tutorial, we will solve this problem with a numerical approach that does not require such simplification.
# Setup
%matplotlib nbagg
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.integrate import odeint
l = 1. # m
g = 9.81 # m/s**2
- This problem can be reformulated to match the standard formulation Ẋ = f(X, t):
- Write the function f in Python:
def f(X, t):
"""
The derivative function
"""
# To be completed
return
Solve the problem with Euler, RK4 and ODEint integrators and compare the results. First assume that the pendulum is released with no speed (θ̇ = 0o/s) at θ = 10o. The time discretization will be as follows:
- duration: 10 s,
- time step: 0.01 s.
def Euler(func, X0, t):
"""
Euler integrator.
"""
dt = t[1] - t[0]
nt = len(t)
X = np.zeros([nt, len(X0)])
X[0] = X0
for i in range(nt-1):
X[i+1] = X[i] + func(X[i], t[i]) * dt
return X
def RK4(func, X0, t):
"""
Runge and Kutta 4 integrator.
"""
dt = t[1] - t[0]
nt = len(t)
X = np.zeros([nt, len(X0)])
X[0] = X0
for i in range(nt-1):
k1 = func(X[i], t[i])
k2 = func(X[i] + dt/2. * k1, t[i] + dt/2.)
k3 = func(X[i] + dt/2. * k2, t[i] + dt/2.)
k4 = func(X[i] + dt * k3, t[i] + dt)
X[i+1] = X[i] + dt / 6. * (k1 + 2. * k2 + 2. * k3 + k4)
return X
# ODEint is preloaded.
# Define the time vector t and the initial conditions X0
Calculate and plot the kinetic energy Ec, the potential energy Ep and the total energy Et = Ec + Ep for all 3 cases, plot it and comment.