Implementation of the Lagrangian formalism using Julia

General Usage
,
1

Hi, I am a newcomer to Julia and I am specifically interested in applying Julia to physics problems. I have already worked through differential equations and plotting and stuff, and I am quite drawn by the language and how fast it is.

I was wondering if there was a package or a method to implement the Lagrangian formalism using Julia.
For example, in Python, one can use SciPy and SymPy to define variables and the Lagrangian (as an equation, i.e., L = … ). We can then define the equations of motion as:

E_1 = sp.diff(L,x1) - sp.diff(sp.diff(L,x1_d),t).simplify()

And so on, for the other coordinates. These systems of equations can then be solved using sp.solve(), odeint() and so on.

From what I have seen, to solve DEs in Julia, a major requirement is that one should be able to write out the differential equation in a separable form, i.e.,
D(x) = \ ... , D being the derivative.
However, in most complicated scenarios equations might be coupled and it may not be possible to write them out in this manner. Is there a way that one can define a Lagrangian for a system and then have Julia compute the derivatives, equations of motions and then solve them?
Thanks in advance!!

2 Likes
2

Yes. Replace SymPy with Symbolics.jl and scipy.odeint with DifferentialEquations.jl.

3 Likes
3

Doesn’t that still have the requirement of being able to cleanly separate the derivatives and bring them on the left?
A second problem will be with problems that are coupled, dx and dy appear in the same equation and can’t be cleanly separated into 2 different equations, one for x and the other for y.

I’ve tried Symbolics.jl though; the github readme clearly stated that it shouldn’t be used for anything important though, since it’s a work in progress.

4

In the Lagrangian formalism you differentiate by one variable at a time. You should never have derivatives of multiple variables in the same equation no?

2 Likes
5

Actually, it is possible to have derivatives of different variables in the same equation.

For example, if you take a system in spherical polar coordinates, we can have terms like r^{2}\sin^{2}\theta (\frac{d\phi}{dt})^{2} in the KE. For the \phi coordinate, differentiating with respect to d\phi/dt will lead to both r and \sin\theta in the same equation. The Euler Lagrange equations are of the form,
\frac{d}{dt}(\frac{\partial L}{\partial \omega_{\phi}}) - \frac{\partial L}{\partial \phi} = 0
(\omega_{\phi} = d\phi/dt)
So, now we can time derivatives of r and \sin \theta, so, multiple derivatives in a single equation.
It is because of the nature of the metric in 3D Cartesian coordinates that we usually do not have this complication. But even then, for complicated enough potentials, this may not hold (for example, something along the lines of x dy/dt, a velocity dependent coupled potential).

6

  • First of all, if that is an option for you, you could directly use the Hamiltonian formalism, which is part of DifferentialEquations.jl: Dynamical, Hamiltonian and 2nd Order ODE Problems · DifferentialEquations.jl

  • How do you exactly couple the systems? Do you have algebraic constraints like g(q_1, q_2) = 0 and then the Lagrangian becomes L = L_1 + L_2 + \lambda g where \lambda is the Lagrangian multiplier?

  • I once also implemented a subset of the Lagrangian formalism for myself with ModelingToolkit.jl. In that case, I also had a complex Lagrangian where mixed derivatives occur, like in your example.

What I did was to help the computations a bit in the following way:
Given variables q_1, \dots, q_n, let me write p_i = \frac{\partial L}{\partial \dot q_i}(t, \mathbf{q}, \dot{ \mathbf{q}} ) then we get

\frac{\mathrm{d} }{\mathrm{d} t} p_i = \frac{\partial p_i}{\partial t} + \sum_{j=1}^n \frac{\partial p_i}{\partial q_j} \frac{\mathrm{d} q_j}{\mathrm{d} t} + \sum_{j=1}^n \frac{\partial p_i}{\partial \dot q_j} \frac{\mathrm{d} \dot q_j}{\mathrm{d} t}.

This term is now linear in the velocity and acceleration variables. So, you could write maybe B = (\frac{\partial p_i}{\partial q_j})_{i,j} \in \mathbb{R}^{n \times n}, \mathbf{F} = (\frac{\partial L}{\partial q_i} - \frac{\partial p_i}{\partial t} )_{i} \in \mathbb{R}^{n} and M = (\frac{\partial p_i}{\partial \dot q_j})_{i,j} \in \mathbb{R}^{n \times n} to obtain the system

