Skip to content

Fix up outer solar system example #898

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
ChrisRackauckas merged 1 commit into from
Aug 23, 2026
+30 −30
Merged

Fix up outer solar system example #898

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions docs/Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -91,6 +91,7 @@ StochasticDiffEqROCK = "db241ea8-0e6b-4abc-8f2d-1adff2294fd9"
StochasticDiffEqRODE = "49714585-0aa1-4f53-b1e6-a9b8c0d5e03f"
StochasticDiffEqWeak = "af2a2fcd-1c36-4cbe-a6d0-5afda784a085"
Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"

[sources]
Expand Down Expand Up @@ -187,5 +188,6 @@ StochasticDiffEqMilstein = "2"
StochasticDiffEqROCK = "2"
StochasticDiffEqRODE = "2"
StochasticDiffEqWeak = "2"
Symbolics = "7.36"
Sundials = "4.11.3, 5, 6.0"
Unitful = "1"
58 changes: 28 additions & 30 deletions docs/src/examples/outer_solar_system.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,39 +2,39 @@

## Data

The chosen units are masses relative to the sun, meaning the sun has mass $1$. We have taken $m_0 = 1.00000597682$ to take account of the inner planets. Distances are in astronomical units, times in earth days, and the gravitational constant is thus $G = 2.95912208286 \cdot 10^{-4}$.
The chosen units are masses relative to the sun, meaning the sun has mass ``1``. We have taken ``m_0 = 1.00000597682`` to take account of the inner planets. Distances are in astronomical units, times in earth days, and the gravitational constant is thus ``G = 2.95912208286 × 10^{-4}``.

| planet | mass | initial position | initial velocity |
|:------- |:--------------------------- |:--------------------------------------- |:-------------------------------------- |
| Jupiter | $m_1 = 0.000954786104043$ | [-3.5023653, -3.8169847, -1.5507963] | [0.00565429, -0.00412490, -0.00190589] |
| Saturn | $m_2 = 0.000285583733151$ | [9.0755314, -3.0458353, -1.6483708] | [0.00168318, 0.00483525, 0.00192462] |
| Uranus | $m_3 = 0.0000437273164546$ | [8.3101420, -16.2901086, -7.2521278] | [0.00354178, 0.00137102, 0.00055029] |
| Neptune | $m_4 = 0.0000517759138449$ | [11.4707666, -25.7294829, -10.8169456] | [0.00288930, 0.00114527, 0.00039677] |
| Pluto | $m_5 = 1/(1.3 \cdot 10^8 )$ | [-15.5387357, -25.2225594, -3.1902382] | [0.00276725, -0.00170702, -0.00136504] |
| planet | mass | initial position | initial velocity |
|:------- |:---------------------------- |:--------------------------------------- |:-------------------------------------- |
| Jupiter | ``m_1 = 0.000954786104043`` | [-3.5023653, -3.8169847, -1.5507963] | [0.00565429, -0.00412490, -0.00190589] |
| Saturn | ``m_2 = 0.000285583733151`` | [9.0755314, -3.0458353, -1.6483708] | [0.00168318, 0.00483525, 0.00192462] |
| Uranus | ``m_3 = 0.0000437273164546`` | [8.3101420, -16.2901086, -7.2521278] | [0.00354178, 0.00137102, 0.00055029] |
| Neptune | ``m_4 = 0.0000517759138449`` | [11.4707666, -25.7294829, -10.8169456] | [0.00288930, 0.00114527, 0.00039677] |
| Pluto | ``m_5 = 1/(1.3 × 10^8)`` | [-15.5387357, -25.2225594, -3.1902382] | [0.00276725, -0.00170702, -0.00136504] |

The data is taken from the book “Geometric Numerical Integration” by E. Hairer, C. Lubich and G. Wanner.

