# Implementation of the Lagrangian formalism using Julia
**URL:** https://discourse.julialang.org/t/implementation-of-the-lagrangian-formalism-using-julia/89502
**Category:** General Usage
**Tags:** question, physics
**Created:** [October 29, 2022, 7:12pm UTC](https://discourse.julialang.org/t/implementation-of-the-lagrangian-formalism-using-julia/89502 "2022-10-29T19:12:48Z")
**Posts on this page:** 15
**Page:** 1
### Author:  [@Sh\_166](https://discourse.julialang.org/u/Sh_166)
#### Post date: [October 29, 2022, 7:12pm UTC](https://discourse.julialang.org/t/implementation-of-the-lagrangian-formalism-using-julia/89502/1 "2022-10-29T19:12:48Z")
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!!
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [October 29, 2022, 10:42pm UTC](https://discourse.julialang.org/t/implementation-of-the-lagrangian-formalism-using-julia/89502/2 "2022-10-29T22:42:20Z")
> [@Sh\_166](#):
>
> 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.
Yes. Replace SymPy with Symbolics.jl and `scipy.odeint` with DifferentialEquations.jl.
---
### Author:  [@Sh\_166](https://discourse.julialang.org/u/Sh_166)
#### Post date: [October 30, 2022, 2:47am UTC](https://discourse.julialang.org/t/implementation-of-the-lagrangian-formalism-using-julia/89502/3 "2022-10-30T02:47:37Z")
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.
---
### Author:  [@dlakelan](https://discourse.julialang.org/u/dlakelan)
#### Post date: [October 30, 2022, 3:17am UTC](https://discourse.julialang.org/t/implementation-of-the-lagrangian-formalism-using-julia/89502/4 "2022-10-30T03:17:00Z")
> [@Sh\_166](#):
>
> 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.
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?
---
### Author:  [@Sh\_166](https://discourse.julialang.org/u/Sh_166)
#### Post date: [October 30, 2022, 3:36am UTC](https://discourse.julialang.org/t/implementation-of-the-lagrangian-formalism-using-julia/89502/5 "2022-10-30T03:36:39Z")
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).
---
### Author:  [@SteffenPL](https://discourse.julialang.org/u/SteffenPL)
#### Post date: [October 30, 2022, 9:56am UTC](https://discourse.julialang.org/t/implementation-of-the-lagrangian-formalism-using-julia/89502/6 "2022-10-30T09:56:21Z")
- 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](https://diffeq.sciml.ai/stable/types/dynamical_types/#Hamiltonian-Problems)
- 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](https://github.com/SciML/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](https://diffeq.sciml.ai/stable/types/dynamical_types/#SciMLBase.SecondOrderODEProblem) where the terms M, B, \mathbf{F} are computed with Symbolics and converted into a fast Julia function via [Symbolics.build\_function](https://symbolics.juliasymbolics.org/stable/tutorials/symbolic_functions/#Building-Functions-1).
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](https://symbolics.juliasymbolics.org/stable/manual/expression_manipulation/#Symbolics.solve_for).
---
### Author:  [@zdenek\_hurak](https://discourse.julialang.org/u/zdenek_hurak)
#### Post date: [October 30, 2022, 11:47am UTC](https://discourse.julialang.org/t/implementation-of-the-lagrangian-formalism-using-julia/89502/7 "2022-10-30T11:47:52Z")
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](https://symbolics.juliasymbolics.org/dev/).
Alternatively, you may want to search for solvers that accept the problem is **mass matrix form** – [such solvers are available in Julia](https://diffeq.sciml.ai/stable/tutorials/dae_example/#Mass-Matrix-Differential-Algebraic-Equations-(DAEs)).
---
### Author:  [@Sh\_166](https://discourse.julialang.org/u/Sh_166)
#### Post date: [October 30, 2022, 12:26pm UTC](https://discourse.julialang.org/t/implementation-of-the-lagrangian-formalism-using-julia/89502/8 "2022-10-30T12:26:06Z")
Thanks a lot for your responses, guys!
> [@SteffenPL](#):
>
> 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?
@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](https://math24.net/double-pendulum.html)). 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!
---
### Author:  [@SteffenPL](https://discourse.julialang.org/u/SteffenPL)
#### Post date: [October 30, 2022, 1:13pm UTC](https://discourse.julialang.org/t/implementation-of-the-lagrangian-formalism-using-julia/89502/9 "2022-10-30T13:13:46Z")
I did a small example, **please check the formulas, they might be wrong**!
```julia
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.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [October 30, 2022, 1:39pm UTC](https://discourse.julialang.org/t/implementation-of-the-lagrangian-formalism-using-julia/89502/10 "2022-10-30T13:39:38Z")
> [@SteffenPL](#):
>
> I once also implemented a subset of the Lagrangian formalism for myself with [ModelingToolkit.jl](https://github.com/SciML/ModelingToolkit.jl). In that case, I also had a complex Lagrangian where mixed derivatives occur, like in your example.
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.
---
### Author:  [@SteffenPL](https://discourse.julialang.org/u/SteffenPL)
#### Post date: [October 30, 2022, 3:48pm UTC](https://discourse.julialang.org/t/implementation-of-the-lagrangian-formalism-using-julia/89502/11 "2022-10-30T15:48:28Z")
I didn’t want to say that there is an issue; I think everything works as it should 🙂
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](https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/connectors/connections/#Translational)
>
> ```julia
> 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())
>
> ```
---
### Author:  [@Sh\_166](https://discourse.julialang.org/u/Sh_166)
#### Post date: [October 30, 2022, 4:43pm UTC](https://discourse.julialang.org/t/implementation-of-the-lagrangian-formalism-using-julia/89502/12 "2022-10-30T16:43:54Z")
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!
---
### Author:  [@Balinus](https://discourse.julialang.org/u/Balinus)
#### Post date: [October 30, 2022, 7:11pm UTC](https://discourse.julialang.org/t/implementation-of-the-lagrangian-formalism-using-julia/89502/13 "2022-10-30T19:11:28Z")
Did you look at DynamicalSystems.jl? Not a specialist but perhaps is has some building blocks?
---
### Author:  [@Sh\_166](https://discourse.julialang.org/u/Sh_166)
#### Post date: [October 31, 2022, 12:13pm UTC](https://discourse.julialang.org/t/implementation-of-the-lagrangian-formalism-using-julia/89502/14 "2022-10-31T12:13:35Z")
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. 😅
---
### Author:  [@finmod](https://discourse.julialang.org/u/finmod)
#### Post date: [October 31, 2023, 11:57am UTC](https://discourse.julialang.org/t/implementation-of-the-lagrangian-formalism-using-julia/89502/15 "2023-10-31T11:57:26Z")
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ᶠ = p_y - w_l - (Ŝ - S)^2
Uʰ = c^(γ) - (L̂ - l)^2 - (M̂ʰ - Mʰ)^2
Zₕ = w_l - p_c - mʰ
Zₚ = p_c - w_l - 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))