Fitting Potentials into `ModelingToolkit` Types

Specific Domains Modelling & Simulations
, ,
1

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))],
        collect(x),
        collect(ω);
        name=name
    )
end
2

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:

https://symbolics.juliasymbolics.org/stable/getting_started/

3 Likes
3

Would components be appropriate for defining potential fields, which can be summed together, and differentiated to get ODEs? Components and Connectors · ModelingToolkit.jl

4

Like modeling in the integral form? I think that could be done.

1 Like
5

I’d like to use the scalar equations for each potential because the scalar form is what is shown in textbooks, etc. I’m bootstrapping equations in Gala — a Python package for galactic dynamics — into forms compatible with SciML & MTK. The end goal for this project is to sum several of these scalar potentials to form a model of the Milky Way, then integrate orbits for massless particles using DiffEqPhysics.HamiltonianProblem. I’m a bit unclear on how I should define each potential as a component. Any guidance you, or others could provide would be much appreciated!

Here is what I have so far for a scalar potential due to a harmonic oscillator. You can see the full list at GalacticPotentials.jl; note though, that the committed potential functions are not components, but are instead a custom AbstractTimeDependenentSystem type called ScalarField. I’m trying to port those into the more standard component interface.

@mtkmodel HarmonicOscillatorPotential begin
    @variables begin
        Φ(t)
        x(t)
    end

    @parameters begin
        ω(t)
    end

    @equations begin
        Φ ~ (1 // 2) * ω * x
    end
end

The code above produces one equation with two unknowns. I’d like to do the following, but I get ERROR: LoadError: ArgumentError: expression should be of the form sys = foo(a, b).

@mtkbuild MilkyWayPotential begin
    @variables begin
        Φ(t)
    end

    @components begin
        disk = HarmonicOscillatorPotential()
        bulge = HarmonicOscillatorPotential()
        halo = HarmonicOscillatorPotential()
    end

    @equations begin
        Φ ~ disk.Φ + bulge.Φ + halo.Φ
    end
end

Is my @mtkmodel block correct? Will the fact that each Potential model (component) has two unknowns and one equation ultimately be an issue?

6

The acausal DSL is not the right DSL for this. You’ll want to develop a different macro which connects via summation like that.

7

Thanks for your replies & time. Understood that the @mtkbuild DSL is not right for this purpose. Would you still recommend using the Model type and @mtkmodel macro to define each individual potential?

8

You’d probably want to build a new macro that directly targets ODESystem.

1 Like