Symbolics & MTK Pipeline - Lagrangians

Specific Domains Modelling & Simulations
,
1

Hello,

I am trying to use Symbolics.jl with ModelingToolkit. I’m confused on how to properly use the two packages together to define a Lagrangian that can be integrated in MTK. I suspect that my approach (broadly) is in the wrong direction.

For context, I haven’t taken a differential equations class. There’s a good chance that I am missing something conceptually here that would help me understand how to use the packages better.

As an exploration, I’ve tried constructing the Lagrangian for a particle with some kinetic energy while in a gravitational potential energy field. That code is shown below. I was able to solve the Euler-Lagrange equations using Symbolics.jl, but I’m stuck at actually using MTK to solve them.

The full markdown code is further down. Here are some screenshots from a Pluto notebook which show the relevant Latex:

image
image1555×1686 142 KB

using ModelingToolkit, OrdinaryDiffEq

### Setting up parameters, variables, and operators

@independent_variables t
@parameters m=1 G=1 M=1 # (I think that using m=1 it gives a value for the ODE solver later?)
@variables r(t)=10 θ(t)=1 ϕ(t)=1 # (I was doing r=1 etc. to try and specify an initial condition w/o u0)
D=Differential(t)

### Defining the Lagrangian

T = 1/2 * m * sum([D(r), D(θ)*r, r*sin(θ)*D(ϕ)].^2)

V = -m * G * M / r

L=T-V

### Euler-Lagrange Functions

euler_lagrange_eq(L, q) = expand_derivatives(D(Differential(D(q))(L)) - Differential(q)(L)) ~ 0

eq_r = euler_lagrange_eq(L, r)
eq_θ = euler_lagrange_eq(L, θ)
eq_ϕ = euler_lagrange_eq(L, ϕ)

eqs=[eq_r, eq_θ, eq_ϕ]

### Defining the System

@named sys=ODESystem(eqs, t)

sys2 = structural_simplify(sys)

unknowns(sys2) # θ got eliminated since it wasn't doing anything

### Setting up the problem

u0 = [1.0, 1.0, 100.0, 1.0, 1.0, 1.0, 0.0, 0.0] 
# Here I had to specify an initial condition for everything
# otherwise the solver wouldn't accept it. This seems wrong.

tspan=(0.0, 10.0)

prob = ODEProblem(sys2, u0, tspan)
# Error: "Initialization system is overdetermined. 2 equations for 0 unknowns." 
# That makes sense since I just specified 8 initial conditons, but if I don't do that
# then I get an error saying that the number of initial conditions isn't matching my unknowns.

Is there a different way I should approach this, or is there a subtlety I missed?

1 Like
2

The error is telling you you specified too many values and the solver can’t set up a system to solve for a consistent initial condition. Changing

u0 = []

Allows the system to initialize and solve, with a couple of warnings about initialization. These warnings are only a concern if you discover it picked a bad set of initial conditions.

What partial set of values did you try that resulted in a failure to initialize?

3

I think my main source of confusion was how to cleanly set up the system. Your comment about using

u0 = []

Helped me understand how complications earlier in the problem caused a ripple-effect that resulted in a very messy post-structural_simplify system.

Here’s the code I ended up with which works pretty well:

(edit: the vr, vtheta, and vphi variables defined at the start aren’t used. They’re fossils of an earlier test. I want to call them fossils because then I get to visualize an asteroid hitting the broken Pluto notebook I made earlier).

I think the most important aspect of this newer version is the line where I define dummy placeholder variables for the acceleration components so that I could use symbolic_linear_solve on the euler equations a line later.

image
image1920×1409 187 KB


using ModelingToolkit, OrdinaryDiffEq, Symbolics, LinearAlgebra, Plots, Nemo, CSV, DataFrames

using ModelingToolkit: t_nounits as t, D_nounits as D

###

@parameters g=9.8 m=1

@variables r(t) θ(t) ϕ(t)

Lagrangian = (1/2 * m * sum( [D(r), D(θ)*r, D(ϕ)*r*sin(θ)] .^ 2)) - (g * r * m)

###

build_euler_lagrange(L, q) = expand_derivatives( D(Differential( D(q) )(L) )+0.1 ~ Differential(q)(L) )

test_equation = build_euler_lagrange(Lagrangian, r)

euler_eqs=build_euler_lagrange.(Lagrangian, [r, θ, ϕ])

begin
	ar = D(D(r))
	aθ = D(D(θ))
	aϕ = D(D(ϕ))
end

eqs_of_motion = simplify.(symbolic_linear_solve.(euler_eqs, [ar, aθ, aϕ])) .~ [ar, aθ, aϕ]

@mtkbuild sys = ODESystem(eqs_of_motion, t)

unknowns(sys)

###

prob = ODEProblem(sys, [ 5.0, 1.0, 1.0, 0.1, 1.0, 3.0 ], (0.0, 100.0))

sol=ModelingToolkit.solve(prob, Tsit5())

df_sol = select(DataFrame(sol), 1:4)

###

df = DataFrame(
	timestamp=df_sol.timestamp,
	x= @.(pr * sin(pθ) * cos(pϕ)),
	y= @.(pr * sin(pθ) * sin(pϕ)),
	z= @.(pr * cos(pθ))
)

plot(df.x, df.y, df.z,
	seriestype=path3d,
	legend=false
	)

Thank you so much for your advice on this. Tomorrow I plan to try out a more complicated Lagrangian in spherical coordinates, hopefully it works just as well.