Degenerate Lagrangian Systems

A Lagrangian $L(x,v)$ is said to be degenerate if the determinant of the Hessian with respect to the velocities vanishes, i.e.,

\[\left\vert \frac{\partial^2 L}{\partial v^i \, \partial v^j} \right\vert = 0 .\]

The Euler-Lagrange equations,

\[\frac{d}{dt} \frac{\partial L}{\partial v} - \frac{\partial L}{\partial x} = 0 ,\]

of such Lagrangians are a set of first-order ordinary differential equations and not second-order differential equations as for non-degenerate Lagrangians.

Of particular interest in various applications are degenerate Lagrangians of the form

\[L(x,v) = \vartheta(x) \cdot v - H (x) .\]

Their Euler-Lagrange equations take the form

\[\frac{d \vartheta}{dt} (x) = \nabla \vartheta(x) \cdot v - \nabla H (x) .\]

We exemplify this with the Lotka-Volterra problem in 2d.

Lotka-Volterra

Before any use, we need to load EulerLagrange:

using EulerLagrange

Next, we generate symbolic variables for a two-dimensional Lagrangian system:

t, x, v = lagrangian_variables(2)
(t, (x(t))[1:2], (v(t))[1:2])

We define a named tuple with typical values for the parameters, e.g.,

params = (
    a₁ = -1.0,
    a₂ = -1.0,
    b₁ = 1.0,
    b₂ = 2.0,
)
(a₁ = -1.0, a₂ = -1.0, b₁ = 1.0, b₂ = 2.0)

We use the function symbolize to generate a symbolic version of the parameters:

sparams = symbolize(params)
(a₁ = a₁ₚ, a₂ = a₂ₚ, b₁ = b₁ₚ, b₂ = b₂ₚ)

Define the Hamiltonian function and the symplectic potential:

H(x, params) = params.a₁ * x[1] + params.a₂ * x[2] + params.b₁ * log(x[1]) + params.b₂ * log(x[2])
ϑ(x, params) = [log(x[2]) / x[1] / 2, - log(x[1]) / x[2] / 2]
ϑ (generic function with 1 method)

The Hamiltonian and the symplectic potential, evaluated on and together with the symbolic variables and parameters are used to construct a DegenerateLagrangianSystem:

lag_sys = DegenerateLagrangianSystem(ϑ(x,sparams), H(x,sparams), t, x, v, sparams)

Degenerate Lagrangian system with

