Workshop 3: Deriving the EOM of a pendulum

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( $ϕ$ ): https://en.wikipedia.org/wiki/Rotation_matrix the zero angle points to the right for this matrix.

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 $M$

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 $l 0$ , the original length of the spring, does not change in time, then the elongation of the spring $d l$ is given by:

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 $k$ and gravity.

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. $ϕ_{1}$ and its velocity.

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 $ϕ (t) = 0$ was set to the right, rather than downward which is more conventional for a pendulum.

# 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 $k : 0$ (or another value) to one of the earlier dictionaries for example does the trick here, but it can only be applied if the goal was to not include a constant stiffness spring in the first place.

Bonus: Obtaining the EOMs for nDOF system#

If your system contains multiple DOFs, you will have to take the derivative of $L$ towards each of these separately, thereby obtaining the EOM for each DOF. Let’s say you obtained the following two EOMs (the EOMs are just random things I entered):

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()

The solution can be found here.