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.
using AstrodynamicalModels, ModelingToolkit
model = R2BSystem()
J = calculate_jacobian(model)
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
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)
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:
12-element Vector{Any}:
disk₊b
disk₊G
halo₊b
halo₊m
disk₊m
bulge₊b
halo₊G
bulge₊m
bulge₊G