# 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) end`

`12-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`