N-body Problem Dynamics

Also known as NBP dynamics!


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$!


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):
Parameters (3):

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