N-body Problem Dynamics

Also known as NBP dynamics!

Overview

In an astrodynamical context, the N-body problem assumes $N$ celestial bodies which move with respect to some common origin. A body $i$ moves due to the cumulative gravity of every other body in the system. This problem is notoriously difficult because it cannot be solved analytically for $N\geq3$!

Examples

All NBSystem calls require the number of bodies to be specified as the first argument, like so. As always, use the stm argument at your leisure. Beware, though! using stm=true for N-body systems with more than 5 bodies may cause NBSystem to compute for a really, really long time! True story – v1.0.1 of this package had a version of these docs which tried to compute NBSystem(30; stm=true). That resulted in GitHub and JuliaHub failing the job on a timeout after several hours – at first I was surprised, until I realized that appending state transition matrix dynamics to a 30-body system results in a total of 32580 states!

If you're curious, the number of states of an N-body system with state transition matrix dynamics appended is equivalent to N*6 + (N*6)^2.

julia> model = NBSystem(2; stm=true)Model NBPWithSTM with 156 equations
Unknowns (156):
  (x(t))[1]
  (y(t))[1]
  (z(t))[1]
  (x(t))[2]
⋮
Parameters (3):
  G
  m[1]
  m[2]

Like other models, we can compute the Jacobian for these dynamics.

julia> using SparseArrays
julia> J = sparse(calculate_jacobian(NBSystem(4)))24×24 SparseArrays.SparseMatrixCSC{Num, Int64} with 156 stored entries:
⎡⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⎤
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⎥
⎢⣿⣿⣿⣿⣿⣿⠀⠀⠀⠀⠀⠀⎥
⎢⣿⣿⣿⣿⣿⣿⠀⠀⠀⠀⠀⠀⎥
⎣⣿⣿⣿⣿⣿⣿⠀⠀⠀⠀⠀⠀⎦

Finally, let's construct a Julia function which implements these dynamics!

julia> f = NBFunction(2)(::ODEFunction{true, SciMLBase.FullSpecialize, ModelingToolkit.var"#f#684"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x93f533b6, 0xd2605d4d, 0x05af5f1a, 0x13d0c335, 0x83c1de29), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xdbb37bc7, 0x7359defc, 0xb1f17fc6, 0x1177c116, 0xcf9de5ee), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.var"#51211#generated_observed#693"{Bool, ODESystem, Dict{Any, Any}, Vector{SymbolicUtils.BasicSymbolic{Real}}}, Nothing, ODESystem, Nothing, Nothing}) (generic function with 1 method)
julia> let u = randn(12), m = randn(2), G = rand(), t = 0
           f(u, [G, m...], t)
       end12-element Vector{Float64}:
 -0.8536588137291672
 -1.2917616554630802
  1.7808916382970057
  1.3093980584111857
  0.7222982383925539
 -1.666276641899703
  2.6783844804729723e-5
  9.21427659972393e-6
 -0.00010357795826493306
  0.0004366924640600646
  0.00015023254436397713
 -0.0016887685150056104