Numerical solution for Hamilton's equations

[YingboMa]

One can compute Hamiltonian from Lagrangian via Legendre transformation, so I propose the following implementation.

function Legendre_transformation(F, w) # Legendre transformation for quadratic function
    wv = a->ForwardDiff.gradient(F, a)
    z = zeros(w)
    M = ForwardDiff.jacobian(wv, z)
    b = wv(z)
    v = bkfact!(M)\(w-b)
    w'v - F(v)
end

function Lagrangian2Hamiltonian(Lagrangian, t, q, p)
    L = q̇ -> Lagrangian(t, q, q̇)
    Legendre_transformation(L, p)
end

Finding the Legendre transformation of an arbitrary function is very hard numerically, but Lagrangian is quadratic in most of mechanics.