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.
Before any use, we need to load EulerLagrange:
using EulerLagrange
Next, 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
\[ \begin{equation}
+\frac{1}{2} \left( v(t)_1^{2} + v(t)_2^{2} \right) - \frac{1}{2} \left( x(t)_1^{2} + x(t)_2^{2} \right)
+\end{equation}
+ \]
This Lagrangian together with the symbolic variables is then used to construct a LagrangianSystem:
lag_sys = LagrangianSystem(L, t, x, v)
+Lagrangian system with
+
+L = (1//2)*((v(t))[1]^2 + (v(t))[2]^2) - (1//2)*((x(t))[1]^2 + (x(t))[2]^2)
+
+and equations of motion
+
+-(x(t))[1] - Differential(t)(Differential(t)((x(t))[1]))
+-(x(t))[2] - Differential(t)(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₀ = [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", (0xfca7c8ae, 0xfa7e27aa, 0x95adf5b8, 0xff38bc48, 0xee2992dd), Expr}(quote
+ #= none:1 =#
+ #= none:5 =#
+ begin
+ #= none:7 =#
+ #= none:7 =# @inbounds begin
+ #= none:9 =#
+ ˍ₋out[1] = getindex(V, 1)
+ #= none:10 =#
+ ˍ₋out[2] = getindex(V, 2)
+ #= none:12 =#
+ nothing
+ end
+ end
+end)
+ f = RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :t, :X, :V, :params), EulerLagrange.var"#_RGF_ModTag", EulerLagrange.var"#_RGF_ModTag", (0xaef78b2c, 0x6cf956c4, 0xc224cee9, 0x7427de99, 0xba543221), Expr}(quote
+ #= none:1 =#
+ #= none:5 =#
+ begin
+ #= none:7 =#
+ #= none:7 =# @inbounds begin
+ #= none:9 =#
+ ˍ₋out[1] = -1 // 1 * getindex(X, 1)
+ #= none:10 =#
+ ˍ₋out[2] = -1 // 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", (0x5703091c, 0x9ce67ff9, 0x290d6f43, 0xe14957ba, 0xf0a2a78e), Expr}(quote
+ #= none:1 =#
+ #= none:5 =#
+ begin
+ #= none:7 =#
+ #= none:7 =# @inbounds begin
+ #= none:9 =#
+ ˍ₋out[1] = (Differential(t))(getindex(V, 1))
+ #= none:10 =#
+ ˍ₋out[2] = (Differential(t))(getindex(V, 2))
+ #= none:12 =#
+ nothing
+ end
+ end
+end)
+
+ Lagrangian: L = RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:t, :X, :V, :params), EulerLagrange.var"#_RGF_ModTag", EulerLagrange.var"#_RGF_ModTag", (0x91bc2205, 0xb9d3a437, 0xaefd9ff7, 0xafa91fbc, 0x18ba0e5e), Expr}(quote
+ #= none:1 =#
+ #= none:5 =#
+ 1 // 2 * (getindex(V, 1) ^ 2 + getindex(V, 2) ^ 2) + -1 // 2 * (getindex(X, 1) ^ 2 + getindex(X, 2) ^ 2)
+end)
+
+ Invariants:
+ GeometricBase.NullInvariants()
+
+ Timespan: (0.0, 10.0)
+ Timestep: 0.01
+
+ Initial conditions:
+ (t = 0.0, q = [1.0, 1.0], p = [0.5, 2.0], λ = [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: q₀ and q₁ in initial guess are identical!
+└ @ GeometricIntegrators.Integrators ~/.julia/packages/GeometricIntegrators/Ptffa/src/initial_guess/hermite.jl:45
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) / 2
\[ \begin{equation}
+\frac{1}{2} \left( v(t)_1^{2} + v(t)_2^{2} \right) - \frac{1}{2} \left( x(t)_1^{2} + x(t)_2^{2} \right) \alpha_p
+\end{equation}
+ \]
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 = (1//2)*((v(t))[1]^2 + (v(t))[2]^2) - (1//2)*((x(t))[1]^2 + (x(t))[2]^2)*αₚ
+
+and equations of motion
+
+-Differential(t)(Differential(t)((x(t))[1])) - (x(t))[1]*αₚ
+-Differential(t)(Differential(t)((x(t))[2])) - (x(t))[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", (0xfca7c8ae, 0xfa7e27aa, 0x95adf5b8, 0xff38bc48, 0xee2992dd), Expr}(quote
+ #= none:1 =#
+ #= none:5 =#
+ begin
+ #= none:7 =#
+ #= none:7 =# @inbounds begin
+ #= none:9 =#
+ ˍ₋out[1] = getindex(V, 1)
+ #= none:10 =#
+ ˍ₋out[2] = getindex(V, 2)
+ #= none:12 =#
+ nothing
+ end
+ end
+end)
+ f = RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :t, :X, :V, :params), EulerLagrange.var"#_RGF_ModTag", EulerLagrange.var"#_RGF_ModTag", (0x1a87a248, 0x2bc6b04e, 0x2dc18435, 0x778f853e, 0x62d0c494), Expr}(quote
+ #= none:1 =#
+ #= none:5 =#
+ begin
+ #= none:7 =#
+ #= none:7 =# @inbounds begin
+ #= none:9 =#
+ ˍ₋out[1] = (-1 // 1 * getindex(X, 1)) * params.α
+ #= none:10 =#
+ ˍ₋out[2] = (-1 // 1 * getindex(X, 2)) * params.α
+ #= none:12 =#
+ nothing
+ end
+ end
+end)
+ g = RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :t, :X, :Λ, :V, :params), EulerLagrange.var"#_RGF_ModTag", EulerLagrange.var"#_RGF_ModTag", (0x5703091c, 0x9ce67ff9, 0x290d6f43, 0xe14957ba, 0xf0a2a78e), Expr}(quote
+ #= none:1 =#
+ #= none:5 =#
+ begin
+ #= none:7 =#
+ #= none:7 =# @inbounds begin
+ #= none:9 =#
+ ˍ₋out[1] = (Differential(t))(getindex(V, 1))
+ #= none:10 =#
+ ˍ₋out[2] = (Differential(t))(getindex(V, 2))
+ #= none:12 =#
+ nothing
+ end
+ end
+end)
+
+ Lagrangian: L = RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:t, :X, :V, :params), EulerLagrange.var"#_RGF_ModTag", EulerLagrange.var"#_RGF_ModTag", (0xe0f4c210, 0x22495478, 0xbbedb468, 0xdf30cd50, 0x25997ada), Expr}(quote
+ #= none:1 =#
+ #= none:5 =#
+ 1 // 2 * (getindex(V, 1) ^ 2 + getindex(V, 2) ^ 2) + (-1 // 2 * (getindex(X, 1) ^ 2 + getindex(X, 2) ^ 2)) * params.α
+end)
+
+ Invariants:
+ GeometricBase.NullInvariants()
+
+ Timespan: (0.0, 10.0)
+ Timestep: 0.01
+
+ Initial conditions:
+ (t = 0.0, q = [1.0, 1.0], p = [0.5, 2.0], λ = [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: q₀ and q₁ in initial guess are identical!
+└ @ GeometricIntegrators.Integrators ~/.julia/packages/GeometricIntegrators/Ptffa/src/initial_guess/hermite.jl:45
