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:
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?