L = ((1//2)*(v(t))[1]*log((x(t))[2])) / (x(t))[1] + ((-1//2)*(v(t))[2]*log((x(t))[1])) / (x(t))[2] - a₁ₚ*(x(t))[1] - a₂ₚ*(x(t))[2] - b₁ₚ*log((x(t))[1]) - b₂ₚ*log((x(t))[2])

and equations of motion

-a₁ₚ + (-(v(t))[2]) / (2(x(t))[1]*(x(t))[2]) + ((-1//2)*(v(t))[1]*log((x(t))[2])) / ((x(t))[1]^2) + (-b₁ₚ) / (x(t))[1] + (-Differential(t)((x(t))[2])) / (2(x(t))[1]*(x(t))[2]) + (((1//2)*log((x(t))[2])) / ((x(t))[1]^2))*Differential(t)((x(t))[1])
-a₂ₚ + Differential(t)((x(t))[1]) / (2(x(t))[1]*(x(t))[2]) + ((1//2)*(v(t))[2]*log((x(t))[1])) / ((x(t))[2]^2) + (-b₂ₚ) / (x(t))[2] + (v(t))[1] / (2(x(t))[1]*(x(t))[2]) + (((-1//2)*log((x(t))[1])) / ((x(t))[2]^2))*Differential(t)((x(t))[2])

The constructor computes the Euler-Lagrange equations and generates the corresponding Julia code.

In the last step, we can now construct a LODEProblem from the LagrangianSystem and some appropriate initial conditions, a time span to integrate over and a time step:

tspan = (0.0, 10.0)
tstep = 0.01

q₀ = [2.0, 1.0]
p₀ = ϑ(q₀, params)

lprob = LODEProblem(lag_sys, tspan, tstep, q₀, p₀; parameters = params)
Geometric Equation Problem for Lagrangian Ordinary Differential Equation (LODE)

 with vector fields
   ϑ = RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :t, :X, :V, :params), EulerLagrange.var"#_RGF_ModTag", EulerLagrange.var"#_RGF_ModTag", (0x7e3af40f, 0xf7247472, 0xe3c9590c, 0x32f9ee27, 0xa81cb7c6), Expr}(quote
    #= none:1 =#
    #= none:5 =#
    begin
        #= none:7 =#
        #= none:7 =# @inbounds begin
                #= none:9 =#
                ˍ₋out[1] = (1 // 2 * log(getindex(X, 2))) / getindex(X, 1)
                #= none:10 =#
                ˍ₋out[2] = (-1 // 2 * log(getindex(X, 1))) / getindex(X, 2)
                #= none:12 =#
                nothing
            end
    end
end)
   f = RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :t, :X, :V, :params), EulerLagrange.var"#_RGF_ModTag", EulerLagrange.var"#_RGF_ModTag", (0x601e816e, 0x5f412bb9, 0xa6618f59, 0xd778f2ee, 0xdee03567), Expr}(quote
    #= none:1 =#
    #= none:5 =#
    begin
        #= none:7 =#
        #= none:7 =# @inbounds begin
                #= none:9 =#
                ˍ₋out[1] = ((-1 * params.a₁ + (-1 * params.b₁) / getindex(X, 1)) + ((-1 // 2 * getindex(V, 1)) * log(getindex(X, 2))) / getindex(X, 1) ^ 2) + (-1 * getindex(V, 2)) / ((2 * getindex(X, 1)) * getindex(X, 2))
                #= none:10 =#
                ˍ₋out[2] = ((-1 * params.a₂ + (-1 * params.b₂) / getindex(X, 2)) + ((1 // 2 * getindex(V, 2)) * log(getindex(X, 1))) / getindex(X, 2) ^ 2) + getindex(V, 1) / ((2 * getindex(X, 1)) * getindex(X, 2))
                #= none:12 =#
                nothing
            end
    end
end)
   g = RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :t, :X, :Λ, :V, :params), EulerLagrange.var"#_RGF_ModTag", EulerLagrange.var"#_RGF_ModTag", (0xf4cd6db1, 0xc8ad859a, 0x83e1263f, 0x3d20d3c7, 0x2ebbb721), Expr}(quote
    #= none:1 =#
    #= none:5 =#
    begin
        #= none:7 =#
        #= none:7 =# @inbounds begin
                #= none:10 =#
                nothing
            end
        #= none:12 =#
        begin
            #= none:14 =#
            #= none:14 =# @inbounds begin
                    #= none:16 =#
                    (ˍ₋out[1])[1] = (Differential(t))((1 // 2 * log(getindex(X, 2))) / getindex(X, 1))
                    #= none:17 =#
                    (ˍ₋out[1])[2] = (Differential(t))((-1 // 2 * log(getindex(X, 1))) / getindex(X, 2))
                    #= none:19 =#
                    nothing
                end
        end
    end
end)

 Lagrangian: L = RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:t, :X, :V, :params), EulerLagrange.var"#_RGF_ModTag", EulerLagrange.var"#_RGF_ModTag", (0xe192e52e, 0xa2be6a58, 0xed9a779a, 0x92db5d3e, 0x526e1366), Expr}(quote
    #= none:1 =#
    #= none:5 =#
    ((((((1 // 2 * getindex(V, 1)) * log(getindex(X, 2))) / getindex(X, 1) + ((-1 // 2 * getindex(V, 2)) * log(getindex(X, 1))) / getindex(X, 2)) + (-1 * getindex(X, 1)) * params.a₁) + (-1 * getindex(X, 2)) * params.a₂) + (-1 * params.b₁) * log(getindex(X, 1))) + (-1 * params.b₂) * log(getindex(X, 2))
end)

 Invariants: 
   GeometricBase.NullInvariants()

 Timespan: (0.0, 10.0) 
 Timestep: 0.01 

 Initial conditions: 
   (t = 0.0, q = [2.0, 1.0], p = [0.0, -0.34657359027997264], λ = [0.0, 0.0])

 Parameters: 
   (a₁ = -1.0, a₂ = -1.0, b₁ = 1.0, b₂ = 2.0)

We can integrate this system using GeometricIntegrators:

using GeometricIntegrators
sol = integrate(lprob, Gauss(1))

using CairoMakie
fig = lines(parent(sol.q[:,1]), parent(sol.q[:,2]);
    axis = (; xlabel = "x₁", ylabel = "x₂", title = "Lotka-Volterra system in 2d"),
    figure = (; size = (800,600), fontsize = 22))