# Symbolics & MTK Pipeline - Lagrangians

**URL:** <https://discourse.julialang.org/t/symbolics-mtk-pipeline-lagrangians/129730>\
**Category:** Modelling & Simulations\
**Tags:** question, package\
**Created:** [June 7, 2025, 11:40pm UTC](https://discourse.julialang.org/t/symbolics-mtk-pipeline-lagrangians/129730 "2025-06-07T23:40:26Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![Gray\_GALLANT](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gray_gallant/32/214253_2.png) [@Gray\_GALLANT](https://discourse.julialang.org/u/Gray_GALLANT)\
**Post date:** [June 7, 2025, 11:40pm UTC](https://discourse.julialang.org/t/symbolics-mtk-pipeline-lagrangians/129730/1 "2025-06-07T23:40:26Z")

</div>

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](https://global.discourse-cdn.com/julialang/original/3X/d/2/d274d3484811fd26c5d93165934784b69751caec.jpeg)

```julia
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?

---

<div class="post-metadata">

**Author:** ![contradict](https://avatars.discourse-cdn.com/v4/letter/c/ac91a4/32.png) [@contradict](https://discourse.julialang.org/u/contradict)\
**Post date:** [June 8, 2025, 8:59pm UTC](https://discourse.julialang.org/t/symbolics-mtk-pipeline-lagrangians/129730/2 "2025-06-08T20:59:53Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![Gray\_GALLANT](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gray_gallant/32/214253_2.png) [@Gray\_GALLANT](https://discourse.julialang.org/u/Gray_GALLANT)\
**Post date:** [June 9, 2025, 4:04am UTC](https://discourse.julialang.org/t/symbolics-mtk-pipeline-lagrangians/129730/3 "2025-06-09T04:04:46Z")

</div>

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

```julia
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](https://global.discourse-cdn.com/julialang/original/3X/d/b/db72f6e99b6e7ef1254ef3aa6f6abadabbc3c1dd.jpeg)

```julia

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.