Getting Started

This notebook demonstrates the generic workflow in the Lagrangian machinery for holonomic systems.

• Decide upon your configuration space: define the Lagrangian coordinates \(\mathbf q\)

• Write the Kinetic energy \(T(\mathbf q, \dot{\mathbf q})\)

• Write the Potential energy \(V(\mathbf q)\)

• Assemble the Lagrangian \(L = T - V\)

• Define the holonomic constraints \(\mathbf F(\mathbf q) = 0\)

• Write the equations of motion using quivira : lagrange_eom(L, F, q, dq)

The example is a bead constrained to move along a circular wire, subject to gravity.

# Our stuff
import heyoka as hy
import quivira as qv

# Third party
import numpy as np

# Plots and visualization
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from IPython.display import HTML

1. We decide upon our configuration space and define the Lagrangian coordinates \(\mathbf q = [x, y]\) as the Cartesian positions of the bead.

# Generalized coordinates: the Cartesian position of the bead.
x, y = hy.make_vars("x" ,"y")
vx, vy = hy.make_vars("vx", "vy")
q = [x, y]
dq = [vx, vy]

2. We write the kinetic energy of the bead as \(T = \frac{1}{2} m (\dot{x}^2 + \dot{y}^2)\).

m = 1. # kg
T = m * (vx**2 + vy**2) / 2

3. We write the Potential energy of the bead as \(V = m g y\).

g = 9.80665 # m/s^2
V = m * g * y

4. Assemble the Lagrangian of the bead as \(L = T - V\).

L = T - V

5. We define the holonomic constraint of the bead as \(x^2 + y^2 = 1.\).

F = [x**2 + y**2 - 1]

6. We assemble the equations of motion for the bead using Quivira amazing capabilities.

eom = qv.lagrange_eom(L, F, q, dq)

We can now use the assembled equations of motion to instantiate heyoka Taylor adaptive integrator

ta = hy.taylor_adaptive(eom)

ta.time = 0.
ta.state[:] = [1.,0.,0.,0.]
tgrid = np.linspace(0.,2.,200)
sol = ta.propagate_grid(tgrid)

states = sol[-1]
bead_x = states[:, 0]
bead_y = states[:, 1]

wire_angle = np.linspace(0.0, 2.0 * np.pi, 300)
wire_x = np.cos(wire_angle)
wire_y = np.sin(wire_angle)

fig, ax = plt.subplots(figsize=(6, 6))
ax.plot(wire_x, wire_y, color="0.35", linewidth=2, label="Wire")
bead, = ax.plot([], [], "o", color="tab:orange", markersize=10, label="Bead")
ax.set_aspect("equal")
ax.set_xlim(-1.2, 1.2)
ax.set_ylim(-1.2, 1.2)
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.grid(True)
ax.legend(loc="upper right")

def animate(frame):
    bead.set_data([bead_x[frame]], [bead_y[frame]])
    return (bead,)


animation = FuncAnimation(
    fig,
    animate,
    frames=len(states),
    interval=30,
    blit=True,
)
display(HTML(animation.to_jshtml()))
plt.close(fig)

Once Loop Reflect

Surely there were easier ways to simulate the bead on the wire, like choosing a single angle as Lagrangian variable and no constraints .. but weel it was fun to do it this way.