M \frac{\mathrm{d^2} \mathbf{q}}{\mathrm{d} t^2} = \mathbf F - B \frac{\mathrm{d} \mathbf{q}}{\mathrm{d} t}.

In the form you can feed it directly into a SecondOrderODEProblem Dynamical, Hamiltonian and 2nd Order ODE Problems · DifferentialEquations.jl where the terms M, B, \mathbf{F} are computed with Symbolics and converted into a fast Julia function via Symbolics.build_function.

Probably one could also get the same result by just solving a few equations in a smart way, e.g. using Expression Manipulation · Symbolics.jl.

3 Likes
7

So, what you are aiming at is a (vector) ODE with explicitly expressed first derivatives of the variables (aka state equations), that is,

\begin{bmatrix} \dot x_1 \\ \vdots \\ \dot x_n \end{bmatrix} = \begin{bmatrix} f_1(\boldsymbol x) \\ \vdots \\ f_n(\boldsymbol x) \end{bmatrix},

while what the Lagrange equation gives you is this

\underbrace{ \begin{bmatrix} m_{11}(\boldsymbol q) & \ldots & m_{1n}(\boldsymbol q)\\ \vdots \\ m_{n1}(\boldsymbol q) & \ldots & m_{nn}(\boldsymbol q) \end{bmatrix}}_{\boldsymbol M(\boldsymbol q)} \begin{bmatrix} \ddot q_1 \\ \vdots \\ \ddot q_n \end{bmatrix} = \ldots

First, reducing the order of these ODEs from 2 to 1 is done by doubling the number of variables through

\boldsymbol x = \begin{bmatrix} q \\ \dot q\end{bmatrix}.

The first-order ODEs are then

\begin{aligned} \dot{\boldsymbol{x}}_1 &= \boldsymbol x_2\\ \boldsymbol M( \boldsymbol x_1)\dot{\boldsymbol{x}}_2 &= \ldots \end{aligned}

Now, if you insist on obtaining the explicit ODEs (aka state equations), you have to solve the second (vector) equation for \dot{\boldsymbol{x}}_2. That is, you have to find the inverse of \boldsymbol M. Some package for symbolic computation might help here, but I am not fluent with them (yet). The crucial thing is that in order to solve the equation(s) symbolically, you have to treat the individual entries of the vector \dot{\boldsymbol{x}}_2 as simple variables (you are really solving for \dot{\boldsymbol{x}}_2 now and not for \boldsymbol{x}_2). I guess this might require some renaming for the purpose of solving the equations, but once again, I am not familiar with Symbolics.jl.

Alternatively, you may want to search for solvers that accept the problem is mass matrix form – such solvers are available in Julia.

2 Likes
8

Thanks a lot for your responses, guys!

@SteffenPL, well, I did not really have any particular problem in mind and was just asking in a very general sense. But what I had in mind with regards to the coupling was more like an explicitly non linear system, like the double pendulum where we have coupling terms in the velocities (ref: here). However, I will keep your suggestions in mind about Symbolics.jl. Thanks a lot!

@zdenek_hurak, yes, that is exactly what I would like to do. I implemented the introducing new variables method earlier, with simpler examples and they worked, but the caveat was that my examples were simple enough (at least in the lower orders) to be written down and separated by hand. For more problematic cases, I would prefer it if the code took care of the entire thing after defining the Lagrangian, so I would check out what you suggested, though I am not familiar with matrix manipulation in Julia yet. I have also tried a bit of the ModelingToolkit.jl, and there might be something in there as well. Thanks a lot for your help!

9

I did a small example, please check the formulas, they might be wrong!

using Symbolics: gradient, jacobian, derivative, scalarize
using OrdinaryDiffEq
using LinearAlgebra

L(t,x,v) = dot(v,v) - dot(x,x)

n = 4
@variables t, x[1:n], v[1:n]

