Lagrange’s Equations

$L = T - V$ $\frac{d}{dt}\left(\frac{\partial L}{\partial \dot{q_i}}\right) - \frac{\partial L}{\partial q_i} = Q_i$

Where:

• $$T$$ is the system’s kinetic energy
• $$V$$ is the system’s potential energy
• $$q_i$$ represents a coordinate in a system
• e.g in cartesian $$\mathbb{R}^3$$, we have $$(q_1, q_2, q_3) = (x, y, z)$$
• $$Q_i$$ are the forces acting on the $$i^{th}$$ coordinate

Example: Projectile

Kinetic Energy, T and Potential Energy, V

$T = \frac{1}{2}M\dot{x}^2 + \frac{1}{2}M\dot{y}^2$ $V = Mgy$

Lagrangian, L

$L = T - V$ $L = \frac{1}{2}M\dot{x}^2 + \frac{1}{2}M\dot{y}^2 - Mgy$

Forces

\begin{aligned} drag &= \frac{1}{2}\rho v^2 C_D A \\ &= \frac{1}{2}k_D v^2 \end{aligned}

Where:

• $$k_D = \rho C_D A$$ is approximately constant

Equations of motion w.r.t x

$L = \frac{1}{2}M\dot{x}^2 + \frac{1}{2}M\dot{y}^2 - Mgy$ $\frac{d}{dt}\left(\frac{\partial L}{\partial \dot{x}}\right) - \frac{\partial L}{\partial x} = F_x$ $\frac{\partial L}{\partial \dot{x}} = M\dot{x}, \frac{d}{dt}\left(\frac{\partial L}{\partial \dot{x}}\right) = M\ddot{x}, \frac{\partial L}{\partial x} = 0$ $F_x = \frac{1}{2}k_D \dot{x}^2$ $\therefore M\ddot{x} = \frac{1}{2}k_D \dot{x}^2$

Equations of motion w.r.t y

$L = \frac{1}{2}M\dot{x}^2 + \frac{1}{2}M\dot{y}^2 - Mgy$ $\frac{d}{dt}\left(\frac{\partial L}{\partial \dot{y}}\right) - \frac{\partial L}{\partial y} = F_y$ $\frac{\partial L}{\partial \dot{y}} = M\dot{y}, \frac{d}{dt}\left(\frac{\partial L}{\partial \dot{y}}\right) = M\ddot{y}, \frac{\partial L}{\partial y} = -Mg$ $F_y = \frac{1}{2}k_D \dot{y}^2$ $\therefore M\ddot{y} + Mg = \frac{1}{2}k_D \dot{y}^2$

State space representation

Rearranging gives:

\begin{aligned} \ddot{x} &= -\frac{1}{2M}k_D\dot{x} \\ \ddot{y} &= -g - \frac{1}{2M}k_D\dot{y} \end{aligned}

We can now represent this in an ode solver:

import scipy.integrate as integrate
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
from matplotlib.animation import FuncAnimation
import numpy as np

class Ball():
def __init__(self, Kd=0.0, M=1.0, g=9.81):
self.M = M
self.Kd = Kd
self.g = g

def _ode(self, X0, t):
[x, dx, y, dy] = X0
drag = self.Kd/(2.0 * self.M)

return [dx, -drag * dx, dy, -self.g - drag * dy]

def simulate(self, x_init, t):
return integrate.odeint(self._ode, x_init, t)

class DrawBall:
def __init__(self, ax):
with plt.style.context("dark_background"):
self.ball, = ax.plot([], [])
self.path, = ax.plot([], [], linestyle='--', alpha=0.7)
self.ax = ax
ax.grid(True, linestyle=":", alpha=0.5)

def draw(self, x, y):
phi = np.linspace(-np.pi, np.pi, 32)
self.ball.set_data(x + np.cos(phi), y + np.sin(phi))
xd, yd = self.path.get_data()
if x == 0.0 and y == 0.0:
xd, yd = [], []
self.path.set_data(np.append(xd, x), np.append(yd, y))

self.ax.set_xlim(x - 40, x + 40)
self.ax.set_ylim(y - 50, y + 50)
return self.ax,

ball = Ball(Kd=10.0, M=10.0)

T0 = 0.0
TN = 7.0
t = np.linspace(T0, TN, num=5000)

x_init = [0.0, 25.0, 0.0, 50.0]
x = ball.simulate(x_init, t)

with plt.style.context("dark_background"):
color1 = plt.rcParams['axes.prop_cycle'].by_key()['color'][0]
color2 = plt.rcParams['axes.prop_cycle'].by_key()['color'][1]
color3 = plt.rcParams['axes.prop_cycle'].by_key()['color'][2]
color4 = plt.rcParams['axes.prop_cycle'].by_key()['color'][3]
color5 = plt.rcParams['axes.prop_cycle'].by_key()['color'][4]

fig = plt.figure(figsize=(10, 5))
gs = GridSpec(2, 2, figure=fig)
ax0 = fig.add_subplot(gs[:, 0])
ax1 = fig.add_subplot(gs[0, 1])
ax11 = ax1.twinx()
ax2 = fig.add_subplot(gs[1, 1])
ax21 = ax2.twinx()

points = np.full(5, None)

ax1.plot(t, x[:, 0], color=color1)

ax1.set_ylabel(r'x (m)', color=color1)
ax11.plot(t, x[:, 1], color=color2)

ax11.set_ylabel(r'$\dot{x}$ ($\frac{m}{s}$)', color=color2)

ax2.plot(t, x[:, 2], color=color3)
ax2.set_ylabel(r'$\theta$ (rad)', color=color3)
ax21.plot(t, x[:, 3], color=color4)
ax21.set_ylabel(r'$\dot{\theta}$ ($\frac{rad}{s}$)', color=color4)

db = DrawBall(ax0)

points[0], = ax1.plot([], [], 'o', color=color1)
points[1], = ax11.plot([], [], 'o', color=color2)
points[2], = ax2.plot([], [], 'o', color=color3)
points[3], = ax21.plot([], [], 'o', color=color4)

def animate(i):
time = (i % 70) * 0.1
xp = np.interp(time, t, x[:, 0])
dxp = np.interp(time, t, x[:, 1])
yp = np.interp(time, t, x[:, 2])
dyp = np.interp(time, t, x[:, 3])

points[0].set_data(time, xp)
points[1].set_data(time, dxp)
points[2].set_data(time, yp)
points[3].set_data(time, dyp)

ax, = db.draw(xp, yp)

return ax, points[0], points[1], points[2], points[3], points[4],

ani = FuncAnimation(fig, animate, interval=100, save_count=70)
# ani.save('projectile.gif', writer='imagemagick', dpi=64)

plt.show()


Categories:

Updated: