# Namespacing Symbols in \`ModelingToolkit.calculate\_jacobian\`

**URL:** <https://discourse.julialang.org/t/namespacing-symbols-in-modelingtoolkit-calculate-jacobian/123337>\
**Category:** Modelling & Simulations\
**Tags:** question, modelingtoolkit, namespaces\
**Created:** [December 2, 2024, 3:07am UTC](https://discourse.julialang.org/t/namespacing-symbols-in-modelingtoolkit-calculate-jacobian/123337 "2024-12-02T03:07:37Z")\
**Posts on this page:** 3\
**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:** [December 2, 2024, 3:07am UTC](https://discourse.julialang.org/t/namespacing-symbols-in-modelingtoolkit-calculate-jacobian/123337/1 "2024-12-02T03:07:37Z")

</div>

I (think I) need to retrieve the Jacobian of an `ODESystem`, but with namespaces for each symbol. For example, calling `calculate_jacobian` on a system below produces a matrix with symbols without namespaces.

```julia
using AstrodynamicalModels, ModelingToolkit

model = R2BSystem()

J = calculate_jacobian(model)

```

```julia
6×6 Matrix{Num}:
                                                                                                                                                                     0 … 1 0 0
                                                                                                                                                                     0 0 1 0
                                                                                                                                                                     0 0 0 1
 (-μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^3) - 3x(t)*((-x(t)*μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^6))*sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t))) 0 0 0
                                                        -3((-y(t)*μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^6))*x(t)*sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t))) 0 0 0
                                                        -3x(t)*((-z(t)*μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^6))*sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t))) … 0 0 0

```

For context, I’m trying to compose an `ODESystem` from components which specify scalar potentials in the integral form. The top-level `ODESystem` takes the gradient of each scalar potential variable with respect to `x`, `y`, and `z`. The line ` grad(sys) = calculate_jacobian...` is where the namespacing goes away. I have tried replacing the `calculate_jacobian` call with `Differential(sys.x)(sys.Φ)`, but this is not allowed because `ODESystem` instances can only have one `Differential` (with respect to the one independent variable).

> **Composition Example**
>
> ```julia
> using ModelingToolkit
> 
> function ScalarField(value, t, u, p; name, simplify = true, kwargs...)
> @variables Φ(t) [output = true]
> 
> return ODESystem(vcat(Φ ~ value), t; name = name, kwargs...)
> end
> 
> function PlummerPotential(; name = :PlummerPotential, kwargs...)
> @independent_variables t
> u = @variables x(t) [input = true] y(t) [input = true] z(t) [input = true]
> p = @parameters G m b
> 
> value = -G * m / sqrt(b^2 + x^2 + y^2 + z^2)
> return ScalarField(value, t, u, p; name = name, kwargs...)
> end
> 
> function MilkyWay(; name = :MilkyWay, kwargs...)
> @named disk = PlummerPotential()
> @named bulge = PlummerPotential()
> @named halo = PlummerPotential()
> 
> @independent_variables t
> @variables Φ(t) [output = true]
> D = Differential(t)
> u = @variables x(t) [input = true] y(t) [input = true] z(t) [input = true]
> du = @variables dx(t) [input = true] dy(t) [input = true] dz(t) [input = true]
> 
> aliases = [
> disk.x => x,
> disk.y => y,
> disk.z => z,
> bulge.x => x,
> bulge.y => y,
> bulge.z => z,
> halo.x => x,
> halo.y => y,
> halo.z => z
> ]
> 
> grad(sys) = calculate_jacobian(sys; simplify = true)[(begin + 1):end] # TODO remove manual indexing
> 
> eqs = vcat(
> Φ ~ disk.Φ + bulge.Φ + halo.Φ,
> D.(u) .~ du,
> D.(du) .~ -(grad(disk) .+ grad(bulge) .+ grad(halo)),
> [alias.first ~ alias.second for alias in aliases]
> )
> 
> return compose(
> ODESystem(
> eqs, t; name = name, defaults = Dict(aliases)),
> disk, bulge, halo
> )
> end
> 
> mw = MilkyWay()
> mws = structural_simplify(mw)
> parameters(mws)
> 
> ```
> 
> ```julia
> 12-element Vector{Any}:
> G
> b
> m
> disk₊b
> disk₊G
> halo₊b
> halo₊m
> disk₊m
> bulge₊b
> halo₊G
> bulge₊m
> bulge₊G
> 
> ```
> 
> The non-namespaces parameters should not be present, like below:
> 
> ```julia
> 12-element Vector{Any}:
> disk₊b
> disk₊G
> halo₊b
> halo₊m
> disk₊m
> bulge₊b
> halo₊G
> bulge₊m
> bulge₊G
> 
> ```

---

<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:** [December 3, 2024, 3:51am UTC](https://discourse.julialang.org/t/namespacing-symbols-in-modelingtoolkit-calculate-jacobian/123337/2 "2024-12-03T03:51:46Z")

</div>

I solved this by computing the gradient inside each subsystem, and storing the result in a new variable. After `structural_simplify` this variable is converted to an observable, so there is no hit to performance. Computing the gradient inside each subsystem preserves the namespaces for each term in the gradient. The commit hash for this change is [c407a55](https://github.com/cadojo/GalacticPotentials.jl/commit/c407a55f81d8df91ae324acf199cdd131d03d7c7), for anyone curious in the future.

---

<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:** [December 3, 2024, 4:41pm UTC](https://discourse.julialang.org/t/namespacing-symbols-in-modelingtoolkit-calculate-jacobian/123337/3 "2024-12-03T16:41:07Z")

</div>

You can also just do this on the flattened system.