p = gradient(L(t,x,v), v)
M = jacobian(p, v)
B = jacobian(p, x)
F = gradient(L(t,x,v),x) - derivative.(p, t)
rhs = scalarize( F - B*v )

M_fnc = eval( build_function(M, [x,v,t])[1] )
rhs_fnc! = eval( build_function(rhs, [x,v,t])[2] ) # in-place variant

function f_ode(ddx, dx, x, p, t)
    rhs_fnc!(ddx, [x,dx,t] )
    ddx .= M_fnc( [x,dx,t] ) \ ddx
end

tspan = (0.0, 10.0)
x0 = [1.0,1,1,1]
v0 = [0.0,0,0,0]
prob = SecondOrderODEProblem(f_ode, v0, x0, tspan)
sol = solve(prob, SymplecticEuler(), dt = 0.1)

One part which is not that great is that one has to solve each step a linear system.
I couldn’t find a way to pass non-constant mass matrices to an ODEProblem; in cases where the mass is constant, one could, of course, pass the mass matrix so some ODE/DAE solver.

2 Likes
10

Is there a current open issue? Yingbo has been working a lot on the state selection algorithms and IIRC partial state selection was the last bit needed for structural_simplify to generate stable forms for this kind of thing.

1 Like
11

I didn’t want to say that there is an issue; I think everything works as it should :slight_smile:

Of course, for plugging mechanical components together, connect + structural_simplify works well, it’s just a bit harder to learn.

e.g. two springs with connect and MKL components

Following this About Acausal Connections · ModelingToolkitStandardLibrary.jl

using ModelingToolkit
using ModelingToolkitStandardLibrary
using OrdinaryDiffEq

TV = ModelingToolkitStandardLibrary.Mechanical.TranslationalPosition

@parameters t

@named center = TV.Fixed(s_0 = 0)
@named spring_1 = TV.Spring(Val(:relative), k=1.0, delta_s_0 = 0.1)
@named mass_1 = TV.Mass(m=1, v_0 = 0, s_0 = 1.0)
@named spring_2 = TV.Spring(Val(:relative), k=4.0, delta_s_0 = -0.1)
@named mass_2 = TV.Mass(m=1, v_0 = 0, s_0 = 2.0)

eqs = [
    connect(center.flange, spring_1.flange_a),
    connect(mass_1.flange, spring_1.flange_b),
    connect(mass_1.flange, spring_2.flange_a),
    connect(mass_2.flange, spring_2.flange_b),
]

@named model = ODESystem(eqs, t; systems=[center, spring_1, mass_1, spring_2, mass_2])

sys = structural_simplify(model)
prob = ODEProblem(sys, [], (0, 10.0), [] )
sol = solve(prob, Tsit5())
1 Like
12

Thanks a ton, @SteffenPL! I will surely check this out once I have some time and let you know if I have any questions. Once more, thanks!

1 Like
13

Did you look at DynamicalSystems.jl? Not a specialist but perhaps is has some building blocks?

1 Like
14

I did go through a bit of the DynamicalSystems.jl documentation, but didn’t find anything that seemed applicable to the problem. Of course, something might have escaped my attention, but as far as I recall, there was nothing that seemed suitable. But, frankly, I have a lot of reading to do before I can say that it isn’t applicable. :sweat_smile:

2 Likes
15

I have tried to use your example above for a Lagrangian D’Alembert dynamical system with nonholonomic constraints below. As you can see, it does the job but it is ugly! Can you help me to shape up this fairly standard Lagrangian problem into a simpler re-usable function?
…Julia
using DifferentialEquations
using OrdinaryDiffEq
using ModelingToolkit
using Plots; gr()

ps = @parameters α β dp̂ γ K₀ l₀ L̂ Mᶠ₀ Mʰ₀ M̂ʰ p₀ S₀ Ŝ w₀ μ₁ μ₂ μ₃ μ₄ μ₅ μ₆ μ₇ μ₈ μ₉ μ₁₀ 𝛌₁ 𝛌₂ 𝛌₃

