Solution (source code)

= Solution

For an <action functional> $S[q]=\int_{t_i}^{t_f}L(q,\dot q,t)dt$, vary a path while fixing its endpoint coordinates. <Integration by parts> gives
$$
\delta S=\left[\sum_i\frac{\partial L}{\partial\dot q_i}\delta q_i\right]_{t_i}^{t_f}
+\int_{t_i}^{t_f}\sum_i\left[\frac{\partial L}{\partial q_i}-\frac d{dt}\frac{\partial L}{\partial\dot q_i}\right]\delta q_i\,dt.
$$
The endpoint term vanishes. Since the interior variations are arbitrary, stationarity gives
$$
\boxed{\frac d{dt}\frac{\partial L}{\partial\dot q_i}-\frac{\partial L}{\partial q_i}=0.}
$$
Define <canonical momentum> $p_i=\partial L/\partial\dot q_i$. For a regular <velocity> Hessian, these relations can be inverted for $\dot q(q,p,t)$, and the <Legendre transform> defines
$$
\boxed{H(q,p,t)=\sum_ip_i\dot q_i-L(q,\dot q,t).}
$$
Its differential is $dH=\sum_i(\dot q_i\,dp_i-L_{q_i}\,dq_i)-L_tdt$. Inserting the <Euler-Lagrange equations> already suggests the canonical equations.

To derive them directly from an <principle of stationary action>, vary the first-order <action> $I[q,p]=\int(\sum_ip_i\dot q_i-H)dt$ with independent $q,p$ and fixed endpoint $q$. Integrating $p_i\delta\dot q_i$ by parts gives
$$
\delta I=\int\sum_i\left[(\dot q_i-H_{p_i})\delta p_i+(-\dot p_i-H_{q_i})\delta q_i\right]dt.
$$
Therefore
$$
\boxed{\dot q_i=\frac{\partial H}{\partial p_i},\qquad
\dot p_i=-\frac{\partial H}{\partial q_i}.}
$$
The endpoint <momenta> need not be fixed. A singular <velocity> Hessian instead requires a constrained <Hamiltonian> formulation; the ordinary Legendre inversion is an assumption, not an identity valid for every <Lagrangian>.