Update solar system

.buildkite/run_tutorial.yml

@@ -1,6 +1,7 @@
 agents:
   queue: "juliaecosystem"
   sandbox.jl: true
+  arch: "x86_64"
 
 # This is a pipeline that weaves a tutorial, then uploads the resultant
 # .PDF and other reports as (buildkite, not Julia) artifacts.  The `coppermind`
@@ -18,7 +19,7 @@ steps:
             - BUILDKITE_S3_SECRET_ACCESS_KEY="U2FsdGVkX1+SPF81nkK7KQ64DsafSl0qq2iG7BsQs1xlTYEtZV3MqQl3l/NWaiocaEywZZFbAB5zpnKPD0xHTQ=="
             - BUILDKITE_S3_DEFAULT_REGION="U2FsdGVkX1/cORlxhXcxhja2JkqC0f8RmaGYxvGBbEg="
       - JuliaCI/julia#v1:
-          version: 1.6
+          version: 1.8
       - staticfloat/sandbox:
           rootfs_url: "https://jc-rootfs-images.s3.amazonaws.com/aws_uploader-2021-11-12.x86_64.tar.gz"
           rootfs_treehash: "986217e5b36efd3b3b91ed90df8e36d628cf543f"
@@ -63,7 +64,7 @@ steps:
           files:
             - .buildkite/ssh_deploy.key
       - JuliaCI/julia#v1:
-          version: 1.7
+          version: 1.8
       - staticfloat/sandbox:
           rootfs_url: "https://jc-rootfs-images.s3.amazonaws.com/aws_uploader-2021-11-12.x86_64.tar.gz"
           rootfs_treehash: "986217e5b36efd3b3b91ed90df8e36d628cf543f"

.buildkite/test_sciml.yml

@@ -1,11 +1,12 @@
 agents:
   queue: "juliaecosystem"
+  arch: "x86_64"
 
 steps:
-  - label: ":julia: Run tests on 1.7"
+  - label: ":julia: Run tests on 1.8"
     plugins:
       - JuliaCI/julia#v1:
-          version: 1.7
+          version: 1.8
       - JuliaCI/julia-test#v1:
     timeout_in_minutes: 20
     artifact_paths:

README.md