@variables t
states = @variables c(t) K(t) l(t) Mᶠ(t) Mʰ(t) p(t) S(t) w(t) k(t) mᶠ(t) mʰ(t) s(t)
@variables dep(t) y(t) Uᶠ(t) Uʰ(t) Zₕ(t) Zₚ(t) Zᵣ(t) 𝛌₁(t) 𝛌₂(t) 𝛌₃(t)

Dₜ = Differential(t)

dep = dp̂K
y = βK^(1 - α) * l^(α)

Uᶠ = py - wl - (Ŝ - S)^2
Uʰ = c^(γ) - (L̂ - l)^2 - (M̂ʰ - Mʰ)^2

Zₕ = wl - pc - mʰ
Zₚ = pc - wl - mᶠ
Zᵣ = y - c - k - s - dep

DLP = [Uᶠ, Uʰ, 𝛌₁*Zₕ, 𝛌₂*Zₚ, 𝛌₃*Zᵣ]

∂DLPᵀ = transpose(Symbolics.jacobian(DLP, states))

eqs = [Dₜ(c) ~ μ₁∂DLPᵀ[1,2] + ∂DLPᵀ[1,3] + ∂DLPᵀ[1,4] + ∂DLPᵀ[1,5],
Dₜ(K) ~ μ₂∂DLPᵀ[2,1] + ∂DLPᵀ[2,3] + ∂DLPᵀ[2,4] + ∂DLPᵀ[9,5],
Dₜ(l) ~ μ₃∂DLPᵀ[3,1] + μ₄∂DLPᵀ[3,2] + ∂DLPᵀ[3,3] + ∂DLPᵀ[3,4] + ∂DLPᵀ[3,5],
Dₜ(Mᶠ) ~ μ₅∂DLPᵀ[4,1] + ∂DLPᵀ[4,3] + ∂DLPᵀ[10,4]+ ∂DLPᵀ[4,5],
Dₜ(Mʰ) ~ μ₆∂DLPᵀ[5,2] + ∂DLPᵀ[11,3]+ ∂DLPᵀ[5,4] + ∂DLPᵀ[5,5],
Dₜ(p) ~ μ₇∂DLPᵀ[6,1] + ∂DLPᵀ[6,3] + ∂DLPᵀ[6,4] + ∂DLPᵀ[6,5],
Dₜ(S) ~ μ₈∂DLPᵀ[7,1] + ∂DLPᵀ[7,3] + ∂DLPᵀ[7,4] + ∂DLPᵀ[12,5],
Dₜ(w) ~ μ₉*∂DLPᵀ[8,1] + ∂DLPᵀ[8,3] + ∂DLPᵀ[8,4] + ∂DLPᵀ[8,5],
Dₜ(K) ~ k,
Dₜ(Mᶠ) ~ mᶠ,
Dₜ(Mʰ) ~ mʰ,
Dₜ(S) ~ s,
0. ~ Zₕ,
0. ~ Zₚ,
0. ~ Zᵣ]

defs = Dict(Uᶠ => 1.0, Uʰ => 1.0, dep => 1.0, y => 1.0, c => β*K₀^(1 - α) * l₀^(α), K => K₀, l => l₀, Mᶠ => Mᶠ₀, Mʰ => Mʰ₀, p => p₀, S => S₀, w => w₀, k => 0.1, mᶠ => 0, mʰ => 0, s => 0, 𝛌₁ => 0.01, 𝛌₂ => 0.01, 𝛌₃ => 0.02,
α => 0.5, β => 1.0, dp̂ => 0.1, γ => 0.5, K₀ => 1.6, l₀ => 1.0, L̂ => 1.0, Mᶠ₀ => 1.05, Mʰ₀ => 1.2, M̂ʰ => 1.0, p₀ => 0.8, S₀ => 0.4, Ŝ => 0.1, w₀ => 1.2, μ₁ => 1.0, μ₂ => 1.0, μ₃ => 1.0, μ₄ => 1.0, μ₅ => 1.0, μ₆ => 1.0, μ₇ => 1.0, μ₈ => 1.0, μ₉ => 0.0, μ₁₀ => 0.0)

@named A1_dae = ModelingToolkit.ODESystem(eqs, t, states, ps; defaults = defs, checks = false)
…Julia

A1_ode = structural_simplify(dae_index_lowering(A1_dae))