Fitting Potentials into `ModelingToolkit` Types

I’m writing a package that will be a library of common galactic dynamics potential models, e.g. the harmonic oscillator potential. The equation for this potential is \Phi = \frac{1}{2} \omega^2 x^2. I want to write this potential as some ModelingToolkit type to allow users to symbolically generate the gradient and hessian, and then hook into numerical simulations through the SciML ecosystem.

Is a NonlinearSystem appropriate here? Really, each potential is more like a function, but I believe there is benefit to be gained by fitting these models into the AbstractSystem interface. For example, models could be composed to form more complicated composite potential models that can be simplified auto-magically through ModelingToolkit.

Here’s my attempt at fitting this into a NonlinearSystem. Is this what you would do? Am I missing a better way to fit these models into the SciML ecosystem?

using ModelingToolkit
using LinearAlgebra

function HarmonicOscillatorPotential(N::Integer; name=:HarmonicOscillatorPotential)
    @parameters ω[1:N] [description = "The frequencies of the harmonic oscillator"]
    @variables x[1:N] [input = true, description = "The N-dimensional state of the oscillator"]
    @variables V [output = true, description = "The potential of the oscillator."]

    return NonlinearSystem(
        [V ~ sum(1 // 2 * ωᵢ^2 * xᵢ^2 for (ωᵢ, xᵢ) in zip(ω, x))],

This doesn’t look like something that requires NonlinearSystem since it’s not a nonlinear rootfinding problem, i.e. it’s not something of the form of:

Instead, this is just straight symbolic calculations. Note that ModelingToolkit is built on Symbolics.jl, so you can simply do the symbolic calculations. What you’re looking for are the symbolic differentiation functions, which are discussed in the main tutorial:

1 Like