```@example outersolarsystem
import Plots, OrdinaryDiffEq as ODE
import ModelingToolkit as MTK
using ModelingToolkit: t_nounits as t, D_nounits as D, @mtkbuild, @variables
Plots.gr()
import OrdinaryDiffEq as ODE
using ModelingToolkit: System, t_nounits as t, D_nounits as D, @mtkbuild, @variables
using Symbolics: gradient
using Plots: plot, plot!

G = 2.95912208286e-4
M = [
1.00000597682,
0.000954786104043,
0.000285583733151,
0.0000437273164546,
0.0000517759138449,
9.54786104043e-4,
2.85583733151e-4,
4.37273164546e-5,
5.17759138449e-5,
1 / 1.3e8,
]
planets = ["Sun", "Jupiter", "Saturn", "Uranus", "Neptune", "Pluto"]

pos = [
0.0 -3.5023653 9.0755314 8.310142 11.4707666 -15.5387357
0.0 -3.5023653 9.0755314 8.310142 11.4707666 -15.5387357
0.0 -3.8169847 -3.0458353 -16.2901086 -25.7294829 -25.2225594
0.0 -1.5507963 -1.6483708 -7.2521278 -10.8169456 -3.1902382
0.0 -1.5507963 -1.6483708 -7.2521278 -10.8169456 -3.1902382
]
vel = [
0.0 0.00565429 0.00168318 0.00354178 0.0028893 0.00276725
Expand All @@ -47,10 +47,10 @@ tspan = (0.0, 200_000.0)
The N-body problem's Hamiltonian is

```math
H(p,q) = \frac{1}{2}\sum_{i=0}^{N}\frac{p_i^T p_i}{m_i} - G\sum_{i=1}^N \sum_{j=0}^{i-1}\frac{m_i m_j}{\left\lVert q_i - q_j \right\rVert}
H(p,q) = \frac{1}{2} ∑_{i=0}^N \frac{p_i^T p_i}{m_i} - G ∑_{i=1}^N ∑_{j=0}^{i-1} \frac{m_i m_j}{\left\| q_i - q_j \right\|}
```

where each ``p_i`` and ``q_i`` is a 3-dimensional vector describing the planet's position and momentum, respectively.
where each ``q_i`` and ``p_i`` is a 3-dimensional vector describing the planet's position and momentum, respectively.

Here, we want to solve for the motion of the five outer planets relative to the sun, namely, Jupiter, Saturn, Uranus, Neptune, and Pluto.

Expand All @@ -71,30 +71,28 @@ potential = -G *
`NBodyProblem` constructs a second order ODE problem under the hood. We know that a Hamiltonian system has the form of

```math
\dot{p} = -\frac{\partial H}{\partial q}, \quad \dot{q} = \frac{\partial H}{\partial p}
\dot{p} = -\frac{∂H}{∂q}, \quad \dot{q} = \frac{∂H}{∂p}
```

For an N-body system, we can simplify this as:

```math
\dot{p} = -\nabla V(q), \quad \dot{q} = M^{-1} p.
\dot{p} = -∇ V(q), \quad \dot{q} = M^{-1} p.
```

Thus, $\dot{q}$ is defined by the masses. We only need to define $\dot{p}$, and this is done internally by taking the gradient of $V$. Therefore, we only need to pass the potential function and the rest is taken care of.
Thus, ``\dot{q}`` is defined by the masses. We only need to define ``\dot{p}``, and this is done internally by taking the gradient of ``V``. Therefore, we only need to pass the potential function and the rest is taken care of.

```@example outersolarsystem
eqs = vec(@. D(D(u))) .~ .-MTK.gradient(potential, vec(u)) ./
repeat(M, inner = 3)
@mtkbuild sys = MTK.System(eqs, t)
eqs = vec(@. D(D(u))) .~ .-gradient(potential, vec(u)) ./ repeat(M, inner = 3)
@mtkbuild sys = System(eqs, t)
prob = ODE.ODEProblem(sys, [vec(u .=> pos); vec(D.(u) .=> vel)], tspan)
sol = ODE.solve(prob, ODE.Tsit5());
```

```@example outersolarsystem
plt = Plots.plot()
plt = plot(xlab = "x", ylab = "y", zlab = "z", title = "Outer solar system")
for i in 1:N
Plots.plot!(plt, sol, idxs = (u[:, i]...,), lab = planets[i])
plot!(plt, sol, idxs = (u[:, i]...,), lab = planets[i])
end
Plots.plot!(plt; xlab = "x", ylab = "y", zlab = "z", title = "Outer solar system")
plt
```

Loading