Lagrangian Systems
The Euler-Lagrange equations, that is the dynamical equations of a Lagrangian system, are given in terms of the Lagrangian $L(x,v)$ by
\[\frac{d}{dt} \frac{\partial L}{\partial v} - \frac{\partial L}{\partial x} = 0 .\]
For regular (i.e. non-degenerate) Lagrangians, this is a set of second-order ordinary differential equations. In many numerical applications, it is advantageous to solve the implicit form of these equations, given by
\[\begin{align*} \frac{d \vartheta}{dt} &= f , & \vartheta &= \frac{\partial L}{\partial v} , & f = \frac{\partial L}{\partial x} . \end{align*}\]
In the following, we show how these equations can be obtained for the example of a particle in a square potential.
Particle in a potential
Before any use, we need to load EulerLagrange:
using EulerLagrangeNext, we generate symbolic variables for a two-dimensional system:
t, x, v = lagrangian_variables(2)(t, (x(t))[1:2], (v(t))[1:2])With those variables, we can construct a Lagrangian
using LinearAlgebra
L = v ⋅ v / 2 - x ⋅ x / 2-(1//2)*LinearAlgebra.dot(x(t), x(t)) + LinearAlgebra.dot(v(t), v(t)) / 2This Lagrangian together with the symbolic variables is then used to construct a LagrangianSystem:
lag_sys = LagrangianSystem(L, t, x, v)
Lagrangian system with
L = ((v(t))[1]^2 + (v(t))[2]^2) / 2 - (1//2)*((x(t))[1]^2 + (x(t))[2]^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₀ = [1.0, 1.0]
p₀ = [0.5, 2.0]
lprob = LODEProblem(lag_sys, tspan, tstep, q₀, p₀)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", (0x46371b34, 0xde7758bd, 0x645c9d07, 0xc0897de9, 0x38bd4dc0), Expr}(:(#= /home/runner/.julia/packages/Symbolics/KL0Tx/src/build_function.jl:456 =# @inbounds begin
#= /home/runner/.julia/packages/Symbolics/KL0Tx/src/build_function.jl:456 =#
begin
#= /home/runner/.julia/packages/SymbolicUtils/HQQQT/src/code.jl:1298 =# @inbounds begin
#= /home/runner/.julia/packages/SymbolicUtils/HQQQT/src/code.jl:1294 =#
ˍ₋out[1] = (getindex)(V, 1)
ˍ₋out[2] = (getindex)(V, 2)
#= /home/runner/.julia/packages/SymbolicUtils/HQQQT/src/code.jl:1296 =#
ˍ₋out
end
end
end))
f = RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :t, :X, :V, :params), EulerLagrange.var"#_RGF_ModTag", EulerLagrange.var"#_RGF_ModTag", (0xe4adc2d4, 0x2f566ab1, 0xbefe602c, 0xde2ceef2, 0x74e3da50), Expr}(:(#= /home/runner/.julia/packages/Symbolics/KL0Tx/src/build_function.jl:456 =# @inbounds begin
#= /home/runner/.julia/packages/Symbolics/KL0Tx/src/build_function.jl:456 =#
begin
begin
var"##cse#1" = -1//1
var"##cse#2" = (*)(var"##cse#1", (getindex)(X, 1))
var"##cse#3" = (*)(var"##cse#1", (getindex)(X, 2))
#= /home/runner/.julia/packages/SymbolicUtils/HQQQT/src/code.jl:1298 =# @inbounds begin
#= /home/runner/.julia/packages/SymbolicUtils/HQQQT/src/code.jl:1294 =#
ˍ₋out[1] = var"##cse#2"
ˍ₋out[2] = var"##cse#3"
#= /home/runner/.julia/packages/SymbolicUtils/HQQQT/src/code.jl:1296 =#
ˍ₋out
end
end
end
end))
g = RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :t, :X, :V, :Λ, :params), EulerLagrange.var"#_RGF_ModTag", EulerLagrange.var"#_RGF_ModTag", (0xb6b215b1, 0xd0bac697, 0x36895769, 0x652dc15f, 0x64f3d2c0), Expr}(:(#= /home/runner/.julia/packages/Symbolics/KL0Tx/src/build_function.jl:456 =# @inbounds begin
#= /home/runner/.julia/packages/Symbolics/KL0Tx/src/build_function.jl:456 =#
begin
#= /home/runner/.julia/packages/SymbolicUtils/HQQQT/src/code.jl:1298 =# @inbounds begin
#= /home/runner/.julia/packages/SymbolicUtils/HQQQT/src/code.jl:1294 =#
ˍ₋out[1] = (getindex)(Λ, 1)
ˍ₋out[2] = (getindex)(Λ, 2)
#= /home/runner/.julia/packages/SymbolicUtils/HQQQT/src/code.jl:1296 =#
ˍ₋out
end
end
end))
Lagrangian: L = RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:t, :X, :V, :params), EulerLagrange.var"#_RGF_ModTag", EulerLagrange.var"#_RGF_ModTag", (0x7aca55e8, 0xd7ab62da, 0x54cf5451, 0x09cbc704, 0x4a25bb77), Expr}(:(#= /home/runner/.julia/packages/Symbolics/KL0Tx/src/build_function.jl:198 =# @inbounds begin
#= /home/runner/.julia/packages/Symbolics/KL0Tx/src/build_function.jl:198 =#
begin
begin
var"##cse#1" = -1//2
var"##cse#2" = 2
var"##cse#3" = (^)((getindex)(X, 1), var"##cse#2")
var"##cse#4" = (^)((getindex)(X, 2), var"##cse#2")
var"##cse#5" = (+)(var"##cse#3", var"##cse#4")
var"##cse#6" = (*)(var"##cse#1", var"##cse#5")
var"##cse#7" = (^)((getindex)(V, 1), var"##cse#2")
var"##cse#8" = (^)((getindex)(V, 2), var"##cse#2")
var"##cse#9" = (+)(var"##cse#7", var"##cse#8")
var"##cse#10" = (/)(var"##cse#9", var"##cse#2")
var"##cse#11" = (+)(var"##cse#6", var"##cse#10")
var"##cse#11"
end
end
end))
Invariants:
GeometricBase.NullInvariants()
Timespan: (0.0, 10.0)
Timestep: 0.01
Initial conditions:
(t = fill(0.0), q = [1.0, 1.0], p = [0.5, 2.0], v = [0.0, 0.0])
Parameters:
GeometricBase.NullParameters()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 = "Particle moving in a square potential"),
figure = (; size = (800,600), fontsize = 22))┌ Warning: Hermite Extrapolation: q's history[1] and history[2] are identical!
└ @ GeometricIntegratorsBase ~/.julia/packages/GeometricIntegratorsBase/bVkej/src/extrapolation/hermite.jl:206
Parameters
We can also include parametric dependencies in the Lagrangian. Consider, for example, a parameter α that determines the strength of the potential.
The easiest way, to account for parameters, is to create a named tuple with typical values for each parameter, e.g.,
params = (α = 5.0,)(α = 5.0,)In the next step, we use the function symbolize to generate a symbolic version of the parameters:
sparams = symbolize(params)(α = αₚ,)Now we modify the Lagrangian to account for the parameter:
L = v ⋅ v / 2 - sparams.α * (x ⋅ x) / 2LinearAlgebra.dot(v(t), v(t)) / 2 - (1//2)*LinearAlgebra.dot(x(t), x(t))*αₚFrom here on, everything follows along the same lines as before, the only difference being that we also need to pass the symbolic parameters sparams to the LagrangianSystem constructor:
lag_sys = LagrangianSystem(L, t, x, v, sparams)
Lagrangian system with
L = ((v(t))[1]^2 + (v(t))[2]^2) / 2 - (1//2)*((x(t))[1]^2 + (x(t))[2]^2)*αₚAnalogously, we need to pass actual parameter values params to the LODEProblem constructor via the parameters keyword argument:
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", (0x46371b34, 0xde7758bd, 0x645c9d07, 0xc0897de9, 0x38bd4dc0), Expr}(:(#= /home/runner/.julia/packages/Symbolics/KL0Tx/src/build_function.jl:456 =# @inbounds begin
#= /home/runner/.julia/packages/Symbolics/KL0Tx/src/build_function.jl:456 =#
begin
#= /home/runner/.julia/packages/SymbolicUtils/HQQQT/src/code.jl:1298 =# @inbounds begin
#= /home/runner/.julia/packages/SymbolicUtils/HQQQT/src/code.jl:1294 =#
ˍ₋out[1] = (getindex)(V, 1)
ˍ₋out[2] = (getindex)(V, 2)
#= /home/runner/.julia/packages/SymbolicUtils/HQQQT/src/code.jl:1296 =#
ˍ₋out
end
end
end))
f = RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :t, :X, :V, :params), EulerLagrange.var"#_RGF_ModTag", EulerLagrange.var"#_RGF_ModTag", (0x1a0a201d, 0x803dccf8, 0xd1987318, 0x0e3a59ca, 0xb87045fe), Expr}(:(#= /home/runner/.julia/packages/Symbolics/KL0Tx/src/build_function.jl:456 =# @inbounds begin
#= /home/runner/.julia/packages/Symbolics/KL0Tx/src/build_function.jl:456 =#
begin
begin
var"##cse#1" = -1//1
var"##cse#2" = (*)((*)(var"##cse#1", (getindex)(X, 1)), params.α)
var"##cse#3" = (*)((*)(var"##cse#1", (getindex)(X, 2)), params.α)
#= /home/runner/.julia/packages/SymbolicUtils/HQQQT/src/code.jl:1298 =# @inbounds begin
#= /home/runner/.julia/packages/SymbolicUtils/HQQQT/src/code.jl:1294 =#
ˍ₋out[1] = var"##cse#2"
ˍ₋out[2] = var"##cse#3"
#= /home/runner/.julia/packages/SymbolicUtils/HQQQT/src/code.jl:1296 =#
ˍ₋out
end
end
end
end))
g = RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :t, :X, :V, :Λ, :params), EulerLagrange.var"#_RGF_ModTag", EulerLagrange.var"#_RGF_ModTag", (0xb6b215b1, 0xd0bac697, 0x36895769, 0x652dc15f, 0x64f3d2c0), Expr}(:(#= /home/runner/.julia/packages/Symbolics/KL0Tx/src/build_function.jl:456 =# @inbounds begin
#= /home/runner/.julia/packages/Symbolics/KL0Tx/src/build_function.jl:456 =#
begin
#= /home/runner/.julia/packages/SymbolicUtils/HQQQT/src/code.jl:1298 =# @inbounds begin
#= /home/runner/.julia/packages/SymbolicUtils/HQQQT/src/code.jl:1294 =#
ˍ₋out[1] = (getindex)(Λ, 1)
ˍ₋out[2] = (getindex)(Λ, 2)
#= /home/runner/.julia/packages/SymbolicUtils/HQQQT/src/code.jl:1296 =#
ˍ₋out
end
end
end))
Lagrangian: L = RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:t, :X, :V, :params), EulerLagrange.var"#_RGF_ModTag", EulerLagrange.var"#_RGF_ModTag", (0xa18ed549, 0x431ea4a1, 0x243771ae, 0x18e23afe, 0x98669bcf), Expr}(:(#= /home/runner/.julia/packages/Symbolics/KL0Tx/src/build_function.jl:198 =# @inbounds begin
#= /home/runner/.julia/packages/Symbolics/KL0Tx/src/build_function.jl:198 =#
begin
begin
var"##cse#1" = -1//2
var"##cse#2" = 2
var"##cse#3" = (^)((getindex)(X, 1), var"##cse#2")
var"##cse#4" = (^)((getindex)(X, 2), var"##cse#2")
var"##cse#5" = (+)(var"##cse#3", var"##cse#4")
var"##cse#6" = (*)((*)(var"##cse#1", var"##cse#5"), params.α)
var"##cse#7" = (^)((getindex)(V, 1), var"##cse#2")
var"##cse#8" = (^)((getindex)(V, 2), var"##cse#2")
var"##cse#9" = (+)(var"##cse#7", var"##cse#8")
var"##cse#10" = (/)(var"##cse#9", var"##cse#2")
var"##cse#11" = (+)(var"##cse#6", var"##cse#10")
var"##cse#11"
end
end
end))
Invariants:
GeometricBase.NullInvariants()
Timespan: (0.0, 10.0)
Timestep: 0.01
Initial conditions:
(t = fill(0.0), q = [1.0, 1.0], p = [0.5, 2.0], v = [0.0, 0.0])
Parameters:
(α = 5.0,)This problem can again be integrated using GeometricIntegrators:
sol = integrate(lprob, Gauss(1))
fig = lines(parent(sol.q[:,1]), parent(sol.q[:,2]);
axis = (; xlabel = "x₁", ylabel = "x₂", title = "Particle moving in a square potential"),
figure = (; size = (800,600), fontsize = 22))┌ Warning: Hermite Extrapolation: q's history[1] and history[2] are identical!
└ @ GeometricIntegratorsBase ~/.julia/packages/GeometricIntegratorsBase/bVkej/src/extrapolation/hermite.jl:206