# Symbolics.jl, ModelingToolkit.jl, DifferentialEquations.jl and the Lagrangian approach

**URL:** <https://discourse.julialang.org/t/symbolics-jl-modelingtoolkit-jl-differentialequations-jl-and-the-lagrangian-approach/89629>\
**Category:** General Usage\
**Tags:** physics, modelingtoolkit, symbolics, differentialequation\
**Created:** [November 1, 2022, 7:23pm UTC](https://discourse.julialang.org/t/symbolics-jl-modelingtoolkit-jl-differentialequations-jl-and-the-lagrangian-approach/89629 "2022-11-01T19:23:41Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![Sh\_166](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sh_166/32/43965_2.png) [@Sh\_166](https://discourse.julialang.org/u/Sh_166)\
**Post date:** [November 1, 2022, 7:23pm UTC](https://discourse.julialang.org/t/symbolics-jl-modelingtoolkit-jl-differentialequations-jl-and-the-lagrangian-approach/89629/1 "2022-11-01T19:23:42Z")

</div>

So, as I mentioned in this [topic](https://discourse.julialang.org/t/implementation-of-the-lagrangian-formalism-using-julia/89502/14), I was trying to implement the Lagrangian approach from scratch in Julia, and inquired for existing ways to do that. As the comments on that discussion point out, nothing of the sort exists yet, but the different packages suggested were `Symbolics.jl`, `ModelingToolkit.jl`, `DifferentialEquations.jl`.

I thought about it a bit and came up with a way to use these packages to implement simpler Lagrangians without much hassle.  
The general idea is to use the `Symbolics.jl` to define 3 kinds of derivatives, d/dt, d/dx, d/d\dot{x}. The Lagrangian can be defined in terms of T and V properly, by using `expand_derivatives()` in conjunction with the above. Then, we can define the Euler - Lagrange equations, use `Symbolics.solve_for()` to get the double derivatives with which we can define our ODESystem.

An example, for the isotropic Harmonic Oscillator in 2D, is given below:

```julia
  using Symbolics, Plots, DifferentialEquations, ModelingToolkit

  @variables t x(t) y(t) 
  @parameters m, k

  d = Differential(t)
  D1 = Differential(x)
  D2 = Differential(y)
  D11 = Differential(d(x))
  D22 = Differential(d(y))

T = 0.5*m*(expand_derivatives(d(x))^2 + expand_derivatives(d(y))^2) #KE
V = 0.5*k*(x^2 + y^2) #Potential energy

L = T - V

E₁ = expand_derivatives(d(D11(L))) - expand_derivatives(D1(L))
E₁₁ = Symbolics.solve_for(E₁ ~ 0, d(d(x)))

E₂ = expand_derivatives(d(D22(L))) - expand_derivatives(D2(L))
E₂₂ = Symbolics.solve_for(E₂ ~ 0, d(d(y)))

system = [d(d(x)) ~ E₁₁,
          d(d(y)) ~ E₂₂]

@named HO = ODESystem(system, t, [x, y], [m, k])

HO = structural_simplify(HO)

X_0 = [d(x) => 0,
       d(y) => 0,
       x => 10,
       y => 10]

M = [ m => 1
      k => 1 ]

tspan = (0.0, 100.0)

prob = ODEProblem(HO, X_0, tspan, M, abstol = 1e-8, reltol = 1e-8)
sol = solve(prob)

```

And, done.  
I have plotted this solution below against an expected one.

 ![Trial](https://global.discourse-cdn.com/julialang/original/3X/f/2/f2a7c9ddd7cdfe7a7c624b0488b28874a57a2b1f.png)

However, this actually doesn’t solve my original problem, because I have coupled derivatives in there (so the equation for d^{2}x/dt^{2} will contain double derivatives of y and z as well, and same for the other 2 coordinates as well). I would appreciate it if someone can think of a way around that issue.

And, if not, this post will hopefully save someone some time in the future. : )

---

<div class="post-metadata">

**Author:** ![SteffenPL](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/steffenpl/32/206270_2.png) [@SteffenPL](https://discourse.julialang.org/u/SteffenPL)\
**Post date:** [November 1, 2022, 10:08pm UTC](https://discourse.julialang.org/t/symbolics-jl-modelingtoolkit-jl-differentialequations-jl-and-the-lagrangian-approach/89629/2 "2022-11-01T22:08:20Z")

</div>

Could you do this?

```julia
E₁₁, E₂₂ = Symbolics.solve_for( [E₁ ~ 0, E₂ ~ 0], [d(d(x)), d(d(y))] )

```

That should work in the mixed case.

---

<div class="post-metadata">

**Author:** ![Sh\_166](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sh_166/32/43965_2.png) [@Sh\_166](https://discourse.julialang.org/u/Sh_166)\
**Post date:** [November 2, 2022, 3:13am UTC](https://discourse.julialang.org/t/symbolics-jl-modelingtoolkit-jl-differentialequations-jl-and-the-lagrangian-approach/89629/3 "2022-11-02T03:13:19Z")

</div>

Thanks for your suggestion, @SteffenPL. I tried it out with my system, but unfortunately I get the error `LinearAlgebra.SingularException(2)`. I am not much familiar with Linear Algebra in Julia, but even then this seems a bit weird considering the fact that I haven’t assigned any values to my parameters yet at this point. I actually have a system with 3 degrees of freedom, and the Euler Lagrange equations are very messy, coupled and nonlinear. Any ideas on what might be causing this?

```julia
using Symbolics, Plots, DifferentialEquations, ModelingToolkit
@variables t θ(t) ϕ(t) r(t) 
@parameters m₁, m₂, l₁, l₂, k, g
d = Differential(t)
D1 = Differential(θ)
D2 = Differential(ϕ)
D3 = Differential(r)
D11 = Differential(d(θ))
D22 = Differential(d(ϕ))
D33 = Differential(d(r))

x₁ = l₁*sin(θ)
x₂ = (l₂ + r)*sin(ϕ) + l₁*sin(θ)
y₁ = -l₁*cos(θ)
y₂ = -(l₂ + r)*cos(ϕ) - l₁*cos(θ)

T = 0.5*m₁*(expand_derivatives(d(x₁))^2 + expand_derivatives(d(y₁))^2) + 0.5*m₂*(expand_derivatives(d(x₂))^2 + expand_derivatives(d(y₂))^2)
T = simplify(T)

V = m₁*g*y₁ + m₂*g*y₂ + 0.5*k*r^2

```

The rest of the code is the same as I used above, albeit with an extra equation to account for the third degree of freedom. And, then, I used this:

```julia
E₁₁, E₂₂, E₃₃ = Symbolics.solve_for([E₁ ~ 0, E₂ ~ 0, E₃ ~ 0], [d(d(θ)), d(d(ϕ)), d(d(r))])

```

At this point, I get the error.  
I had worked out the case analytically, quite some time back, and I am pretty certain Julia returns the proper equations of motion after differentiating the Lagrangian.

Edit: On closer inspection, I noticed something weird. One is that the equations of motion contain terms like dl\_1/dt (should be 0 as lengths don’t change) and the other is that none of the 3 equations contain a d^{2}\phi/dt^{2} term which is extremely weird (and prolly causes the error in the first place). The first one may be because I defined lengths as parameters (though I am not sure how to work around it; plus this didn’t happen with the harmonic oscillator above), but I have no idea about how the second one is cropping up. Any help would be appeciated!

---

<div class="post-metadata">

**Author:** ![Sh\_166](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sh_166/32/43965_2.png) [@Sh\_166](https://discourse.julialang.org/u/Sh_166)\
**Post date:** [November 2, 2022, 1:40pm UTC](https://discourse.julialang.org/t/symbolics-jl-modelingtoolkit-jl-differentialequations-jl-and-the-lagrangian-approach/89629/4 "2022-11-02T13:40:57Z")

</div>

I managed to figure out the problem.  
Apparently, this line was causing issues:  
`E₁ = expand_derivatives(d(D11(L))) - expand_derivatives(D1(L))`  
For the case of the Harmonic oscillator, this works fine, but for a more complicated system, for some reason it runs into issues (it completely ignored d\phi/dt and tried to differentiate l\_1). I used the alternative:

```julia
E₁ = expand_derivatives(d(expand_derivatives(D11(L)))) - expand_derivatives(D1(L)) 

```

This worked fine. @SteffenPL’s suggestion worked as well. Unfortunately though, the kernel just hung up after a point and I had to reload the page. :")  
I wonder how Python managed to deal with this without hanging; is there a way to make the code run better? Since I am new to this, I expect I am writing a lot of memory heavy stuff which might be the reason for this.  
To be specific, it is after I run this piece of code that everything hangs:

```julia
system = [d(d(θ)) ~ E₁₁,
          d(d(ϕ)) ~ E₂₂,
          d(d(r)) ~ E₃₃]

```

**EDIT:** Apparently, nothing had hanged, the math rendering of the output had slowed everything down and that was the reason everything was unresponsive. After that though, everything worked quite fast and smooth and I managed to do what I wanted to. Even then, the general question of efficiency and speed remains. Thanks to everyone for their help!

---

<div class="post-metadata">

**Author:** ![Balinus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/balinus/32/243_2.png) [@Balinus](https://discourse.julialang.org/u/Balinus)\
**Post date:** [November 2, 2022, 4:39pm UTC](https://discourse.julialang.org/t/symbolics-jl-modelingtoolkit-jl-differentialequations-jl-and-the-lagrangian-approach/89629/5 "2022-11-02T16:39:47Z")

</div>

Good material here for a blog post somewhere 🙂

---

<div class="post-metadata">

**Author:** ![ziolai](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ziolai/32/23422_2.png) [@ziolai](https://discourse.julialang.org/u/ziolai)\
**Post date:** [January 4, 2023, 3:30pm UTC](https://discourse.julialang.org/t/symbolics-jl-modelingtoolkit-jl-differentialequations-jl-and-the-lagrangian-approach/89629/6 "2023-01-04T15:30:03Z")

</div>

I agree to the blog post!

---

<div class="post-metadata">

**Author:** ![Daniel\_Lewis](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/daniel_lewis/32/31250_2.png) [@Daniel\_Lewis](https://discourse.julialang.org/u/Daniel_Lewis)\
**Post date:** [January 12, 2023, 10:36am UTC](https://discourse.julialang.org/t/symbolics-jl-modelingtoolkit-jl-differentialequations-jl-and-the-lagrangian-approach/89629/7 "2023-01-12T10:36:57Z")

</div>

hi,  
i’m trying to simulate a cart pole using the code example taken from this topic.  
my code looks like this:

```julia
using Symbolics, Plots, DifferentialEquations, ModelingToolkit

n = 2 # number of general coordinates

vars = @variables t (q(t))[1:n] # q[1] = x, q[2] = θ
params = @parameters m₁ m₂ l g F

Q = [F*cos(t), 0]# general forces

t = vars[1]
qs = vars[2]

dt = Differential(t)
dqs = [Differential(q) for q in qs]
dq̇s = [Differential(dt(q)) for q in qs]

q̇s = [dt(q) for q in qs]

v = q̇s[2]*l

T = 0.5*(m₁+m₂)*q̇s[1]^2 + 0.5*m₂*v^2 + m₂*l*q̇s[1]*q̇s[2]*cos(qs[2])
U = -m₂*g*l*cos(qs[2])

L = T - U

Eᵢ = Vector{Num}()
Eᵢᵢ = Vector{Num}()

for (ind, (q, dq̇, dq)) in enumerate(zip(qs, dq̇s, dqs))
    push!(Eᵢ, expand_derivatives(dt(expand_derivatives(dq̇(L)))) - expand_derivatives(dq(L)))
    push!(Eᵢᵢ, Symbolics.solve_for(Eᵢ[ind] ~ Q[ind], dt(dt(q))))
end

system = [dt(dt(q)) ~ Eᵢᵢ[ind] for (ind, q) in enumerate(qs)]

@named HO = ODESystem(system, t, qs, params)

HO = structural_simplify(HO) # HERE IS THE PROBLEM : Model HO with 5 equations

X_0 = [qs[1] => 0, 
        qs[2] => 0,
        q̇s[1] => 0,
        q̇s[2] => 0]

M = [ m₁ => 30, # kg
        m₂ => 10, # kg
        l => 30, # m
        g => 9.81, # m/s^2
        F => 0] # N

tspan = (0.0, 100.0)

prob = ODEProblem(HO, X_0, tspan, M, abstol = 1e-8, reltol = 1e-8)
sol = solve(prob)

```

my problem is this: when I call the ` structural_simplify(HO)` I get 5 equations which souled be 4.  
what am I missing and what did I do wrong?

thanks in ahead.