Solution (Workshop 2): Linearizing the Equations of Motion (Solution)#
In this tutorial you will learn to linearize the Equations of Motion of the pendulum as covered in the lecture.
# for time comparison later
import time
start_time = time.time()
Part 1: Kinematic equations#
We first start by defining the variables. Let’s start with the kinematic relations of a pendulum hanging from the point \((x_1, z_1)\). Note that the horizontal coordinate x is defined with a cosine, indicating that the angle is defined counter-clockwise starting from the horizontal. A hanging pendulum thus has \(\phi_1 = \frac{-\pi}{2} = \frac{3\pi}{2}\).
from sympy import *
var("t x1 z1 r phi0") # independent variables
phi1 = Function("phi1")(t)
x2 = x1 + r*cos(phi1)
z2 = z1 + r*sin(phi1)
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[2], line 1
----> 1 from sympy import *
2 var("t x1 z1 r phi0") # independent variables
3 phi1 = Function("phi1")(t)
ModuleNotFoundError: No module named 'sympy'
The velocities can then be obtained using:
xdot = diff(x2, t)
zdot = diff(z2, t)
Part 2: Energy equations#
Kinetic energy:#
The in this system is \(m\), such that the kinetic energy is given by:
var("m")
T = 0.5*m*(xdot**2 + zdot**2)
T.evalf()
This expression can be simplified to:
T = simplify(T) # Linearization
T.evalf()
Potential energy:#
The potential energy on the pendulum is due to gravity:
var("g")
V = m*g*z2
V.evalf()
Step 3: Construct the Lagrangian#
L = T - V
L.evalf()
Step 4: Obtaining the EoM#
In order to obtain the EoMs we have to take derivatives w.r.t. \(\phi_1\) and its velocity.
EOM_phi = diff( diff(L, diff(phi1, t)), t) - diff(L, phi1)
# ! LHS needs derivative to get "velocity"
EOM_phi = simplify(EOM_phi)
EOM_phi.evalf()
The equation is nonlinear since it depends on the cosine of the angle. We want to obtain a linear function, so let’s linearize. We consider a perturbation around the point we want to linearize (\(\phi_0\)). We linearize with the following function substitution: \(\phi_1(t) = \phi_0 + \epsilon \psi(t)\). The \(\epsilon\) is added in order to apply the Taylor series expansion.
On a technical note: Using SymPy for function substitution can be tricky, hence we will use temporary symbols called tmp1 and tmp2 in order to do the substitution.
var("phi0 epsilon")
psi = Function("psi")(t) # perturbation function
tmp1 = symbols("tmp1")
EOM_psi = EOM_phi.evalf(subs={phi1: tmp1})
EOM_psi = EOM_psi.evalf(subs={tmp1: phi0 + epsilon*psi})
print(EOM_psi)
EOM_psi.evalf()
m*r*(g*cos(epsilon*psi(t) + phi0) + 1.0*r*Derivative(epsilon*psi(t) + phi0, (t, 2)))
Now, we can apply the Taylor series expansion using the function \(\epsilon\) as a variable.
EOM_lin = series(EOM_psi, epsilon, n=2)
EOM_lin.evalf()
Note that we know that \(\frac{d^2\phi_0}{dt^2}=0\), so we can redefine the derivative in terms of \(\psi\) only.
tmp2 = symbols("tmp2")
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()
EOM_lin = series(EOM_psi2, epsilon, n=2)
EOM_lin.evalf()
Then, we obtain the linearized EOM by setting \(\epsilon = 1\).
EOM_lin = EOM_lin.removeO().evalf(subs={epsilon: 1})
EOM_lin.evalf()
We see that we get an expression that only depends on (linearly) on the perturbation \(\psi(t)\). Isolating the second derivative with respect to time of the perturbation, we get the final expression of the linearized system.
EOM_lin_iso = solve(EOM_lin, diff(psi, (t, 2)))
EOM_lin_iso[0].evalf()
In this problem, our initial angle is \(\phi_0 = 3*\frac{\pi}{2}\), then the final Equation of Motion will be:
a = EOM_lin_iso[0].evalf(subs={phi0: 3*pi/2})
a
end_time = time.time()
Part 3: Solve the equation#
Now we can solve the equation using an ODE solver
tsym, psisym = symbols("tsym psisym")
print(a)
a = a.evalf(subs={t: tsym, psi: psisym, g: 9.81, r: 1.0})
print(a)
def qdot(t,q):
vt = q[1]
at = a.evalf(subs={psisym: q[0], tsym: t})
return [vt,at]
-g*psi(t)/r
-9.81*psisym
# Just for demonstration
qdot(0,[0,0])
[0, 0]
from scipy.integrate import solve_ivp
sol = solve_ivp(fun=qdot,t_span=[0,10],y0=[1,0])
# Note: The initial "angle" of 1 would not be considered small, so the linearization would in practice not be physically accurate.
import matplotlib.pyplot as plt
plt.plot(sol.t,sol.y[0])
plt.xlabel("Time [s]")
plt.ylabel("Excursion [rad]")
plt.title("Pendulum motion from Lagrangian equations");
Part 4: Performance comparison#
print(end_time-start_time)
6.6967291831970215
Depending on the machine you are running on the execution of this notebook takes anywhere from 1 tin a few seconds. For the development of this tutorial a TU Delft issue HP Zbook G5 was used, carrying an i7-8750H CPU at 2.2 GHz with 16 GB of 64 bit RAM. The notebook took 1.64 seconds to execute.
An alternative way of doing symbolic calculations is through the use of Maplesoft. Below you can see an example of a Maple math file that performs the same calculations as in this notebook from the start_time to the end_time. This gives, as shown below, an execution time of 2.256. While Maple can be more intuitive for symbolic manipulations, it becomes clear that the computional cost of the sympy package is rather good.
Part 5: Exercise#
Now it is time to applied what you have learned. Using the tutorial above, try to construct, and then solve the following problem:
A mass-spring-damper system with a point mass of 200 kg at the end.
A spring with spring-stiffness 10 kN/m.
A damper with damping force 1 kN s/m.
An external force defined as (in kN) \(F(t) = 10 \sin(\pi t)\).
Hint: How would you implement the energy dissipation from the dashpot?
Anwer: Prefered way, in the Lagrangian derivation: \(\frac{d}{dt}\frac{dL}{d \dot{q}} - \frac{dL}{d \dot{q}} = \frac{dD}{d \dot{q}}\) with \(D = -0.5 c \dot{x}^2\) the dissipative energy.
# independent variables
var("t x0 x")
# position
x = Function("x")(t)
# velocity
xdot = diff(x,t)
# kinetic energy
var("m")
T = 0.5*m*xdot**2
# potential energy
var("k")
V = 0.5*k*(x-x0)**2
# lagrangian
L = T - V
# damping component
var("c")
D = 0.5*c*xdot**2
# obtaining the EOM, including the external force
Qi = var("Qi")
EOM_x = diff( diff(L, xdot), t) + diff(D, xdot) - diff(L, x) - Qi
# Does it make sense? Is it Linear?
simplify(EOM_x).evalf()
# extract the acceleration
xdotdot = solve(EOM_x, diff(x, (t,2)))
a = xdotdot[0].evalf()
# Solving the system
# Defining constants
m_model = 200
c_model = 1e3
k_model = 10e3
x0_model = 0
import numpy as np
def F_ext(t):
return 10e3*np.sin(np.pi * t)
t0 = 0
tend = 10
tsym, xsym, xdotsym = symbols("tsym xsym xdotsym")
a_model = a.evalf(subs={t: tsym, x: xsym, Derivative(x, t): xdotsym, m: m_model, c: c_model, k: k_model, x0: x0_model})
def qdot(t,q):
vt = q[1]
at = F_ext(t)/m_model + a_model.evalf(subs={xsym: q[0], tsym: t, xdotsym: q[1], Qi: F_ext(t)})
return [vt,at]
sol = solve_ivp(fun=qdot,t_span=[0,tend],y0=[x0_model,0])
plt.plot(sol.t,sol.y[0])
plt.xlabel("Time [s]")
plt.ylabel("Excursion [m]")
plt.title("MCK-system motion from Lagrangian equations");
