Lagrangian and Hamiltonian Systems

[jonniedie]

Based on this here's a first crack at how to handle Lagrangian systems.

using Symbolics, ModelingToolkit, DifferentialEquations

∂_∂(x) = y -> expand_derivatives(Differential(x)(y))

function LagrangianSystem(L, q̇, q, p, t;
                     Q = zeros(length(q)),
                     default_u0 = [q̇..., q...] .=> 0.0,
                     kwargs...)
    
    Q_vals = Q
    inds = 1:length(q)
    @variables q̈[inds](t) Q[inds](t)
    d_dt = ∂_∂(t)

    lagrange_eqs = map(q̇, q, Q) do q̇ᵢ, qᵢ, Qᵢ
        d_dt(∂_∂(q̇ᵢ)(L)) - ∂_∂(qᵢ)(L) ~ Qᵢ
    end

    try
        subs = d_dt.([q̇; q]) .=> [q̈; q̇]
        new_eqs = [substitute(eq.lhs, subs) ~ eq.rhs for eq in lagrange_eqs]
        lagrange_eqs = d_dt.(q̇) .~ Symbolics.solve_for(new_eqs, q̈)
    catch e
        @warn "Could not reduce equations due to: $e"
    end

    deriv_eqs = d_dt.(q) .~ q̇

    force_eqs = Q .~ Q_vals

    eqs = [lagrange_eqs; deriv_eqs; force_eqs]

    return structural_simplify(ODESystem(eqs, t, [q̇..., q..., Q...], p; default_u0, kwargs...))
end

When @YingboMa gets this stuff implemented, I can open up a PR along these lines.

[cadojo]

@jonniedie @YingboMa @ChrisRackauckas Are LagrangianSystem and HamiltonianSystem still desired in ModelingToolkit? I've been working on integrating galactic potential equations (in the integral / scalar form), and I created a new ScalarField subtype of AbstractTimeDependentSystem for this purpose. I think HamiltonianSystem would be more appropriate, and I'd be interested in helping to contribute this type if it's still desired in the SciML ecosystem.

[ChrisRackauckas]

I think it would be nice to have these alternative system forms, making it easier to generate the equations of motion. I'd be happy to review PRs around it but never got to it myself.