@@ -52,6 +52,7 @@ SciMLTutorials.open_notebooks()
   - [Tutorial on using spatial SSAs in DiffEqJump](http://tutorials.sciml.ai/html/jumps/spatial.html)
   - [Kepler Problem Orbit](http://tutorials.sciml.ai/html/models/05-kepler_problem.html)
   - [Spiking Neural Systems](http://tutorials.sciml.ai/html/models/08-spiking_neural_systems.html)
+  - [Outer Solar System](http://tutorials.sciml.ai/html/models/05-outer_solar_system.html)
 - Advanced ODE Features
   - [Feagin's Order 10, 12, and 14 Methods](http://tutorials.sciml.ai/html/ode_extras/01-feagin.html)
   - [Finding Maxima and Minima of DiffEq Solutions](http://tutorials.sciml.ai/html/ode_extras/02-ode_minmax.html)

src/SciMLTutorials.jl

@@ -8,12 +8,12 @@ latexfile = joinpath(@__DIR__, "..", "templates", "julia_tex.tpl")
 default_builds = (:script,:github)
 
 function weave_file(folder,file,build_list=default_builds)
-  target = joinpath(repo_directory, "tutorials", folder, file)
+  target = joinpath(folder, file)
   @info("Weaving $(target)")
 
-  if isfile(joinpath(repo_directory, "tutorials", folder, "Project.toml"))
+  if isfile(joinpath(folder, "Project.toml"))
     @info("Instantiating", folder)
-    Pkg.activate(joinpath(repo_directory,"tutorials", folder))
+    Pkg.activate(joinpath(folder))
     Pkg.instantiate()
     Pkg.build()
 
@@ -66,7 +66,7 @@ function weave_all(build_list=default_builds)
 end
 
 function weave_folder(folder,build_list=default_builds)
-  for file in readdir(joinpath(repo_directory,"tutorials",folder))
+  for file in readdir(joinpath(folder))
     # Skip non-`.jmd` files
     if !endswith(file, ".jmd")
       continue

tutorials/models/05-outer_solar_system.jmd

@@ -0,0 +1,84 @@
+---
+title: The Outer Solar System
+author: Yingbo Ma, Chris Rackauckas
+---
+
+## Data
+
+
+The chosen units are: masses relative to the sun, so that 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}$.
+
+| planet | mass | initial position | initial velocity |
+| --- | --- | --- | --- |
+| Jupiter | $m_1 = 0.000954786104043$ | <ul><li>-3.5023653</li><li>-3.8169847</li><li>-1.5507963</li></ul> | <ul><li>0.00565429</li><li>-0.00412490</li><li>-0.00190589</li></ul>
+| Saturn | $m_2 = 0.000285583733151$ | <ul><li>9.0755314</li><li>-3.0458353</li><li>-1.6483708</li></ul> | <ul><li>0.00168318</li><li>0.00483525</li><li>0.00192462</li></ul>
+| Uranus | $m_3 = 0.0000437273164546$ | <ul><li>8.3101420</li><li>-16.2901086</li><li>-7.2521278</li></ul> | <ul><li>0.00354178</li><li>0.00137102</li><li>0.00055029</li></ul>
+| Neptune | $m_4 = 0.0000517759138449$ | <ul><li>11.4707666</li><li>-25.7294829</li><li>-10.8169456</li></ul> | <ul><li>0.00288930</li><li>0.00114527</li><li>0.00039677</li></ul>
+| Pluto | $ m_5 = 1/(1.3 \cdot 10^8 )$ | <ul><li>-15.5387357</li><li>-25.2225594</li><li>-3.1902382</li></ul> | <ul><li>0.00276725</li><li>-0.00170702</li><li>-0.00136504</li></ul>
+
+The data is taken from the book "Geometric Numerical Integration" by E. Hairer, C. Lubich and G. Wanner.
+
+```julia
+using Plots, OrdinaryDiffEq, ModelingToolkit
+gr()
+
+G = 2.95912208286e-4
+M = [1.00000597682, 0.000954786104043, 0.000285583733151, 0.0000437273164546, 0.0000517759138449, 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.8169847  -3.0458353  -16.2901086  -25.7294829  -25.2225594
+        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
+       0.0  -0.0041249   0.00483525  0.00137102  0.00114527  -0.00170702
+       0.0  -0.00190589  0.00192462  0.00055029  0.00039677  -0.00136504]
+tspan = (0.0, 200_000.0)
+```
+
+The N-body problem's Hamiltonian is
+
+$$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}$$
+
+Here, we want to solve for the motion of the five outer planets relative to the sun, namely, Jupiter, Saturn, Uranus, Neptune and Pluto.
+
+```julia
+const ∑ = sum
+const N = 6
+@variables t u(t)[1:3, 1:N]
+u = collect(u)
+D = Differential(t)
+potential = -G*∑(i->∑(j->(M[i]*M[j])/√(∑(k->(u[k, i]-u[k, j])^2, 1:3)), 1:i-1), 2:N)
+```
+
+## Hamiltonian System
+
+`NBodyProblem` constructs a second order ODE problem under the hood. We know that a Hamiltonian system has the form of
+
+$$\dot{p} = -H_{q}(p,q)\quad \dot{q}=H_{p}(p,q)$$
+
+For an N-body system, we can symplify this as:
+
+$$\dot{p} = -\nabla{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.
+
+```julia
+eqs = vec(@. D(D(u))) .~ .- ModelingToolkit.gradient(potential, vec(u)) ./ repeat(M, inner=3)
+@named sys = ODESystem(eqs, t)
+ss = structural_simplify(sys)
+prob = ODEProblem(ss, [vec(u .=> pos); vec(D.(u) .=> vel)], tspan)
+sol = solve(prob, Tsit5());
+```
+
+```julia
+plt = plot()
+for i in 1:N
+    plot!(plt, sol, idxs=(u[:, i]...,), lab = planets[i])
+end
+plot!(plt; xlab = "x", ylab = "y", zlab = "z", title = "Outer solar system")
+```
+
+```julia, echo = false, skip="notebook"
+using SciMLTutorials
+SciMLTutorials.tutorial_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])
+```