# Fitting Potentials into \`ModelingToolkit\` Types

**URL:** <https://discourse.julialang.org/t/fitting-potentials-into-modelingtoolkit-types/105213>\
**Category:** Modelling & Simulations\
**Tags:** question, plotting, modelingtoolkit\
**Created:** [October 20, 2023, 3:51am UTC](https://discourse.julialang.org/t/fitting-potentials-into-modelingtoolkit-types/105213 "2023-10-20T03:51:03Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![cadojo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cadojo/32/25328_2.png) [@cadojo](https://discourse.julialang.org/u/cadojo)\
**Post date:** [October 20, 2023, 3:51am UTC](https://discourse.julialang.org/t/fitting-potentials-into-modelingtoolkit-types/105213/1 "2023-10-20T03:51:03Z")

</div>

I’m writing a [package](https://github.com/cadojo/GalacticPotentials.jl) that will be a library of common galactic dynamics potential models, e.g. the [harmonic oscillator potential](http://gala.adrian.pw/en/latest/api/gala.potential.potential.HarmonicOscillatorPotential.html#gala.potential.potential.HarmonicOscillatorPotential. __call__ ). 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?

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [October 21, 2023, 2:19pm UTC](https://discourse.julialang.org/t/fitting-potentials-into-modelingtoolkit-types/105213/2 "2023-10-21T14:19:45Z")

</div>

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:

> **[Symbolic System Definition and Acceleration via ModelingToolkit ·...](https://docs.sciml.ai/NonlinearSolve/stable/tutorials/modelingtoolkit/)**
>
> Documentation for NonlinearSolve.jl.

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/](https://symbolics.juliasymbolics.org/stable/getting_started/)

---

<div class="post-metadata">

**Author:** ![cadojo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cadojo/32/25328_2.png) [@cadojo](https://discourse.julialang.org/u/cadojo)\
**Post date:** [June 19, 2024, 5:08am UTC](https://discourse.julialang.org/t/fitting-potentials-into-modelingtoolkit-types/105213/3 "2024-06-19T05:08:50Z")

</div>

Would components be appropriate for defining potential fields, which can be summed together, and differentiated to get ODEs? [Components and Connectors · ModelingToolkit.jl](https://docs.sciml.ai/ModelingToolkit/stable/basics/MTKModel_Connector/)

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 19, 2024, 11:37am UTC](https://discourse.julialang.org/t/fitting-potentials-into-modelingtoolkit-types/105213/4 "2024-06-19T11:37:08Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![cadojo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cadojo/32/25328_2.png) [@cadojo](https://discourse.julialang.org/u/cadojo)\
**Post date:** [June 22, 2024, 9:50pm UTC](https://discourse.julialang.org/t/fitting-potentials-into-modelingtoolkit-types/105213/5 "2024-06-22T21:50:31Z")

</div>

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](http://gala.adrian.pw/en/latest/) — 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](http://gala.adrian.pw/en/latest/tutorials/define-milky-way-model.html), then integrate orbits for massless particles using [`DiffEqPhysics.HamiltonianProblem`](https://docs.sciml.ai/DiffEqDocs/stable/types/dynamical_types/#Hamiltonian-Problems). 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`](https://github.com/cadojo/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.

```julia
@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)`.

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 22, 2024, 11:15pm UTC](https://discourse.julialang.org/t/fitting-potentials-into-modelingtoolkit-types/105213/6 "2024-06-22T23:15:42Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![cadojo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cadojo/32/25328_2.png) [@cadojo](https://discourse.julialang.org/u/cadojo)\
**Post date:** [June 23, 2024, 4:28am UTC](https://discourse.julialang.org/t/fitting-potentials-into-modelingtoolkit-types/105213/7 "2024-06-23T04:28:28Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 23, 2024, 6:55am UTC](https://discourse.julialang.org/t/fitting-potentials-into-modelingtoolkit-types/105213/8 "2024-06-23T06:55:31Z")

</div>

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