Lagrangian and Hamiltonian mechanics

Lagrangian and Hamiltonian mechanics#

Added in version 4.0.0.

heyoka.py offers direct support for Lagrangian and Hamiltonian mechanics. This means that you can supply heyoka.py a Lagrangian or Hamiltonian, and heyoka.py will take care of formulating the equations of motion for you.

Let us explore this feature with a couple of simple examples.

Pendulum on a movable support#

Consider the motion of a simple pendulum attached to a massive support that can freely move horizontally:

../_images/movable_pendulum.svg

(credits to Wikipedia for the image)

We can choose \(x\) and \(\theta\) as the generalised coordinates for this system. Let us introduce a few initial quantities:

import heyoka as hy

# Generalised coordinates.
x, th = hy.make_vars("x", "theta")

# Generalised velocities.
vx, vth = hy.make_vars("vx", "vtheta")

# System parameters.
M = hy.par[0]
m = hy.par[1]
l = hy.par[2]
g = hy.par[3]

The masses \(M\) and \(m\), the pendulum length \(l\) and the gravitational acceleration \(g\) have all been introduced as runtime parameters.

The Cartesian position of the pendulum’s bob in terms of the generalised coordinates can be written as:

\[\begin{split} \begin{cases} x_{\mathrm{pend}} & = x + l\sin\left( \theta \right)\\ y_{\mathrm{pend}} & = -l\cos\left( \theta \right) \end{cases}\,. \end{split}\]

In order to compute the kinetic energy of the system, we must first compute the cartesian velocity of the bob as a function of the generalised coordinates and velocities. We do not have to do this by hand, rather we can let heyoka.py’s symbolic differentiation capabilities do the work for us (see the tutorial on computing derivatives):

import numpy as np

# Cartesian coordinates of the bob.
xpend = x + l * hy.sin(th)
ypend = -l * hy.cos(th)

# Cartesian velocities of the bob.
vxpend, vypend = hy.diff_tensors([xpend, ypend], diff_args=[x, th]).jacobian @ [vx, vth]

We can now formulate the kinetic energy, which reads:

\[ T = \frac{1}{2}M\dot{x}^2 + \frac{1}{2}m\left( \dot{x}_\mathrm{pend}^2 + \dot{y}_\mathrm{pend}^2 \right). \]
T = 0.5 * M * vx**2 + 0.5 * m * (vxpend**2 + vypend**2)

The gravitational potential of the bob instead reads:

\[ V=mgy_\mathrm{pend}. \]
V = m * g * ypend

We can now assemble the Lagrangian \(L = T - V\):

L = T - V

The equations of motion can now be formulated via the lagrangian() function. This function takes 3 mandatory arguments in input:

  • the Lagrangian,

  • the list of generalised coordinates,

  • the list of generalised velocities.

lagrangian() will return a list of differential equations that can be directly used to create a taylor_adaptive() integrator.

# Generate the system of differential equations.
sys = hy.lagrangian(L, [x, th], [vx, vth])

We can now assign values to the parameters of the system, and create an integrator object:

# Parameter values.
Mval, mval, lval, gval = [1.0, 0.1, 0.7, 0.8]

# Integrator object.
ta = hy.taylor_adaptive(sys, [0.0, 0.3, 0.0, 0.0], pars=[Mval, mval, lval, gval])

Let us numerically integrate the system for a few time units:

res = ta.propagate_grid(np.linspace(0, 50, 200))

We are now ready to produce a nice animation:

%%capture
%matplotlib inline

import matplotlib.pyplot as plt
from matplotlib import animation, rc
from IPython.display import HTML

hist = res[-1]

fig = plt.figure(figsize=(10, 10))
ax = plt.subplot(111)

# Init the graphical elements.
pend = plt.Circle((0.2, 0.0), 0.03, ec="black", fc="black", zorder=3)
supp_width = 0.4
supp = plt.Rectangle((0.1, 0.0), supp_width, 0.2)
(ln_pend,) = ax.plot([], [], "k--")
y_axis = plt.hlines(0, -2, 2, colors="gray", linewidths=4.0)

ax.add_artist(pend)
ax.add_artist(supp)


def init():
    ax.set_xlim((-0.4, 0.4))
    ax.set_ylim((-0.8, 0.4))
    ax.set_aspect("equal")
    return (pend, supp)


def animate(i):
    cur_x, cur_th = hist[i][0:2]
    pend_x, pend_y = cur_x + lval * np.sin(cur_th), -lval * np.cos(cur_th)
    supp_x = cur_x - supp_width / 2

    pend.set_center((pend_x, pend_y))

    supp.set_xy((supp_x, 0.0))

    ln_pend.set_data([pend_x, cur_x], [pend_y, 0.0])

    return (pend, supp)


anim = HTML(
    animation.FuncAnimation(
        fig, animate, init_func=init, frames=200, interval=50, blit=True
    ).to_jshtml()
)

anim