From 55be11ca50fb6787cbefee46a2b1b16606990f7d Mon Sep 17 00:00:00 2001 From: abhro <5664668+abhro@users.noreply.github.com> Date: Sun, 23 Aug 2026 23:50:01 +0600 Subject: [PATCH] Fix up outer solar system example MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Use backticks instead of dollar signs for TeX consistency - Use gradient from Symbolics.jl instead of relying on ModelingToolkit.jl to re-export it - Fix qᵢ/pᵢ description - Explicitly import the names being used, don't need namespace prefixes --- docs/Project.toml | 2 + docs/src/examples/outer_solar_system.md | 58 ++++++++++++------------- 2 files changed, 30 insertions(+), 30 deletions(-) diff --git a/docs/Project.toml b/docs/Project.toml index a25d5e292..e855fa4f3 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -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] @@ -187,5 +188,6 @@ StochasticDiffEqMilstein = "2" StochasticDiffEqROCK = "2" StochasticDiffEqRODE = "2" StochasticDiffEqWeak = "2" +Symbolics = "7.36" Sundials = "4.11.3, 5, 6.0" Unitful = "1" diff --git a/docs/src/examples/outer_solar_system.md b/docs/src/examples/outer_solar_system.md index 774752f46..49da2ad8e 100644 --- a/docs/src/examples/outer_solar_system.md +++ b/docs/src/examples/outer_solar_system.md @@ -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 @@ -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. @@ -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 ``` -