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

[Sh_166]

So, as I mentioned in this topic, 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:

  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.

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. : )

[SteffenPL]

Could you do this?

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

That should work in the mixed case.

[Sh_166]

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?

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:

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!

[Sh_166]

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:

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:

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!

[Balinus]

Good material here for a blog post somewhere

[ziolai]

I agree to the blog post!

[Daniel_Lewis]

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

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.