Workshop 3: Deriving the EOM of a pendulum#
In this tutorial you will learn to derive the Equations of Motion of the pendulum as covered in the lecture. That is, a pendulum attached to the floor with a spring. We also consider an external horizontal force acting on the mass.

Part 1: Kinematic equations#
Using the following rotation matrix R(
We first start by defining the variables. Assuming that x1, z1 and r1 do not change in time:
import numpy as np
from sympy import *
var("t s x1 z1 r1")
phi1 = Function("phi1")(t)
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[1], line 2
1 import numpy as np
----> 2 from sympy import *
4 var("t s x1 z1 r1")
5 phi1 = Function("phi1")(t)
ModuleNotFoundError: No module named 'sympy'
Problem: Define the kinematic relations
Hint: How do you write the position of the pendulum (x2, z2) as a function of x1, z1, r1 and phi1?
# Define the kinematic relations here
The velocities can then be obtained using:
Problem: Define the velocities
Hint: Use derivatives.
# Compute/define the velocities here
Part 2: Energy equations#
Kinetic energy:#
The only mass in this system is the point mass located at P2 with mass
Problem: Define the kinetic energy
var("M")
# Define the kinetic energy here (T)
This expression can be simplified to:
T = simplify(T)
T.evalf()
Potential energy:#
Assuming again that
Problem: Define the spring elongation dl
var("dx dz l0 x3 z3")
# Define the spring elongation here dl(l0,x2,z2,x3,z3)
The work done by the spring between P2 and P3 with stiffness
Problem: Define the potential energy V.
var("k g")
# Define the potential energy here (V)
Work by external force#
The work done by an external force working in the horizontal direction:
Problem: Define the external work.
Fx = Function("Fx")(t)
# Define your external work here (W)
Step 3: Construct the Lagrangian#
Problem: Use T, V and W to find the lagrangian L. Simplify it using L.evalf()
# Define your Lagrangian here (L)
Step 4: Obtaining the EoM#
In order to obtain the EoMs we have to take derivatives w.r.t.
Problem: Use the Lagrangian expression to obtain the equation of motion.
Hint: Put all terms on the left hand side
# Compute the EOM here
Now we isolate it for the acceleration.
These cells can in practice take quite some time to execute. For this illustration the added spring stiffness k is set to 0, but it can take any value.
var("phi0 epsilon")
psi = Function("psi")(t) # perturbation function
tmp1 = symbols("tmp1") # the acceleration, to be linearized as phi0+epsilon*psi
tmp2 = symbols("tmp2") # the displacement
EOM_psi2 = EOM_phi.evalf(subs={diff(phi1, (t, 2)): tmp2, phi1: tmp1})
EOM_psi2 = EOM_psi2.evalf(subs={tmp2: diff(phi0 + epsilon*psi, (t, 2)), tmp1: phi0 + epsilon*psi})
EOM_psi2.evalf()
# Linearize by doing series expansion up to order 2 (see last term of output)
EOM_lin = series(EOM_psi2, epsilon, n=2)
EOM_lin.evalf()
EOM_lin2 = EOM_lin.evalf(subs={k: 0}) # For faster execution
EOM_lin2.evalf()
# Simplify and say that the "epsilon" was just a temorary value for perturbations
EOM_lin3 = EOM_lin2.removeO().evalf(subs={epsilon: 1})
EOM_lin3.evalf()
# Isolate the acceleration
EOM_lin_iso = solve(EOM_lin3, diff(psi, (t, 2)))
EOM_lin_iso[0].evalf()
The EOM of an actual pendulum can be recovered by removing the spring and the external force. Note that this returns a cos because
# Setting k=0 (done before) and phi0 straight down, you should recover the pendulum equation found in tutorial 2_1
EOM_lin_iso2 = EOM_lin_iso[0].evalf(subs={phi0: 3*pi/2, k: 0})
print(EOM_lin_iso2)
In the code above the high length of the expressions is a consequence of machine precision. This makes the execution times a multitude larger as well. Depending on the application it becomes more and more usefulf to simplify as early as possible in hte process. adding
Bonus: Obtaining the EOMs for nDOF system#
If your system contains multiple DOFs, you will have to take the derivative of
var("m c k J d q")
u = Function("u")(t)
w = Function("w")(t)
eom1 = m*cos(w)*diff(u, (t, 2)) + m*sin(u)*diff(w, (t, 2)) - c*diff(w, t) - c*diff(w, t)**2*u*ln(sqrt(t))**2*exp(u*w) - k*u
eom1.evalf()
eom2 = J*w**2*diff(u, (t, 2)) + J*diff(w, t)*diff(w, (t, 2)) - d*diff(w, t) + q*u
eom2.evalf()
You can collect both EOMs into matrix form using the following code:
# linear_eq_to_matrix only accepts symbols, therefore we have to substitute the second derivative functions with symbols a1 and a2
a1, a2 = symbols("a1 a2") # acceler u and acceler w
eom1 = eom1.evalf(subs={diff(u, (t, 2)): a1, diff(w, (t, 2)): a2})
eom2 = eom2.evalf(subs={diff(u, (t, 2)): a1, diff(w, (t, 2)): a2})
MTRX = linear_eq_to_matrix([eom1, eom2], [a1, a2])
Resulting in the following system of equations:
M = MTRX[0]
M.evalf()
And the following right-hand-side:
F = MTRX[1]
F.evalf()